# -*- coding: utf-8 -*-

"""
/***************************************************************************
 ForestRoads
                                 A QGIS plugin
 Create a network of forest roads based on zones to access, roads to connect
 them to, and a cost matrix.
 The code of the plugin is based on the "LeastCostPath" plugin available on
 https://github.com/Gooong/LeastCostPath. We thank their team for the template.
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 10-07-2019
        copyright            : (C) 2019 by Clement Hardy
        email                : clem.hardy@outlook.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.                                   *
 *                                                                         *
 ***************************************************************************/
 This script describes the algorithm used to determine the wood flux
 in the network.
"""

__author__ = 'clem.hardy@outlook.fr'
__date__ = 'Currently in work'
__copyright__ = '(C) 2019 by Clement Hardy'

# We load every function necessary from the QIS packages.
import random
import queue
import math
import processing
from .kdtree import KDTree
import numpy as np
from PyQt5.QtCore import QCoreApplication, QVariant
from PyQt5.QtGui import QIcon
from qgis.core import (
    QgsFeature,
    QgsGeometry,
    QgsPoint,
    QgsPointXY,
    QgsField,
    QgsFields,
    QgsWkbTypes,
    QgsProcessing,
    QgsFeatureSink,
    QgsProcessingException,
    QgsProcessingAlgorithm,
    QgsProcessingParameterFeatureSource,
    QgsProcessingParameterFeatureSink,
    QgsProcessingParameterRasterLayer,
    QgsProcessingParameterField,
    QgsProcessingParameterBand,
    QgsProcessingParameterBoolean,
    QgsProcessingParameterNumber,
    QgsProcessingParameterEnum
)

