# -*- coding: utf-8 -*-
"""
/***************************************************************************
 FragScape
                                 A QGIS plugin
 Computes ecological continuities based on environments permeability
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                             -------------------
        begin                : 2018-04-12
        git sha              : $Format:%H$
        copyright            : (C) 2018 by IRSTEA
        email                : mathieu.chailloux@irstea.fr
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""

from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import QgsProcessing, QgsProcessingAlgorithm, QgsProcessingException
from qgis.core import (Qgis,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterExpression,
                       QgsProcessingParameterFeatureSink,
                       QgsProcessingProvider,
                       QgsProcessingParameterMultipleLayers,
                       QgsProcessingUtils,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterMatrix,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterCrs,
                       QgsProcessingParameterVectorDestination,
                       QgsProcessingParameterFile,
                       QgsProcessingParameterString,
                       QgsProcessingParameterEnum,
                       QgsProperty,
                       QgsWkbTypes,
                       QgsCoordinateReferenceSystem,
                       QgsProcessingMultiStepFeedback)
from qgis.core import QgsField, QgsFields, QgsFeature, QgsFeatureSink

from processing.algs.gdal.rasterize import rasterize
import xml.etree.ElementTree as ET

from ..qgis_lib_mc import utils, qgsTreatments, qgsUtils, feedbacks
from ..steps import params
            
NB_DIGITS = 5
            
class FragScapeVectorAlgorithm(qgsUtils.BaseProcessingAlgorithm):
    
    def group(self):
        return "Vector"
    
    def groupId(self):
        return "fsVect"
        
    def displayName(self):
        assert(False)
        
    def shortHelpString(self):
        assert(False)
            

class RasterizeFixAllTouch(qgsUtils.BaseProcessingAlgorithm,rasterize):

    ALG_NAME = 'rasterizefixalltouch'
        
    def displayName(self):
        return self.tr('Rasterize (with ALL_TOUCH fix)')
        
    def group(self):
        return "Auxiliary algorithms"
        
    def groupId(self):
        return 'aux'
        
    def shortHelpString(self):
        return self.tr('Wrapper for gdal:rasterize algorithm allowing to use ALL_TOUCH option (every pixel touching input geometry are rasterized).')

    def initAlgorithm(self, config=None):
        super().initAlgorithm(config)
        self.addParameter(
            QgsProcessingParameterBoolean(
                self.ALL_TOUCH,
                description = 'ALL_TOUCH option',
                defaultValue=False,
                optional=True))
    
# Apply rasterization on field 'field' of vector layer 'in_path'.
# Output raster layer in 'out_path'.
# Resolution set to 25 if not given.
# Extent can be given through 'extent_path'. If not, it is extracted from input layer.
# Output raster layer is loaded in QGIS if 'load_flag' is True.
def applyRasterizationFixAllTouch(in_path,out_path,extent,resolution,
                   field=None,burn_val=None,out_type=Qgis.Float32,
                   nodata_val=qgsTreatments.nodata_val,all_touch=False,overwrite=False,
                   context=None,feedback=None):
    TYPES = ['Byte', 'Int16', 'UInt16', 'UInt32', 'Int32', 'Float32',
         'Float64', 'CInt16', 'CInt32', 'CFloat32', 'CFloat64']
    if overwrite:
        qgsUtils.removeRaster(out_path)
    extra_param_name = 'EXTRA'
    if hasattr(rasterize,extra_param_name):
        res = qgsTreatments.applyRasterization(in_path,out_path,extent,resolution,
                field,burn_val,out_type,nodata_val,all_touch,overwrite,
                context,feedback)
    else:
        parameters = { 'ALL_TOUCH' : True,
                   'BURN' : burn_val,
                   'DATA_TYPE' : out_type,
                   'EXTENT' : extent,
                   'FIELD' : field,
                   'HEIGHT' : resolution,
                   'INPUT' : in_path,
                   'NODATA' : nodata_val,
                   'OUTPUT' : out_path,
                   'UNITS' : 1, 
                   'WIDTH' : resolution }
        res = qgsTreatments.applyProcessingAlg("FragScape","rasterizefixalltouch",parameters,context,feedback)
    return res

class PrepareLanduseAlgorithm(FragScapeVectorAlgorithm):

    ALG_NAME = "prepareLanduse"

    INPUT = "INPUT"
    CLIP_LAYER = "CLIP_LAYER"
    SELECT_EXPR = "SELECT_EXPR"
    OUTPUT = "OUTPUT"
        
    def displayName(self):
        return self.tr("Prepare land cover data")
        
    def shortHelpString(self):
        return self.tr("This algorithms prepares land cover data by applying selection (from expression) and dissolving geometries")

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr("Input layer"),
                [QgsProcessing.TypeVectorAnyGeometry]))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.CLIP_LAYER,
                description=self.tr("Clip layer"),
                types=[QgsProcessing.TypeVectorPolygon],
                optional=True))
        self.addParameter(
            QgsProcessingParameterExpression(
                self.SELECT_EXPR,
                self.tr("Selection expression"),
                "",
                self.INPUT))
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr("Output layer")))
                
    def processAlgorithm(self,parameters,context,feedback):
        # Dummy function to enable running an alg inside an alg
        # def no_post_process(alg, context, feedback):
            # pass
        input = self.parameterAsVectorLayer(parameters,self.INPUT,context)
        feedback.pushDebugInfo("input = " + str(input))
        if input is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))
        qgsUtils.normalizeEncoding(input)
        feedback.pushDebugInfo("input ok")
        clip_layer = self.parameterAsVectorLayer(parameters,self.CLIP_LAYER,context)
        expr = self.parameterAsExpression(parameters,self.SELECT_EXPR,context)
        if clip_layer is None:
            clipped = input
        else:
            clipped_path = params.mkTmpLayerPath('landuseClipped.gpkg')
            qgsTreatments.applyVectorClip(input,clip_layer,clipped_path,context,feedback)
            clipped = qgsUtils.loadVectorLayer(clipped_path)
            feedback.pushDebugInfo("clipped  = " + str(clipped))
        selected_path = params.mkTmpLayerPath('landuseSelection.gpkg')
        qgsTreatments.selectGeomByExpression(clipped,expr,selected_path,'landuseSelection')
        # selected = qgsUtils.loadVectorLayer(selected_path)
        # selected = qgsTreatments.extractByExpression(
           # clipped,expr,'memory:',
           # context=context,feedback=feedback)
        feedback.pushDebugInfo("selected = " + str(selected_path))
        output = parameters[self.OUTPUT]
        dissolved = qgsTreatments.dissolveLayer(selected_path,output,context=context,feedback=feedback)
        dissolved = None
        return {self.OUTPUT : dissolved}
        
        
class PrepareFragmentationAlgorithm(FragScapeVectorAlgorithm):

    ALG_NAME = "prepareFragm"

    INPUT = "INPUT"
    CLIP_LAYER = "CLIP_LAYER"
    SELECT_EXPR = "SELECT_EXPR"
    BUFFER = "BUFFER_EXPR"
    NAME = "NAME"
    OUTPUT = "OUTPUT"
        
    def displayName(self):
        return self.tr("Prepare vector data")
        
    def shortHelpString(self):
        return self.tr("This algorithm prepares a vector layer by applying clip, selection and buffer")

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                description=self.tr("Input layer"),
                types=[QgsProcessing.TypeVectorLine]))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.CLIP_LAYER,
                description=self.tr("Clip layer"),
                types=[QgsProcessing.TypeVectorPolygon],
                optional=True))
        self.addParameter(
            QgsProcessingParameterExpression(
                self.SELECT_EXPR,
                description=self.tr("Selection expression"),
                parentLayerParameterName=self.INPUT,
                optional=True))
        self.addParameter(
            QgsProcessingParameterExpression(
                self.BUFFER,
                description=self.tr("Buffer expression"),
                parentLayerParameterName=self.INPUT,
                optional=True))
        self.addParameter(
            QgsProcessingParameterString(
                self.NAME,
                description=self.tr("Identifier"),
                optional=True))
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                description=self.tr("Output layer")))
                
    def processAlgorithm(self,parameters,context,feedback):
        # Parameters
        feedback.pushDebugInfo("parameters = " + str(parameters))
        input = self.parameterAsVectorLayer(parameters,self.INPUT,context)
        if input is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT))
        qgsUtils.normalizeEncoding(input)
        clip = self.parameterAsVectorLayer(parameters,self.CLIP_LAYER,context)
        clip_flag = (clip is None)
        select_expr = self.parameterAsExpression(parameters,self.SELECT_EXPR,context)
        feedback.pushDebugInfo("select_expr : " + str(select_expr))
        buffer_expr = self.parameterAsExpression(parameters,self.BUFFER,context)
        name = self.parameterAsString(parameters,self.NAME,context)
        if not name:
            name = 'fragm'
        feedback.pushDebugInfo("buffer_expr : " + str(buffer_expr))
        if buffer_expr == "" and input.geometryType() != QgsWkbTypes.PolygonGeometry:
           feedback.pushDebugInfo("Empty buffer with non-polygon layer")
        output = parameters[self.OUTPUT]
        if clip is None:
            clipped = input
        else:
            clipped_path = params.mkTmpLayerPath(name + 'Clipped.gpkg')
            qgsTreatments.applyVectorClip(input,clip,clipped_path,context,feedback)
            clipped = qgsUtils.loadVectorLayer(clipped_path)
        if select_expr == "":
            selected = clipped
        else:
            selected_path = params.mkTmpLayerPath(name + 'Selected.gpkg')
            qgsTreatments.selectGeomByExpression(clipped,select_expr,selected_path,name)
            selected = selected_path
        if buffer_expr == "":
            buffered = selected
        else:
            buffer_expr_prep = QgsProperty.fromExpression(buffer_expr)
            buffered = qgsTreatments.applyBufferFromExpr(selected,buffer_expr_prep,output,context,feedback)
        if buffered == input:
            buffered = qgsUtils.pathOfLayer(buffered)
        return {self.OUTPUT : buffered}
        

        
class ApplyFragmentationAlgorithm(FragScapeVectorAlgorithm):

    ALG_NAME = "applyFragm"

    LANDUSE = "LANDUSE"
    FRAGMENTATION = "FRAGMENTATION"
    CRS = "CRS"
    OUTPUT = "OUTPUT"
        
    def displayName(self):
        return self.tr("Integrates vector data to land cover")
        
    def shortHelpString(self):
        return self.tr("This algorithm builds a layer of patches from a land cover layer and fragmentation layers. Overlaying geometries are removed and remaining ones are cast to single geometry type.")

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.LANDUSE,
                self.tr("Land cover layer"),
                [QgsProcessing.TypeVectorPolygon]))
        self.addParameter(
            QgsProcessingParameterMultipleLayers(
                self.FRAGMENTATION,
                self.tr("Fragmentation layers"),
                QgsProcessing.TypeVectorPolygon))
        self.addParameter(
            QgsProcessingParameterCrs(
                self.CRS,
                description=self.tr("Output CRS"),
                defaultValue=params.defaultCrs))
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr("Output layer")))
                
    def processAlgorithm(self,parameters,context,feedback):
        # Parameters
        landuse = self.parameterAsVectorLayer(parameters,self.LANDUSE,context)
        if landuse is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.LANDUSE))
        qgsUtils.normalizeEncoding(landuse)
        fragm_layers = self.parameterAsLayerList(parameters,self.FRAGMENTATION,context)
        output = self.parameterAsOutputLayer(parameters,self.OUTPUT,context)
        crs = self.parameterAsCrs(parameters,self.CRS,context)
        # output = parameters[self.OUTPUT]
        # Merge fragmentation layers
        fragm_path = params.mkTmpLayerPath("fragm.gpkg")
        fragm_layer = qgsTreatments.mergeVectorLayers(fragm_layers,crs,fragm_path)
        feedback.pushDebugInfo("fragm_layer = " + str(fragm_layer))
        if fragm_layer is None:
            raise QgsProcessingException("Fragmentation layers merge failed")
        # Apply difference
        diff_path = params.mkTmpLayerPath("diff.gpkg")
        diff_layer = qgsTreatments.applyDifference(
            landuse,fragm_layer,diff_path,
            context=context,feedback=feedback)
        if fragm_layer is None:
            raise QgsProcessingException("Difference landuse/fragmentation failed")
        # Multi to single part
        singleGeomLayer = qgsTreatments.multiToSingleGeom(
            diff_layer,output,
            context=context,feedback=feedback)
        if fragm_layer is None:
            raise QgsProcessingException("Multi to single part failed")
        return {self.OUTPUT : singleGeomLayer}
        
 
class MeffAlgUtils:

    NB_DIGITS = 5
    
    INPUT = "INPUT"
    REPORTING = "REPORTING"
    CLASS = "CLASS"
    CRS = "CRS"
    INCLUDE_CBC = "INCLUDE_CBC"
    UNIT = "UNIT"
    OUTPUT = "OUTPUT"
    OUTPUT_VAL = "OUTPUT_VAL"
    
    SUM_AI = "sum_ai"
    SUM_AI_SQ = "sum_ai_sq"
    SUM_AI_SQ_CBC = "sum_ai_sq_cbc"
    DIVISOR = "divisor"
    
    SUFFIX = "suffix"
    CLIP_FLAG = "clip_flag"

    # Output layer fields
    ID = "fid"
    NB_PATCHES = "patches"
    REPORT_AREA = "At"
    INTERSECTING_AREA = "sum_Ai"
    # Main measures
    MESH_SIZE = "MSIZ"
    CBC_MESH_SIZE = "CBC_MSIZ"
    DIVI = "DIVI"
    SPLITTING_INDEX = "SPLI"
    # Auxiliary measures
    COHERENCE = "COH"
    SPLITTING_DENSITY = "SDEN"
    NET_PRODUCT = "NPRO"
    CBC_NET_PRODUCT = "CBC_NPRO"
    
    FIELDS_V1 = {
        ID : "fid",
        NB_PATCHES : "nb_patches",
        REPORT_AREA : "report_area",
        INTERSECTING_AREA : "intersecting_area",
        MESH_SIZE : "effective_mesh_size",
        CBC_MESH_SIZE : "CBC_effective_mesh_size",
        DIVI : "landscape_division",
        SPLITTING_INDEX : "splitting_index",
        COHERENCE : "coherence",
        SPLITTING_DENSITY : "splitting_density",
        NET_PRODUCT : "net_product",
        CBC_NET_PRODUCT : "CBC_net_product" }
    FIELDS_V1_SHP = {
        ID : "fid",
        NB_PATCHES : "nb_patches",
        REPORT_AREA : "report_are",
        INTERSECTING_AREA : "intersecti",
        MESH_SIZE : "effective_",
        CBC_MESH_SIZE : "CBC_effect",
        DIVI : "landscape_",
        SPLITTING_INDEX : "splitting_",
        COHERENCE : "coherence",
        SPLITTING_DENSITY : "splittin_1",
        NET_PRODUCT : "net_produc",
        CBC_NET_PRODUCT : "CBC_net_pr" }
    FIELDS_V2 = {
        ID : "fid",
        NB_PATCHES : "patches",
        REPORT_AREA : "At",
        INTERSECTING_AREA : "sum_Ai",
        MESH_SIZE : "MSIZ",
        CBC_MESH_SIZE : "CBC_MSIZ",
        DIVI : "DIVI",
        SPLITTING_INDEX : "SPLI",
        COHERENCE : "COH",
        SPLITTING_DENSITY : "SDEN",
        NET_PRODUCT : "NPRO",
        CBC_NET_PRODUCT : "CBC_NPRO" }
    
    UNIT_DIVISOR = [1, 100, 10000, 1000000]
    
    DEFAULT_CRS = QgsCoordinateReferenceSystem("epsg:2154")
    
    def getFieldSet(self,fields):
        if "report_area" in fields:
            return self.FIELDS_V1
        elif "report_are" in fields:
            return self.FIELDS_V1_SHP
        elif "At" in fields:
            return self.FIELDS_V2
        else:
            raise QgsProcessingException("Unexpected fields : " + str(fields))
    
    def getUnitOptions(self):
        return [self.tr("m² (square meters)"),
            self.tr("dm² (square decimeters / ares)"),
            self.tr("hm² (square hectometers / hectares)"),
            self.tr("km² (square kilometers)")]
            
    def mkReportFields(self,include_cbc=False):
        report_id_field = QgsField(self.ID, QVariant.Int)
        mesh_size_field = QgsField(self.MESH_SIZE, QVariant.Double)
        nb_patches_field = QgsField(self.NB_PATCHES, QVariant.Int)
        report_area_field = QgsField(self.REPORT_AREA, QVariant.Double)
        intersecting_area_field = QgsField(self.INTERSECTING_AREA, QVariant.Double)
        div_field = QgsField(self.DIVI, QVariant.Double)
        split_index_field = QgsField(self.SPLITTING_INDEX, QVariant.Double)
        coherence_field = QgsField(self.COHERENCE, QVariant.Double)
        split_density_field = QgsField(self.SPLITTING_DENSITY, QVariant.Double)
        net_product_field = QgsField(self.NET_PRODUCT, QVariant.Double)
        unit_divisor_field = QgsField(self.DIVISOR, QVariant.Int)
        if include_cbc:
            cbc_mesh_size_field = QgsField(self.CBC_MESH_SIZE, QVariant.Double)
            cbc_net_product_field = QgsField(self.CBC_NET_PRODUCT, QVariant.Double)
        output_fields = QgsFields()
        # output_fields.append(report_id_field)
        if include_cbc:
            output_fields.append(cbc_mesh_size_field)
        output_fields.append(mesh_size_field)
        output_fields.append(nb_patches_field)
        output_fields.append(report_area_field)
        output_fields.append(intersecting_area_field)
        output_fields.append(div_field)
        output_fields.append(split_index_field)
        output_fields.append(coherence_field)
        output_fields.append(split_density_field)
        output_fields.append(net_product_field)
        if include_cbc:
            output_fields.append(cbc_net_product_field)
        output_fields.append(unit_divisor_field)
        return output_fields
                
    def mkResFeat(self,include_cbc):
        if not self.report_layer or self.report_layer.featureCount() == 0:
            raise QgsProcessingException("Invalid reporting layer")
        for f in self.report_layer.getFeatures():
            report_feat = f
        output_fields = self.mkReportFields(include_cbc)
        res_feat = QgsFeature(output_fields)
        res_feat.setGeometry(report_feat.geometry())
        # res_feat[self.ID] = report_feat.id()
        return res_feat
        
    def fillResFeat(self,res_feat,res_dict):
        divisor = float(res_dict[self.DIVISOR])
        report_area = float(res_dict[self.REPORT_AREA]) / divisor
        report_area_sq = report_area * report_area
        sum_ai = float(res_dict[self.SUM_AI])  / divisor
        sum_ai_sq = float(res_dict[self.SUM_AI_SQ]) / (divisor * divisor)
        utils.debug("divisor = " + str(divisor))
        utils.debug("sum_ai = " + str(sum_ai))
        utils.debug("sum_ai_sq = " + str(sum_ai_sq))
        utils.debug("report_area = " + str(report_area))
        utils.debug("report_area_sq = " + str(report_area_sq))
        res_feat[self.NB_PATCHES] = res_dict[self.NB_PATCHES]
        # Metrics
        res_feat[self.NET_PRODUCT] = round(sum_ai_sq,NB_DIGITS)
        res_feat[self.REPORT_AREA] = report_area
        res_feat[self.INTERSECTING_AREA] = sum_ai
        res_feat[self.COHERENCE] = sum_ai_sq / report_area_sq if report_area_sq > 0 else 0
        res_feat[self.SPLITTING_DENSITY] = report_area / sum_ai_sq if sum_ai_sq > 0 else 0
        res_feat[self.MESH_SIZE] = round(sum_ai_sq / report_area, NB_DIGITS) if report_area > 0 else 0
        res_feat[self.SPLITTING_INDEX] = report_area_sq / sum_ai_sq if sum_ai_sq > 0 else 0
        res_feat[self.DIVI] = 1 - res_feat[self.COHERENCE]
        res_feat[self.DIVISOR] = divisor
        # CBC Metrics
        if self.SUM_AI_SQ_CBC in res_dict:
            sum_ai_sq_cbc = float(res_dict[self.SUM_AI_SQ_CBC]) / (divisor * divisor)
            res_feat[self.CBC_NET_PRODUCT] = round(sum_ai_sq_cbc,NB_DIGITS)
            res_feat[self.CBC_MESH_SIZE] = round(sum_ai_sq_cbc / report_area,NB_DIGITS) if report_area > 0 else 0
            
        
    def mkResSink(self,parameters,res_feat,context,include_cbc=False):
        report_fields = self.mkReportFields(include_cbc)
        wkb_type = self.report_layer.wkbType()
        report_crs = self.report_layer.sourceCrs()
        sink, dest_id = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            report_fields,
            wkb_type,
            report_crs)
        sink.addFeature(res_feat)
        return dest_id
        
    def mkOutputs(self,parameters,res_dict,context):
        if self.report_layer:
            nb_feats = self.report_layer.featureCount()
            if nb_feats != 1:
                raise QgsProcessingException("Report layer has "
                    + str(nb_feats) + " features but only 1 was expected")
            include_cbc = self.SUM_AI_SQ_CBC in res_dict
            utils.debug("include_cbc = " + str(include_cbc))
            res_feat = self.mkResFeat(include_cbc)
            self.fillResFeat(res_feat,res_dict)
            dest_id = self.mkResSink(parameters,res_feat,context,include_cbc)
            res_layer = dest_id
            if include_cbc:
                res_val = res_feat[self.CBC_MESH_SIZE]# / res_dict[self.DIVISOR]
            else:
                res_val = res_feat[self.MESH_SIZE]# / res_dict[self.DIVISOR]
        else:
            res_layer = None
            if res_dict[self.REPORT_AREA] > 0:
                res_val = (res_dict[self.SUM_AI_SQ] / res_dict[self.REPORT_AREA]) / res_dict[self.DIVISOR]
            else:
                res_val = 0
        res_val = round(res_val, self.NB_DIGITS)
        return (res_layer, res_val)


class ResultsDiffAlgorithm(MeffAlgUtils,qgsUtils.BaseProcessingAlgorithm):

    ALG_NAME = 'diffResults'
    
    LAYER_A = "LAYER_A"
    LAYER_B = "LAYER_B"
    
    PREFIX = "B_"
    VARIATION = 'variation'
        
    def displayName(self):
        return self.tr('Compare results layer')
        
    def shortHelpString(self):
        msg = "Compare 2 results layers produced by FragScape (step 4)"
        msg += " by applying difference between indicators values (layer_b - layer_a)"
        return self.tr(msg)
        
    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.LAYER_A,
                self.tr("Layer A"),
                [QgsProcessing.TypeVectorPolygon]))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.LAYER_B,
                self.tr("Layer B"),
                [QgsProcessing.TypeVectorPolygon]))
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr("Output layer")))
                
    def mkDiff(self,a_val,b_val,a_divi,b_divi,feedback,sq_factor=False):
        if a_divi == b_divi:
            feedback.pushDebugInfo("a_divi == b_divi")
            a_factor, b_factor, divi = 1, 1, a_divi
        elif a_divi < b_divi:
            a_factor, b_factor, divi = (b_divi / a_divi), 1, b_divi
        else:
            a_factor, b_factor, divi = 1, (a_divi / b_divi), a_divi
        if sq_factor:
            a_factor = a_factor * a_factor
            b_factor = b_factor * b_factor
        a_val_round = round(a_val / a_factor, self.NB_DIGITS)
        b_val_round = round(b_val / b_factor, self.NB_DIGITS)
        feedback.pushDebugInfo("a_factor = " + str(a_factor))
        feedback.pushDebugInfo("b_factor = " + str(b_factor))
        feedback.pushDebugInfo("a_val_round = " + str(a_val_round))
        feedback.pushDebugInfo("b_val_round = " + str(b_val_round))
        diff_val = b_val_round - a_val_round
        sum_val = a_val_round + b_val_round
        return (diff_val, sum_val, divi)
                
    def processAlgorithm(self,parameters,context,feedback):
        layer_a = self.parameterAsVectorLayer(parameters,self.LAYER_A,context)
        layer_b = self.parameterAsVectorLayer(parameters,self.LAYER_B,context)
        if layer_a is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.LAYER_A))
        if layer_b is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.LAYER_B))
        a_crs, b_crs = layer_a.sourceCrs(), layer_b.sourceCrs()
        if a_crs.authid() != b_crs.authid():
            raise QgsProcessingException("Different CRS for layers (" +
                str(a_crs) + " vs " + str(b_crs))
        nb_feats_a, nb_feats_b = layer_a.featureCount(), layer_b.featureCount()
        if nb_feats_a != nb_feats_b:
            raise QgsProcessingException("Layers do not have same number of features ("
                                         + str(nb_feats_a) + " vs " + str(nb_feats_b) + ")")
        a_fields, b_fields = layer_a.fields().names(), layer_b.fields().names()
        a_field_set, b_field_set = self.getFieldSet(a_fields), self.getFieldSet(b_fields)
        
        include_cbc_a = a_field_set[self.CBC_MESH_SIZE] in a_fields
        include_cbc_b = a_field_set[self.CBC_MESH_SIZE] in b_fields
        # include_cbc = include_cbc_a or include_cbc_b
        if include_cbc_a != include_cbc_b:
            feedback.pushInfo("Comparing CBC with CUT")
            # raise QgsProcessingException("Do not compare results produced in CBC (Cross-Boundary Connection) mode with results produced in CUT mode (default)")
        # Fields
        if self.DIVISOR not in a_fields or self.DIVISOR not in b_fields:
            raise QgsProcessingException("Missing field 'divisor'")
        qgs_fields = self.mkReportFields()
        diff_fields = [self.NB_PATCHES]
        diff_fields_divisor = [self.REPORT_AREA, self.INTERSECTING_AREA, self.NET_PRODUCT]
        # diff_fields_divisor_sq = [self.NET_PRODUCT]
        # same_fields = [self.ID]
        same_fields = []
        if include_cbc_a or include_cbc_b:
            qgs_fields.remove(4)
            qgs_fields.remove(4)
            qgs_fields.remove(4)
            qgs_fields.remove(4)
        else:
            diff_fields += [self.DIVI,self.SPLITTING_INDEX,self.COHERENCE,self.SPLITTING_DENSITY]
        if include_cbc_a:
            a_meff_field, a_np_field = self.CBC_MESH_SIZE, self.CBC_NET_PRODUCT
        else:
            a_meff_field, a_np_field = self.MESH_SIZE, self.NET_PRODUCT
        if include_cbc_b:
            b_meff_field, b_np_field = self.CBC_MESH_SIZE, self.CBC_NET_PRODUCT
        else:
            b_meff_field, b_np_field = self.MESH_SIZE, self.NET_PRODUCT
        # if include_cbc:
            # diff_fields_divisor.append(self.CBC_MESH_SIZE)
            # diff_fields_divisor_sq.append(self.CBC_NET_PRODUCT)
        # fields_names = [f.name() for f in qgs_fields]
        fields_names = []
        # Join layers A and B
        predicates = [2] # 2 <=> join on equal geometries
        joined_path = QgsProcessingUtils.generateTempFilename("joined.gpkg")
        # if qgsUtils.getMaxFid(layer_a) == 1:
            # qgsUtils.normFids(layer_a)
        # if qgsUtils.getMaxFid(layer_b) == 1:
            # qgsUtils.normFids(layer_b)
        joined = qgsTreatments.joinByLoc(layer_a,layer_b,predicates=predicates,
            out_path=joined_path,fields=fields_names,prefix=self.PREFIX,
            context=context,feedback=feedback)
        joined_layer = qgsUtils.loadVectorLayer(joined_path)
        if nb_feats_a != joined_layer.featureCount():
            raise QgsProcessingException("Join by location failed, geometries do not match exactly")
        # Output computation
        qgs_fields.append(QgsField(self.VARIATION,QVariant.Double))
        wkb_type = layer_a.wkbType()
        sink, dest_id = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            qgs_fields,
            wkb_type,
            a_crs)
        for feat in joined_layer.getFeatures():
            new_feat = QgsFeature(qgs_fields)
            new_feat.setGeometry(feat.geometry())
            for fname in same_fields:
                fname_fix = a_field_set[fname]
                new_feat[fname_fix] = feat[fname_fix]
            for fname in diff_fields:
                a_fname = a_field_set[fname]
                b_fname = self.PREFIX + b_field_set[fname]
                a_val, b_val = feat[a_fname], feat[b_fname]
                new_feat[fname] = b_val - a_val
            # Divisor
            if self.DIVISOR in a_fields and self.DIVISOR in b_fields:
                a_divi = feat[self.DIVISOR]
                b_divi = feat[self.PREFIX + self.DIVISOR]
            else:
                a_divi = 1
                b_divi = 1
            # factor = a_divi / b_divi
            # feedback.pushDebugInfo("factor = " + str(factor))
            # factor_sq = factor * factor
            for fname in diff_fields_divisor:
                feedback.pushDebugInfo("fname = " + str(fname))
                a_val = feat[a_field_set[fname]]
                b_val = feat[self.PREFIX + b_field_set[fname]]
                diff_val, sum_vals, divi = self.mkDiff(a_val,b_val,a_divi,b_divi,feedback)
                new_feat[fname] = diff_val
                new_feat[self.DIVISOR] = divi
            # NET PRODUCT
            a_np_val = feat[a_field_set[a_np_field]]
            b_np_val = feat[self.PREFIX + b_field_set[b_np_field]]
            diff_np_val, sum_vals, divi = self.mkDiff(a_np_val,b_np_val,a_divi,b_divi,
                feedback,sq_factor=True)
            new_feat[self.NET_PRODUCT] = diff_np_val
            # MESH SIZE
            a_meff_val = feat[a_field_set[a_meff_field]]
            b_meff_val = feat[self.PREFIX + b_field_set[b_meff_field]]
            diff_meff_val, sum_vals, divi = self.mkDiff(
                a_meff_val,b_meff_val,a_divi,b_divi,feedback)
            new_feat[self.MESH_SIZE] = diff_meff_val
            var_val = diff_meff_val / sum_vals if sum_vals != 0 else 0
            new_feat[self.VARIATION] = round(var_val, self.NB_DIGITS)
            sink.addFeature(new_feat)
        return { self.OUTPUT : dest_id }
        
                
class FragScapeMeffVectorAlgorithm(FragScapeVectorAlgorithm,MeffAlgUtils):
    
    def initAlgorithm(self,config=None,clip_flag=False,out_sink=True):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr("Input layer"),
                [QgsProcessing.TypeVectorPolygon]))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.REPORTING,
                self.tr("Reporting layer"),
                [QgsProcessing.TypeVectorPolygon]))
        self.addParameter(
            QgsProcessingParameterCrs(
                self.CRS,
                description=self.tr("Output CRS"),
                defaultValue=params.defaultCrs))
        self.addParameter(
            QgsProcessingParameterBoolean(
                self.INCLUDE_CBC,
                self.tr("Include Cross-boundary connection metrics")))
        self.addParameter(
            QgsProcessingParameterEnum(
                self.UNIT,
                description=self.tr("Report areas unit"),
                options=self.getUnitOptions()))
        if clip_flag:
            self.addParameter(
                QgsProcessingParameterBoolean(
                    self.CLIP_FLAG,
                    self.tr("Clip input layer at reporting boundaries"),
                    defaultValue=False))
        if out_sink:
            self.addParameter(
                QgsProcessingParameterFeatureSink(
                    self.OUTPUT,
                    self.tr("Output layer")))
        else:
            self.addParameter(
                QgsProcessingParameterVectorDestination(
                    self.OUTPUT,
                    self.tr("Output layer")))
            
                
    def prepareInputs(self,parameters,context,feedback):
        feedbacks.setSubText("Prepare inputs")
        input = self.parameterAsVectorLayer(parameters, self.INPUT, context)
        if input is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.REPORTING))
        reporting = self.parameterAsVectorLayer(parameters,self.REPORTING,context)
        if reporting is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.REPORTING))
        self.crs = self.parameterAsCrs(parameters,self.CRS,context)
        unit = self.parameterAsEnum(parameters,self.UNIT,context)
        self.include_cbc = self.parameterAsBool(parameters,self.INCLUDE_CBC,context)
        self.unit_divisor = self.UNIT_DIVISOR[unit]
        suffix = parameters[self.SUFFIX] if self.SUFFIX in parameters else ""
        # CRS reprojection
        source_crs = input.crs().authid()
        reporting_crs = reporting.crs().authid()
        feedback = feedbacks.ProgressMultiStepFeedback(4,feedback)
        feedback.pushDebugInfo("source_crs = " + str(source_crs))
        feedback.pushDebugInfo("boundary_crs = " + str(reporting_crs))
        feedback.pushDebugInfo("crs = " + str(self.crs.authid()))
        source_name = input.sourceName()
        if source_crs != self.crs.authid():
            source_path = params.mkTmpLayerPath(source_name
                + "_input_reproj" + suffix + ".gpkg")
            qgsTreatments.applyReprojectLayer(input,self.crs,source_path,context,feedback)
            input = qgsUtils.loadVectorLayer(source_path)
            feedback.setCurrentStep(1)
        if reporting_crs != self.crs.authid():
            reporting_path = params.mkTmpLayerPath(source_name
                + "_reporting_reproj" + suffix + ".gpkg")
            qgsTreatments.applyReprojectLayer(reporting,self.crs,reporting_path,context,feedback)
            reporting = qgsUtils.loadVectorLayer(reporting_path)
            feedback.setCurrentStep(2)
        # Dissolve
        if isinstance(self,MeffVectorGlobal) and reporting.featureCount() > 1:
            reporting_path = params.mkTmpLayerPath(source_name + "_reporting_dissolved.gpkg")
            qgsTreatments.dissolveLayer(reporting,reporting_path,context,feedback)
            reporting = qgsUtils.loadVectorLayer(reporting_path)
            feedback.setCurrentStep(3)
        self.report_layer = reporting
        # Clip by boundary
        out_path = params.mkTmpLayerPath(source_name
            + "_prepared" + suffix + ".gpkg")
        if self.include_cbc:
            qgsTreatments.selectIntersection(input,reporting,context,feedback)
            qgsTreatments.saveSelectedFeatures(input,out_path,context,feedback)
        elif self.CLIP_FLAG in parameters and not parameters[self.CLIP_FLAG]:
            clipped_path = params.mkTmpLayerPath(source_name
                + "_intersected" + suffix + ".gpkg")
            qgsTreatments.applyIntersection(input,reporting,clipped_path,context,feedback)
            qgsTreatments.multiToSingleGeom(clipped_path,out_path,context,feedback)
        else:
            clipped_path = params.mkTmpLayerPath(source_name
                + "_clipped" + suffix + ".gpkg")
            qgsTreatments.applyVectorClip(input,reporting,clipped_path,context,feedback)
            qgsTreatments.multiToSingleGeom(clipped_path,out_path,context,feedback)
        feedback.setCurrentStep(4)
        input = qgsUtils.loadVectorLayer(out_path)
        return (input, reporting)


class MeffVectorGlobal(FragScapeMeffVectorAlgorithm):

    ALG_NAME = "meffVectorGlobal"
    
    # OUTPUT_GLOBAL_MEFF = "GLOBAL_MEFF"
        
    def displayName(self):
        return self.tr("Vector Effective Mesh Size (Global)")
        
    def shortHelpString(self):
        return self.tr("Computes effective mesh size from patch layer and boundary of reporting layer (features are dissolved if needed)")
                
    def processAlgorithm(self,parameters,context,feedback):
        feedback.pushDebugInfo("Start " + str(self.name()))
        # Parameters
        source, boundary = self.prepareInputs(parameters, context, feedback)
        output = parameters[self.OUTPUT]
        # CRS reprojection
        # source_crs = source.crs().authid()
        # boundary_crs = boundary.crs().authid()
        # feedback.pushDebugInfo("source_crs = " + str(source_crs))
        # feedback.pushDebugInfo("boundary_crs = " + str(boundary_crs))
        # feedback.pushDebugInfo("crs = " + str(self.crs.authid()))
        # source_name = boundary.sourceName()
        # if source_crs != self.crs.authid():
            # source_path = params.mkTmpLayerPath(source_name + "_source_reproject.gpkg")
            # qgsTreatments.applyReprojectLayer(source,self.crs,source_path,context,feedback)
            # source = qgsUtils.loadVectorLayer(source_path)
        # if boundary_crs != self.crs.authid():
            # boundary_path = params.mkTmpLayerPath(source_name + "_boundary_reproject.gpkg")
            # qgsTreatments.applyReprojectLayer(boundary,self.crs,boundary_path,context,feedback)
            # boundary = qgsUtils.loadVectorLayer(boundary_path)
        # Clip by boundary
        # intersected_path = params.mkTmpLayerPath(source_name + "_source_intersected.gpkg")
        # if self.include_cbc:
            # qgsTreatments.selectIntersection(source,boundary,context,feedback)
            # qgsTreatments.saveSelectedFeatures(source,intersected_path,context,feedback)
        # selected_path = intersected_path
        # source = qgsUtils.loadVectorLayer(selected_path)
        # Dissolve
        # if boundary.featureCount() > 1:
            # dissolved_path = params.mkTmpLayerPath(source_name + "_boundary_dissolved.gpkg")
            # qgsTreatments.dissolveLayer(boundary,dissolved_path,context,feedback)
            # boundary = qgsUtils.loadVectorLayer(dissolved_path)
            # self.report_layer = boundary
        # Algorithm
        # progress step
        nb_feats = source.featureCount()
        feedbacks.setSubText("Iterating on features")
        feedback.pushDebugInfo("nb_feats = " + str(nb_feats))
        if nb_feats == 0:
            utils.warn("Empty input layer : " + qgsUtils.pathOfLayer(source))
        else:
            feedback = feedbacks.ProgressMultiStepFeedback(nb_feats,feedback)
        # Reporting area
        for report_feat in boundary.getFeatures():
            report_geom = report_feat.geometry()
        report_area = report_geom.area()
        sum_ai = 0
        feedback.pushDebugInfo("report_area = " + str(report_area))
        if report_area == 0:
            raise QgsProcessingException("Empty reporting area")
        else:
            feedback.pushDebugInfo("ok")
        net_product = 0
        cbc_net_product = 0
        intersecting_area = 0
        feats_cpt = 0
        for cpt, f in enumerate(source.getFeatures()):
            f_geom = f.geometry()
            if f_geom.intersects(report_geom):
                intersection = f_geom.intersection(report_geom)
                f_area = f_geom.area()
                # sum_ai += f_area
                for part_geom in intersection.parts():
                    feats_cpt += 1
                    part_area = part_geom.area()
                    intersecting_area += part_area
                    net_product += part_area * part_area
                    cbc_net_product += f_area * part_area
            feedback.setProgress(cpt)
        report_area_sq = report_area * report_area
        # Outputs
        res_dict = { self.REPORT_AREA : report_area,
            self.SUM_AI : intersecting_area,
            self.SUM_AI_SQ : net_product,
            self.NB_PATCHES : feats_cpt,
            self.DIVISOR : self.unit_divisor,
        }
        if self.include_cbc:
            res_dict[self.SUM_AI_SQ_CBC] = cbc_net_product
        res_layer, res_val = self.mkOutputs(parameters,res_dict,context)
        return {self.OUTPUT: output, self.OUTPUT_VAL : res_val}

   

class MeffVectorReport(FragScapeMeffVectorAlgorithm):

    ALG_NAME = "meffVectorReport"
        
    def displayName(self):
        return self.tr("Vector Effective Mesh Size per feature")
        
    def shortHelpString(self):
        return self.tr("Computes effective mesh size from patch layer for each feature of reporting layer.")
                
    
    def initAlgorithm(self,config=None):
        super().initAlgorithm(clip_flag=True,out_sink=False)
                
    def processAlgorithm(self,parameters,context,feedback):
        source, reporting = self.prepareInputs(parameters, context, feedback)
        output = parameters[self.OUTPUT]
        out = self.parameterAsOutputLayer(parameters,self.OUTPUT,context)
        # report_fields = self.mkReportFields(self.include_cbc)
        # wkb_type = reporting.wkbType()
        # report_crs = reporting.sourceCrs()
        # sink, dest_id = self.parameterAsSink(
            # parameters,
            # self.OUTPUT,
            # context,
            # report_fields,
            # wkb_type,
            # report_crs)
        # Algorithm
        # progress step
        nb_feats = reporting.featureCount()
        feedback.pushDebugInfo("nb_feats = " + str(nb_feats))
        if nb_feats == 0:
            raise QgsProcessingException("Empty layer")
        multi_feedback = feedbacks.ProgressMultiStepFeedback(nb_feats, feedback)
        # Global
        # dissolved_path = params.mkTmpLayerPath('reportingDissolved.gpkg')
        # qgsTreatments.dissolveLayer(reporting,dissolved_path,context,feedback)
        # if nb_feats == 1:
            # out_path = output
        # elif self.include_cbc:
            # out_path = params.mkTmpLayerPath("reportingComputedGlobalCBC"".gpkg")
        # else:
            # out_path = params.mkTmpLayerPath("reportingComputedGlobalCUT"".gpkg")
        # parameters = { MeffVectorGlobal.INPUT : source,
                   # MeffVectorGlobal.REPORTING : dissolved_path,
                   # MeffVectorGlobal.CRS : self.crs,
                   # MeffVectorGlobal.INCLUDE_CBC : self.include_cbc,
                   # MeffVectorGlobal.UNIT : parameters[self.UNIT],
                   # MeffVectorGlobal.OUTPUT : out_path,
                   # self.SUFFIX : "Global" }
        # global_res = qgsTreatments.applyProcessingAlg('FragScape', 
            # MeffVectorGlobal.ALG_NAME,parameters, context, multi_feedback,
            # onlyOutput=False)
        # multi_feedback.setCurrentStep(1)
        # Iterate on features
        report_layers = []
        parameters[MeffVectorGlobal.INPUT] = source
        source_nb_feats = source.featureCount()
        feedback.pushDebugInfo("source_nb_feats = " + str(source_nb_feats))
        if source_nb_feats >= 1:
            all_ids = set()
            for count, report_feat in enumerate(reporting.getFeatures(),start=2):
                report_id = report_feat.id()
                if report_id in all_ids:
                    raise QgsProcessingException("Identifier of reporting layer not unique : "
                        + str(report_id))
                reporting.selectByIds([report_id])
                select_path = params.mkTmpLayerPath("reportingSelection"
                    + str(report_feat.id()) + ".gpkg")
                qgsTreatments.saveSelectedFeatures(reporting,select_path,context,multi_feedback)
                report_computed_path = params.mkTmpLayerPath("reportingComputed"
                    + str(report_feat.id()) + ".gpkg")
                parameters[MeffVectorGlobal.REPORTING] = select_path
                parameters[MeffVectorGlobal.OUTPUT] = report_computed_path
                parameters[self.SUFFIX] = str(report_id)
                qgsTreatments.applyProcessingAlg('FragScape',
                                                 MeffVectorGlobal.ALG_NAME,
                                                 parameters,context,multi_feedback)
                report_layers.append(report_computed_path)
                multi_feedback.setCurrentStep(count)
            # merged_path = params.mkTmpLayerPath("reportingMerged.gpkg")
            feedback.pushDebugInfo("output = " + str(output))
            feedback.pushDebugInfo("output class = " + str(output.__class__))
            qgsTreatments.mergeVectorLayers(report_layers,self.crs,out)
            # if isinstance(output,str):
                # out_layer = qgsUtils.loadVectorLayer(output)
                # out_layer.dataProvider().deleteAttributes([0])
            # else:
            feedback.pushDebugInfo("out_layer = " + str(out))
            out_layer = qgsUtils.loadVectorLayer(out)
            out_layer.dataProvider().deleteAttributes([0])
                
        return { self.OUTPUT: output }#, self.OUTPUT_VAL : global_res[self.OUTPUT_VAL] }