# The algorithm class heritates from the algorithm class of QGIS.
# There, it can register different parameter during initialization
# that can be put into variables using "
class woodFluxAlgorithm(QgsProcessingAlgorithm):
    """
    Class that described the algorithm, that will be launched
    via the provider, itself launched via initialization of
    the plugin.
    """

    # Constants used to refer to parameters and outputs. They will be
    # used when calling the algorithm from another algorithm, or when
    # calling from the QGIS console.

    INPUT_ROAD_NETWORK = 'INPUT_ROAD_NETWORK'

    INPUT_POLYGONS_CUTTED = 'INPUT_POLYGONS_CUTTED'

    POINTS_RESOLUTION = 'POINTS_RESOLUTION'

    WOOD_ATTRIBUTE = 'WOOD_ATTRIBUTE'

    DENSITY_OR_VOLUME = 'DENSITY_OR_VOLUME'

    WOOD_DENSITY = 'WOOD_DENSITY'

    INPUT_ENDING_POINTS = 'INPUT_ENDING_POINTS'

    RETURN_EMPTY_ROADS = 'RETURN_EMPTY_ROADS'

    OUTPUT = 'OUTPUT'

    def initAlgorithm(self, config):
        """
        Here we define the inputs and output of the algorithm, along
        with some other properties. Theses will be asked to the user.
        """

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT_ROAD_NETWORK,
                self.tr('Forest Road network whose road flux must be determined'),
                [QgsProcessing.TypeVectorLine]
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT_POLYGONS_CUTTED,
                self.tr('Polygons where the wood was cut'),
                [QgsProcessing.TypeVectorPolygon]
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.POINTS_RESOLUTION,
                self.tr('The density of points (for 1 crs units) that will be treated as a source of wood.'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=1,
                optional=False,
                minValue=0.0000000001
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.WOOD_ATTRIBUTE,
                self.tr('Attribute field that indicates quantities of wood in each polygon (density or volume)'),
                parentLayerParameterName=self.INPUT_POLYGONS_CUTTED,
                type=QgsProcessingParameterField.Numeric,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.DENSITY_OR_VOLUME,
                self.tr('Is the chosen attribute field a measure of density or volume of wood ?'),
                ['Density', 'Volume'],
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.WOOD_DENSITY,
                self.tr('If no attribute field is selected, indicate an arbitrary density of wood for all polygons'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=100,
                optional=False,
                minValue=0.0000000001
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT_ENDING_POINTS,
                self.tr('The end points of the network (to where the wood must go)'),
                [QgsProcessing.TypeVectorPoint]
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.RETURN_EMPTY_ROADS,
                self.tr('Should we save a road if its flux is 0 ?'),
                ['Yes', 'No'],
                defaultValue = 1
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr('Output of the algorithm')
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        """
        Here is where the processing itself takes place.
        """

        # First, we register all of the parameters entered by the user into variables.
        road_network = self.parameterAsVectorLayer(
            parameters,
            self.INPUT_ROAD_NETWORK,
            context
        )

        polygons_cutted = self.parameterAsVectorLayer(
            parameters,
            self.INPUT_POLYGONS_CUTTED,
            context
        )

        points_resolution = self.parameterAsDouble(
            parameters,
            self.POINTS_RESOLUTION,
            context
        )

        # If the user entered an attribute field for information on the wood present in the polygons, we read it
        # and we indicate that we will not use the arbitrary density; if not, we do the opposite.
        if self.parameterAsString(
                    parameters,
                    self.WOOD_ATTRIBUTE,
                    context
                ) is not None:
            wood_attribute_index = polygons_cutted.fields().lookupField(
                self.parameterAsString(
                    parameters,
                    self.WOOD_ATTRIBUTE,
                    context
                )
            )
            density_or_volume = self.parameterAsString(
                parameters,
                self.DENSITY_OR_VOLUME,
                context
            )
            # If the user indicated an attribute field, but did not indicate if it corresponds to a volume or a density
            # of wood, we raise an exception.
            if density_or_volume is None or density_or_volume == '':
                raise QgsProcessingException(self.tr("Please, describe if the attribute field you have chosen for the"
                                                     " polygons corresponds to density or volume of wood."))
            wood_density = None

        else:
            wood_density = self.parameterAsDouble(
                parameters,
                self.WOOD_DENSITY,
                context
            )
            wood_attribute_index = None
            density_or_volume = None

        ending_points = self.parameterAsVectorLayer(
            parameters,
            self.INPUT_ENDING_POINTS,
            context
        )

        if self.parameterAsString(
            parameters,
            self.RETURN_EMPTY_ROADS,
            context
        ) == '0':
            return_empty_roads = True
        else:
            return_empty_roads =  False

        # If source was not found, throw an exception to indicate that the algorithm
        # encountered a fatal error.
        if road_network is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_ROAD_NETWORK))
        if polygons_cutted is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_POLYGONS_CUTTED))
        if ending_points is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_ENDING_POINTS))

        # We try to see if there are divergence between the CRSs of the inputs
        if road_network.crs() != ending_points.sourceCrs() or road_network.crs() != polygons_cutted.sourceCrs()\
                or ending_points.crs() != polygons_cutted.sourceCrs():
            raise QgsProcessingException(self.tr("ERROR: The input layers have different CRSs."))

        # We initialize the "sink", an object that will make use able to create an output.
        # First, we create the fields for the attributes of our lines as outputs.
        # They will only have one field :
        sink_fields = WoodFluxHelper.create_fields()
        # We indicate that our output will be a line, stored in WBK format.
        output_geometry_type = QgsWkbTypes.LineString
        # Finally, we create the field object and register the destination ID of it.
        # This sink will be based on the OUTPUT parameter we registered during initialization,
        # will have the fields and the geometry type we just created, and the same CRS as the cost raster.
        (sink, dest_id) = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            fields=sink_fields,
            geometryType=output_geometry_type,
            crs=road_network.crs(),
        )

        # If sink was not created, throw an exception to indicate that the algorithm
        # encountered a fatal error. The exception text can be any string, but in this
        # case we use the pre-built invalidSinkError method to return a standard
        # helper text for when a sink cannot be evaluated
        if sink is None:
            raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

        # Before all, thanks to the polygon extent and the wood density parameter, we will define the "points" inside the
        # polygons from which the wood flux will come from. We do it first so that if it fails due to the input point
        # resolution, the user can try another one quickly.
        feedback.pushInfo(self.tr("Analyzing harvested areas..."))
        feedback.pushInfo(
            self.tr("If this step takes too long, it might be that the point resolutions you gave as input"
                    " is too high, resulting in a very large number of points as source of wood generated."))

        if density_or_volume == '1':
            isItVolume = True
            # feedback.pushInfo(self.tr("Volume detected."))
        else:
            isItVolume = False
            # feedback.pushInfo(self.tr("Density detected."))

        pointsOfWoodGeneration, woodVolumePerPointDictionary = WoodFluxHelper.polygonsToPoints(polygons_cutted,
                                                                                               points_resolution,
                                                                                               wood_density,
                                                                                               isItVolume,
                                                                                               wood_attribute_index,
                                                                                               feedback)

        # Then, because this algorithm needs the line network to be cut at each intersection., we will cut it.
        # The user might have done it already; but just in case, we'll do it again thanks to one a
        # QGIS processing algorithm.
        feedback.pushInfo(self.tr("Splitting lines..."))
        splittedLinesResult = processing.run("native:splitwithlines",
                                      {'INPUT': road_network,
                                       'LINES': road_network,
                                       'OUTPUT': 'memory:'})
        splittedLines = splittedLinesResult['OUTPUT']

        feedback.pushInfo(self.tr("Scanning lines..."))
        # For each line feature in the road network, we initialize an object of class "line"
        multi_line = set(splittedLines.getFeatures())
        setOfRoadsAsLines = set()
        ID = 0
        # We are careful to treat each individual line of the given layer;
        # for that, we have to treat the case of a multiline geometry, or a simple line.
        for featureLine in multi_line:
            featureLineGeo = featureLine.geometry()

            # Case of multiline
            if featureLineGeo.wkbType() == QgsWkbTypes.MultiLineString:
                multi_line = featureLineGeo.asMultiPolyline()
                for line in multi_line:
                    setOfRoadsAsLines.add(LineForAlgorithm(line, ID))
                    ID += 1

            # Case of lines
            elif featureLineGeo.wkbType() == QgsWkbTypes.LineString:
                line = featureLineGeo.asPolyline()
                setOfRoadsAsLines.add(LineForAlgorithm(line, ID))
                ID += 1

            if feedback.isCanceled():
                raise QgsProcessingException(self.tr("ERROR: Operation was cancelled."))

        # Then, we transform the ending points of the network into a list of QgsPointXY we can compare to our lines.
        # If the user have input a layer that is not a multipoint layer, we will first promote it to multipoint.
        endingPointsAsMultipointsResult = processing.run("native:promotetomulti",
                                             {'INPUT': ending_points,
                                              'OUTPUT': 'memory:'})
        endingPointsAsMultipoints = endingPointsAsMultipointsResult['OUTPUT']
        tempSetOfEndingPoints = endingPointsAsMultipoints.getFeatures()
        setOfMultiEndingPoints = [endingPoint.geometry().asMultiPoint() for endingPoint in list(tempSetOfEndingPoints)]
        setOfEndingPoints = set()
        for multiEndingPoints in setOfMultiEndingPoints:
            for point in multiEndingPoints:
                setOfEndingPoints.add(point)

        # Then, we update our line objects to reflect network structure. See function in Line class.
        feedback.pushInfo(self.tr("Scanning network structure..."))

        # First, we got to initialize some k-d tree for some quick nearest-point searches, and the dictionary
        # relating extremities of lines to lines.
        endingPointsOfLines = set()
        linesEndingDictionary = dict()
        for line in setOfRoadsAsLines:
            endingPointsOfLines.add(line.ending1)
            if line.ending1 not in linesEndingDictionary:
                linesEndingDictionary[line.ending1] = list()
                linesEndingDictionary[line.ending1].append(line)
            else:
                linesEndingDictionary[line.ending1].append(line)
            endingPointsOfLines.add(line.ending2)
            if line.ending2 not in linesEndingDictionary:
                linesEndingDictionary[line.ending2] = list()
                linesEndingDictionary[line.ending2].append(line)
            else:
                linesEndingDictionary[line.ending2].append(line)

        endingPointsOfLines = list(endingPointsOfLines)
        numpyArrayOfLinesEndingPoints = np.array(endingPointsOfLines)
        kdTreeForLinesEnding = KDTree(numpyArrayOfLinesEndingPoints, leafsize=20)
        numpyArrayOfEndingPoints = np.array(list(setOfEndingPoints))
        kdTreeForEndingPoints = KDTree(numpyArrayOfEndingPoints, leafsize=20)

        feedbackProgress = 0
        for line in setOfRoadsAsLines:
            line.initializeRelationToNetwork(endingPointsOfLines, kdTreeForLinesEnding,
                                             linesEndingDictionary, kdTreeForEndingPoints)
            if feedback.isCanceled():
                raise QgsProcessingException(self.tr("ERROR: Operation was cancelled."))
            feedbackProgress += 1
            feedback.setProgress(100 * (feedbackProgress / len(setOfRoadsAsLines)))

        # Once it's done, we have to give a direction to each line object. Again, see function in Line class.
        feedback.pushInfo(self.tr("Scanning network direction..."))
        feedbackProgress = 0
        for line in setOfRoadsAsLines:
            line.initializeDownstreamDirection(setOfEndingPoints, feedback)
            if feedback.isCanceled():
                raise QgsProcessingException(self.tr("ERROR: Operation was cancelled."))
            feedbackProgress += 1
            feedback.setProgress(100 * (feedbackProgress / len(setOfRoadsAsLines)))

        # Then, for each of the wood source points, we find the closest line; we add units of timber flux to this line.
        feedback.pushInfo(self.tr("Generating wood flux from harvested areas..."))
        WoodFluxHelper.generateWoodFlux(setOfRoadsAsLines, pointsOfWoodGeneration, woodVolumePerPointDictionary, feedback)

        # Now, we can compute the fluxes that go through the network.
        # We initialize the objects for the loop.
        frontier = queue.Queue()
        setOfLinesToReturn = set()
        for line in setOfRoadsAsLines:
            frontier.put(line)

        feedback.pushInfo(self.tr("Commencing flux algorithm..."))
        # We loop : The loop will end when no more line is in an "open" status.
        # WARNING : if the directions of the lines are not correct, this will become an infinite loop !
        linesAlreadyConsidered = set()
        while not frontier.empty():

            if feedback.isCanceled():
                raise QgsProcessingException(self.tr("ERROR: Operation was cancelled."))

            line = frontier.get()
            # If the line have be fluxed by all of its upstream neighbors,
            # we make it flux to its downstream neighbors, and we put the neighbors in the
            # frontier
            if line.checkIfAllUpstreamHaveFluxed():
                line.fluxToDownstreamNeighbors(feedback)
                setOfLinesToReturn.add(line)
                linesAlreadyConsidered.add(line)
                neighbors = line.getNeighborsDownstream()
                if len(neighbors) != 0:
                    for neighbour in [neighbour for neighbour in neighbors if neighbour not in linesAlreadyConsidered]:
                        # We open the neighbour downstream if they are not already considered/fluxed.
                        frontier.put(neighbour)
            # If not, we re-add the line in the frontier. It will have to be considered later, once
            # that all of its neighbors are fluxed.
            # WARNING : This is where the infinite loop can generate if a culvert is created.
            else:
                frontier.put(line)

            feedback.setProgress(100 * (len(setOfLinesToReturn) / len(setOfRoadsAsLines)))

        feedback.pushInfo(self.tr("Preparing outputs..."))
        # Once that all of this is done, we prepare the output.
        # For every path we create, we save it as a line and put it into the sink !
        for line in setOfLinesToReturn:
            # We only return the line if it's flux is superior to 0, depending on the
            # instruction given by the user
            if return_empty_roads or (not return_empty_roads and line.flux > 0):
                # With the total cost which is the last item in our accumulated cost list,
                # we create the PolyLine that will be returned as a vector.
                path_feature = WoodFluxHelper.create_path_feature_from_points(line.lineFeature,
                                                                              line.uniqueID,
                                                                              line.flux,
                                                                              line.getNeighborsUpstreamID(),
                                                                              line.getNeighborsDownstreamID(),
                                                                              sink_fields, )
                # Into the sink that serves as our output, we put the PolyLines from the list of lines we created
                # one by one
                sink.addFeature(path_feature, QgsFeatureSink.FastInsert)

        # When all is done, we return our output that is linked to the sink.
        feedback.pushInfo("Lines detected in input = " + str(len(setOfRoadsAsLines))
                          + "; Lines in output  " + str(len(setOfLinesToReturn)))
        return {self.OUTPUT: dest_id}

    # Here are different functions used by QGIS to name and define the algorithm
    # to the user.
    def name(self):
        """
        Returns the algorithm name, used for identifying the algorithm. This
        string should be fixed for the algorithm, and must not be localised.
        The name should be unique within each provider. Names should contain
        lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        """
        return 'Wood Flux Determination'

    def displayName(self):
        """
        Returns the translated algorithm name, which should be used for any
        user-visible display of the algorithm name.
        """
        return self.tr(self.name())

    def group(self):
        """
        Returns the name of the group this algorithm belongs to. This string
        should be localised.
        """
        return self.tr(self.groupId())

    def groupId(self):
        """
        Returns the unique ID of the group this algorithm belongs to. This
        string should be fixed for the algorithm, and must not be localised.
        The group id should be unique within each provider. Group id should
        contain lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        """
        # Not enough algorithms in this plugin for the need to make groups of algorithms
        return ''

    # Function used for translation. Called every time something needs to be
    # Translated
    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return woodFluxAlgorithm()

    def helpUrl(self):
        return 'https://github.com/Klemet/ForestRoadNetworkPluginForQGIS'

    def shortHelpString(self):
        return self.tr("""
        This algorithm determine how the wood flux will go throughout a forest road network. To that end, it generates "wood source points" inside polygons where wood have been cut, and "flows" the wood toward the nearest exit point of the network.
        
        **Parameters:**
          
          Please ensure all the input layers have the same CRS.
        
          - Forest road network whose wood flux must be determined : straight-forward.
          
          - Polygons where the wood was cut : Polygons that correspond to the zones where timber will or have been harvested.
          
          - Density of points : The number of "wood-source-points" that will be generated inside the polygons for 1 unit of CRS. If the algorithm doesn't succeed in generating at least one point in a given polygon, it will trying again with a higher resolution until it succeeds. If the density is too low, some smaller roads in the network might have a flux of 0.
          
          - Attribute field for volume or density : An attribute field in the polygon layer where the wood was cut that contains information on the quantity of wood that is inside each polygon.
          
          - Density or volume of wood : Indicate if the measure in the previously given attribute field is a measure of volume of wood, or density of wood. If it is density, it will be considered that it is expressed in CRS units.
          
          - Arbitrary density : If no attribute field with information about the quantity of wood in the polygons has been previously given, this value will be used as a default density of wood for each of the polygons where the wood was cut. It is expressed in CRS units. 
         
          - Ending points : points that correspond EXACTLY to the ending of the network, meaning its connection to the main road network. WARNING : If those points do not correspond exactly with the end of a line or the end of your network, or if there is not at least one such point per "branch" of your network, the algorithm will have problems to complete. To generate them, you can extract the nodes of your lines with the "Extract vertices" tool in QGIS; then, you can create a small buffer around the areas where the timber must go to, and use this buffer to select the nodes inside of it. These nodes will be your ending points.
            
          - Return empty roads : If a road is calculated to have a wood flux of 0, should it be saved in the output or be put aside ?        
        """)

    def shortDescription(self):
        return self.tr('Determine the way the timber will "flow" through a forest road network.')

    # Path to the icon of the algorithm
    def svgIconPath(self):
        return '.icon.png'

    def tags(self):
        return ['wood', 'flux', 'flow', 'harvest', 'roads', 'analysis', 'road', 'network', 'forest']


# The "Line" object used in the algorithm
class LineForAlgorithm:
    """Class used for the algorithm. To initialize, a line feature and a unique ID must be provided"""

    def __init__(self, lineFeature, ID):
        """Constructor for the line class"""
        self.uniqueID = ID  # Used to identify the line.
        self.lineFeature = lineFeature  # As Polyline gives a list of points
        # Length of the line, which is the sum of the distance between the ordered points
        self.lineLength = 0
        for i in range(0, len(lineFeature) - 1):
            self.lineLength = self.lineLength + lineFeature[i].distance(lineFeature[i+1])
        # Ending 1 and 2 are QgsPointXY.
        self.ending1 = self.lineFeature[0]
        self.ending2 = self.lineFeature[-1]
        # These sets are temporary, and replaced by Downstream/upstream once they have been identified.
        self.linesConnectedToEnding1 = set()
        self.linesConnectedToEnding2 = set()
        # Just boxes to put the above sets in, but properly identified.
        self.linesConnectedDownstream = set()
        self.linesConnectedUpstream = set()
        # The flux of wood going through our road
        self.flux = 0
        # Indicators to know if this line is a leaf or a root of the network
        self.isALeafOfTheNetwork = False
        self.isARootOfTheNetwork = False
        # Indicator to know if the line have already given its flux to its downstream neighbors, so that it
        # doesn't do it again by mistake.
        self.haveFluxed = False

    def __lt__(self, other):
        """"Function needed to compare two lines, if there priority in a priority queue are equal. Usefull for the
        dijkstra search below."""
        return self.uniqueID < other.uniqueID

    def initializeRelationToNetwork(self, listOfLineEndings, kdTreeForLinesEnding, linesEndingDictionary, kdTreeForEndingPoints):
        """Initialize the lines that are connected to this line, and we also check if the line is a "leaf" or not.
        Carefull : the linesEndingDictionary must contain a list for each entry. This list can contain several
        lines, as two lines can share the same exact ending point."""

        # We get all of the points at a distance of 1 CRS units around the ending 1 of the line.
        indexOfEndingsOfLinesNearby = kdTreeForLinesEnding.query_ball_point(self.ending1, 1)
        if len(indexOfEndingsOfLinesNearby) > 0:
            for indexOfEndingPointOfLine in indexOfEndingsOfLinesNearby:
                otherLines = linesEndingDictionary[listOfLineEndings[indexOfEndingPointOfLine]]
                for otherLine in otherLines:
                    if otherLine.uniqueID != self.uniqueID:
                        self.linesConnectedToEnding1.add(otherLine)

        # We do the same for the second ending.
        indexOfEndingsOfLinesNearby = kdTreeForLinesEnding.query_ball_point(self.ending2, 1)
        if len(indexOfEndingsOfLinesNearby) > 0:
            for indexOfEndingPointOfLine in indexOfEndingsOfLinesNearby:
                otherLines = linesEndingDictionary[listOfLineEndings[indexOfEndingPointOfLine]]
                for otherLine in otherLines:
                    if otherLine.uniqueID != self.uniqueID:
                        self.linesConnectedToEnding2.add(otherLine)

        # To check if the line is a leaf (an outward ending of the network) or a root, we check if it has an ending
        # without neighbors, and if this ending coincide with an end point. If it does not for any ending point,
        # then the line is a leaf of the network.

        # For that, we use another k-d tree to search for the ending points that could be nearby.
        if len(self.linesConnectedToEnding1) == 0:
            endingPointsNearby = kdTreeForEndingPoints.query_ball_point(self.ending1, 1)
            if len(endingPointsNearby) > 0:
                self.isARootOfTheNetwork = True
                return
            else:
                self.isALeafOfTheNetwork = True

        elif len(self.linesConnectedToEnding2) == 0:
            endingPointsNearby = kdTreeForEndingPoints.query_ball_point(self.ending2, 1)
            if len(endingPointsNearby) > 0:
                self.isARootOfTheNetwork = True
                return
            else:
                self.isALeafOfTheNetwork = True

    def initializeDownstreamDirection(self, endingPoints, feedback):
        """This function will find and register which ending is downstream, and which is upstream.
           We use a kind of dijkstra search for that."""
        # First, we make the case of our line being a leaf or a root. If that's the case, it's relatively
        # easy to deduce where is the downstream or the upstream.
        if self.isALeafOfTheNetwork:
            if len(self.linesConnectedToEnding1) == 0:
                self.linesConnectedDownstream = self.linesConnectedToEnding2
                self.linesConnectedUpstream = self.linesConnectedToEnding1
            else:
                self.linesConnectedDownstream = self.linesConnectedToEnding1
                self.linesConnectedUpstream = self.linesConnectedToEnding2
        elif self.isARootOfTheNetwork:
            if len(self.linesConnectedToEnding1) == 0:
                self.linesConnectedDownstream = self.linesConnectedToEnding1
                self.linesConnectedUpstream = self.linesConnectedToEnding2
            else:
                self.linesConnectedDownstream = self.linesConnectedToEnding2
                self.linesConnectedUpstream = self.linesConnectedToEnding1
        else:
            # Now for the case of a line which is not a root, or a leaf.
            # To find upstream and downstream, we randomly take one of the endings of our line (the ending1).
            # Then, we start a kind of dijkstra search. If we find the ending point at only one ending of the line,
            # Then this ending is downward. If not, the downward ending will be the one for which the least-cost path
            # will be the shortest.

            # We first search the least-cost path from the ending1
            resultsForEnding1 = self._FindLeastCostPathFromEnding(self.linesConnectedToEnding1, feedback)

            # We then do the same with ending2
            resultsForEnding2 = self._FindLeastCostPathFromEnding(self.linesConnectedToEnding2, feedback)

            # We check the different cases if the ending points have been found
            if(resultsForEnding1["found"] == True and resultsForEnding2["found"] == True):
                if(resultsForEnding1["distance"] < resultsForEnding2["distance"]):
                    self.linesConnectedDownstream = self.linesConnectedToEnding1
                    self.linesConnectedUpstream = self.linesConnectedToEnding2
                elif(resultsForEnding1["distance"] > resultsForEnding2["distance"]):
                    self.linesConnectedDownstream = self.linesConnectedToEnding2
                    self.linesConnectedUpstream = self.linesConnectedToEnding1
                # Rare case of if the distances where the same. In that case, ending1 will become the downstream one by default.
                else:
                    self.linesConnectedDownstream = self.linesConnectedToEnding1
                    self.linesConnectedUpstream = self.linesConnectedToEnding2
            elif((resultsForEnding1["found"] == True and resultsForEnding2["found"] == False)):
                self.linesConnectedDownstream = self.linesConnectedToEnding1
                self.linesConnectedUpstream = self.linesConnectedToEnding2
            elif((resultsForEnding1["found"] == False and resultsForEnding2["found"] == True)):
                self.linesConnectedDownstream = self.linesConnectedToEnding2
                self.linesConnectedUpstream = self.linesConnectedToEnding1
            else:
                raise QgsProcessingException(
                    "ERROR: For a certain line, Upstream and Downstream could not be determined." +
                    "No root of the network have been found." +
                    "Problematic line is between those two points : " + str(self.ending1)
                    + " and " + str(self.ending2) + ". Please check if an ending point exist for every branch"
                                                    "of your network..")

            # If the function ended without use being able to determine at least a upstream or a downstream,
            # this could be a problem. We raise an exception with details about the line.
            if len(self.linesConnectedDownstream) == 0 and len(self.linesConnectedUpstream) == 0:
                raise QgsProcessingException(
                    "ERROR: For a certain line, Upstream and Downstream could not be determined."
                    "Problematic line is between those two points : " + str(self.ending1)
                    + " and " + str(self.ending2) + ". Please check if an ending point exist for every branch"
                                                    "of your network..")
            if len(self.linesConnectedDownstream) == 0 and not self.isARootOfTheNetwork:
                raise QgsProcessingException(
                    "ERROR: For a certain line that is not a root, Downstream could not be determined."
                    "Problematic line is between those two points : " + str(self.ending1)
                    + " and " + str(self.ending2) + ". Please check if an ending point exist for every branch"
                                                    "of your network..")


    def _FindLeastCostPathFromEnding(self, linesConnectedToEnding, feedback):
        """This function is here to help with the previous function to do the pathfinding to a ending point of
        the network."""

        # We make a dictionary that will contain our results; and one that will contain the predecessors of the lines
        # for the dijkstra search
        results = dict()
        predecessors = dict()
        costSoFar = dict()
        costSoFar[self] = self.lineLength
        # We initialize a set that will allow us to avoid looping in the search.
        linesAlreadyConsidered = set()
        linesAlreadyConsidered.add(self)

        frontier = queue.PriorityQueue()
        # We start by putting the direct neighbors of the line in the queue
        for line in linesConnectedToEnding:
            predecessors[line] = self
            costSoFar[line] = line.lineLength + costSoFar[self]
            # If one of the direct neighbors is a root, we can stop right there.
            if line.isARootOfTheNetwork:
                results["found"] = True
                results["distance"] = line.lineLength + self.lineLength
                return results
            else:
                frontier.put((line.uniqueID, line))

        while not frontier.empty():
            curentLine = frontier.get()
            curentLine = curentLine[1]
            linesAlreadyConsidered.add(curentLine)
            neighbors = curentLine.linesConnectedToEnding1 | curentLine.linesConnectedToEnding2

            # We will only considered the lines that have not been considered before.
            for neighbour in [neighbour for neighbour in neighbors if neighbour not in linesAlreadyConsidered]:
                # If the neighbour is the root, we can stop the search.
                if neighbour.isARootOfTheNetwork:
                    results["found"] = True
                    predecessors[neighbour] = curentLine
                    results["distance"] = neighbour.lineLength + costSoFar[curentLine]
                    return results
                # If not, we keep the search going on and add the neighbour to be looked at.
                else:
                    # But first, we gotta check if it's the best choice to put the current line as predecessor
                    # of the neighbor if he already has one.
                    if (neighbour in predecessors):
                        if costSoFar[neighbour] > neighbour.lineLength + costSoFar[curentLine]:
                            predecessors[neighbour] = curentLine
                            costSoFar[neighbour] = neighbour.lineLength + costSoFar[curentLine]

                    else:
                        predecessors[neighbour] = curentLine
                        costSoFar[neighbour] = neighbour.lineLength + costSoFar[curentLine]
                        frontier.put((costSoFar[neighbour], neighbour))

        # If we didn't found the root during the search, we return empty results.
        results["found"] = False
        return results

    def getNeighborsDownstream(self):
        """For this function to work, the network initialization must be done. Also, it must be called during
        the algorithm so as force the fact that one of the neighbors must have a smaller flux.
        The downstream neighbors are going to be the one who are connect to this line by their upstream point !"""
        downStreamNeighbors = set()
        for neighbour in self.linesConnectedDownstream:
            if self in neighbour.linesConnectedUpstream:
                downStreamNeighbors.add(neighbour)
        return downStreamNeighbors

    def getNeighborsUpstream(self):
        """For this function to work, the network initialization must be done. Also, it must be called during
        the algorithm so as force the fact that one of the neighbors must have a smaller flux.
        The downstream neighbors are going to be the one who are connect to this line by their upstream point !"""
        upStreamNeighbors = set()
        for neighbour in self.linesConnectedUpstream:
            if self in neighbour.linesConnectedDownstream:
                upStreamNeighbors.add(neighbour)
        return upStreamNeighbors

    def fluxToDownstreamNeighbors(self, feedback):
        """This function adds the flux coming from this line to the neighbours downstream"""

        if not self.isARootOfTheNetwork:

            if not self.haveFluxed:
                downStreamNeighbors = self.getNeighborsDownstream()
                for neighbour in downStreamNeighbors:
                    neighbour.flux = neighbour.flux + (self.flux / len(downStreamNeighbors))
                    # If their is a negative flux, something is really wrong. We raise an exception.
                    if neighbour.flux < 0:
                        raise QgsProcessingException("ERROR: A line flux became negative ! We were trying to add a flux of"
                                                     "quantity " + str((self.flux / len(downStreamNeighbors)))
                                                     + "from the line that is between the points " + str(self.ending1)
                                                     + " and " + str(self.ending2) + " and to the line that is between the"
                                                     " points " + str(neighbour.ending1) + " and " + str(neighbour.ending2))

        self.haveFluxed = True

    def checkIfAllUpstreamHaveFluxed(self):
        """This function checks if the upstream neighbors have all fluxed into this line so that
        we can flux this one in the right order to its own downstream neighbors."""

        # Case of the lines that are leaves
        if self.isALeafOfTheNetwork:
            return True
        # If not a leaf, we look at all of the upstream neighbors and check if they have fluxed.
        else:
            theyHaveAllFluxed = True
            upStreamNeighbors = self.getNeighborsUpstream()
            for neighbour in upStreamNeighbors:
                if not neighbour.haveFluxed:
                    theyHaveAllFluxed = False
                    break
            return theyHaveAllFluxed

    def getNeighborsUpstreamID(self):
        """Get the ID of all of the upstream neighbors lines as a string. For debugging purposes."""
        upStreamNeighbors = set()
        for neighbour in self.linesConnectedUpstream:
            if self in neighbour.linesConnectedDownstream:
                upStreamNeighbors.add(neighbour)

        return ' '.join(["{}".format(i.uniqueID) for i in upStreamNeighbors])

    def getNeighborsDownstreamID(self):
        """Get the ID of all of the downstream neighbors lines as a string. For debugging purposes."""
        downStreamNeighbors = set()
        for neighbour in self.linesConnectedDownstream:
            if self in neighbour.linesConnectedUpstream:
                downStreamNeighbors.add(neighbour)

        return ' '.join(["{}".format(i.uniqueID) for i in downStreamNeighbors])


# Methods to help the algorithm; all static, do not need to initialize an object of this class.
class WoodFluxHelper:

    @staticmethod
    def polygonsToPoints(polygonLayer, pointResolution, arbitraryDensity, isItVolume, wood_attribute_index, feedback):
        """This function transforms a given polygon layer into a set
        of QGS points according to a given resolution of points."""

        # We take all of the polygons from the polygon layer
        featuresOfLayer = list(polygonLayer.getFeatures())
        setOfPolygons = list()
        # We will need a dictionary to register the attributes of the polygons concerning wood, if the user has indicated
        # such an attribute in the polygons
        volumeForPolygonDictionary = dict()

        # First of all, for each single polygon where the wood was cut, we will calculate the wood volume contained
        # in this polygon. The volume will be calculated either on a volume attribute for the polygon, a density attribute
        # for the polygon, or an arbitrary density for all of the polygons.
        for given_feature in featuresOfLayer:
            if given_feature.hasGeometry():
                given_feature_geo = given_feature.geometry()

                # Case of multipolygons
                if given_feature_geo.wkbType() == QgsWkbTypes.MultiPolygon:
                    multi_polygon = given_feature_geo.asMultiPolygon()
                    for polygon in multi_polygon:
                        setOfPolygons.append(polygon)

                        polygonAsGeometry = QgsGeometry.fromPolygonXY(polygon)
                        multiPolygonAsGeometry = QgsGeometry.fromMultiPolygonXY(multi_polygon)
                        if wood_attribute_index is not None:
                            # If the attribute is a volume, it needs to be divided for each small polygon relative to
                            # the area of this polygon, since the attribute is for the multi-polygon.
                            if isItVolume:
                                volumeForPolygonDictionary[str(polygon)] = (polygonAsGeometry.area()/(multiPolygonAsGeometry.area())) \
                                                               * given_feature.attributes()[wood_attribute_index]
                            # If not, we calculate the volume based on the density : it is simply the density of the
                            # polygon multiplied by its area
                            else:
                                volumeForPolygonDictionary[str(polygon)] = given_feature.attributes()[wood_attribute_index] \
                                                                           * polygonAsGeometry.area()
                        else:
                            volumeForPolygonDictionary[str(polygon)] = arbitraryDensity * polygonAsGeometry.area()


                # Case of single polygons
                elif given_feature_geo.wkbType() == QgsWkbTypes.Polygon:
                    polygon = given_feature_geo.asPolygon()
                    setOfPolygons.append(polygon)
                    if wood_attribute_index is not None:
                        if isItVolume:
                            volumeForPolygonDictionary[str(polygon)] = given_feature.attributes()[wood_attribute_index]
                        else:
                            volumeForPolygonDictionary[str(polygon)] = given_feature.attributes()[wood_attribute_index] \
                                                                       * QgsGeometry.fromPolygonXY(polygon).area()
                    else:
                        volumeForPolygonDictionary[str(polygon)] = arbitraryDensity * QgsGeometry.fromPolygonXY(polygon).area()

        setOfPointsToReturn = set()
        woodVolumePerPointDictionary = dict()
        progress = 0
        feedback.setProgress(0)

        # Then, for each polygon, we will determine points that will correspond to a source
        # of wood, according to the resolution given by the user. If no point was determined, we try again with
        # a higher point resolution until some points are made inside the polygon.
        for polygon in setOfPolygons:
            setOfPointsInPolygon = set()
            resolutionCorrection = 0

            while len(setOfPointsInPolygon) == 0:
                # For that, we take the extent of the polygon layer;
                extentOfPolygon = QgsGeometry.fromPolygonXY(polygon).boundingBox()

                # Then, we loop around x and y coordinates according to the wood density parameter
                leftX = extentOfPolygon.xMinimum()
                rightX = extentOfPolygon.xMaximum()
                bottomY = extentOfPolygon.yMinimum()
                upperY = extentOfPolygon.yMaximum()

                # Then, we loop around an incrementation of x and y coordinates in the extent of the polygon,
                # with the resolution given by the user. For each set of coordinates, we check if it is a point
                # inside the polygon. If it is, we add it as a point that will be a source of wood.
                x = leftX
                while x <= rightX:
                    y = bottomY
                    while y <= upperY:
                        if feedback.isCanceled():
                            raise QgsProcessingException("ERROR: Process canceled.")
                        point = QgsPointXY(x, y)
                        if QgsGeometry.fromPolygonXY(polygon).contains(point):
                            setOfPointsToReturn.add(point)
                            setOfPointsInPolygon.add(point)
                        y += (1 / (pointResolution + resolutionCorrection))
                    x += (1 / (pointResolution + resolutionCorrection))

                resolutionCorrection += 0.5 * pointResolution

            # For each point generated in the polygon, we add it in a dictionary to know what volume of wood is
            # associated to that point. It's simply the volume for the polygon divided by the number of source points
            # generated in this polygon.
            for point in setOfPointsInPolygon:
                woodVolumePerPointDictionary[point] = volumeForPolygonDictionary[str(polygon)] \
                                                      / float(len(setOfPointsInPolygon))

            progress += 1
            feedback.setProgress(100 * progress / (len(setOfPolygons)))

        return setOfPointsToReturn, woodVolumePerPointDictionary

    @staticmethod
    def generateWoodFlux(setOfRoadsAsLines, pointsOfWoodGeneration, woodVolumePerPointDictionary, feedback):
        """This function links the points from where the wood will come from in the network
        to the lines created in the algorithm"""

        progress = 0
        feedback.setProgress(0)

        # First, we will create a dictionary that will allow us to retrieve the line to which a point
        # is associated, and the list of points for the k-d tree analysis.
        pointToLineDictionnary = dict()
        setOfPointsFromTheLines = set()
        for line in setOfRoadsAsLines:
            for pointOfLine in line.lineFeature:
                pointToLineDictionnary[pointOfLine] = line
                setOfPointsFromTheLines.add(pointOfLine)

        # Next, we create the k-d search tree with all of the points.
        numpyArrayOfRoadPoints = np.array(list(setOfPointsFromTheLines))
        spatialKDTREEForDistanceSearch = KDTree(numpyArrayOfRoadPoints, leafsize=20)

        # Then, for each point generated
        for point in pointsOfWoodGeneration:
            # We find the closest point, and put it into the dictionary to extract the closest line
            closestPointIndex = spatialKDTREEForDistanceSearch.query(point)[1]
            closestLine = pointToLineDictionnary[list(setOfPointsFromTheLines)[closestPointIndex]]
            # Then, we add a wood flux from the point to this line
            closestLine.flux += woodVolumePerPointDictionary[point]
            progress += 1
            # feedback.pushInfo("Added 1 flux to line : " + str(closestLine.uniqueID))
            feedback.setProgress(100 * progress / (len(pointsOfWoodGeneration)))

    # Function to create the fields for the attributes that we register with the lines.
    @staticmethod
    def create_fields():
        # Create an ID field to identify the lines
        id_field = QgsField("id", QVariant.Int, "integer", 10, 1)
        # Create the field containing the flux value
        flux_field = QgsField("flux", QVariant.Double, "double", 20, 1)
        # Create the field containing the list of neighbors or the line (for debugging purposes)
        upstream_neighbor_field = QgsField("upstream_neighbors", QVariant.String, "text", 100, 1)
        # Create the field containing the list of neighbors or the line (for debugging purposes)
        downstream_neighbor_field = QgsField("downstream_neighbors", QVariant.String, "text", 100, 1)
        # Then, we create a container of multiple fields
        fields = QgsFields()
        # We add the fields to the container
        fields.append(id_field)
        fields.append(flux_field)
        fields.append(upstream_neighbor_field)
        fields.append(downstream_neighbor_field)
        # We return the container with our fields.
        return fields

    # Function to create a polyline with the list of qgs.pointXY
    @staticmethod
    def create_path_feature_from_points(path_points, ID, flux, upstream_neighbors, downstream_neighbors, fields):
        # We create the geometry of the polyline
        polyline = QgsGeometry.fromPolylineXY(path_points)
        # We retrieve the fields and add them to the feature
        feature = QgsFeature(fields)
        id_index = feature.fieldNameIndex("id")
        feature.setAttribute(id_index, ID)
        flux_index = feature.fieldNameIndex("flux")
        feature.setAttribute(flux_index, flux)
        upstream_neighbor_index = feature.fieldNameIndex("upstream_neighbors")
        feature.setAttribute(upstream_neighbor_index, upstream_neighbors)
        downstream_neighbor_index = feature.fieldNameIndex("downstream_neighbors")
        feature.setAttribute(downstream_neighbor_index, downstream_neighbors)
        # We add the geometry to the feature
        feature.setGeometry(polyline)
        return feature




