# -*- coding: utf-8 -*-
"""
Algorithme de traitement QGIS — Anamorphose thématique (cartogramme radial)
Déplacement itératif des sommets le long de rayons issus du centroïde global,
avec pondération en √(D_cible / D_moy) pour respecter la proportionnalité des aires.

Installation : copier ce fichier dans un dossier accessible depuis l'éditeur de
scripts de la Boîte à outils de traitement, ou le charger via
    iface.processingToolbox().loadAlgorithm(CartogrammeRadialAlgorithm())
"""

from qgis.core import (
    QgsProcessing,
    QgsProcessingAlgorithm,
    QgsProcessingParameterFeatureSource,
    QgsProcessingParameterField,
    QgsProcessingParameterFileDestination,
    QgsProcessingParameterNumber,
    QgsProcessingOutputVectorLayer,
    QgsProcessingException,
    QgsFeature,
    QgsGeometry,
    QgsPointXY,
    QgsVectorLayer,
    QgsVectorFileWriter,
    QgsSpatialIndex,
    QgsWkbTypes,
    QgsProject,
    QgsCoordinateTransformContext,
)
from qgis.PyQt.QtCore import QCoreApplication
import math, copy
from collections import defaultdict


# Précision d'arrondi pour l'identification des sommets partagés (nb de décimales).
# Adapter si les coordonnées sont en degrés (SNAP_DIGITS=7) ou en mètres (SNAP_DIGITS=2).
SNAP_DIGITS = 2      # Lambert-93 / coordonnées métriques → 2 décimales (cm)
                     # Coordonnées géographiques (degrés) → 7
MAX_LOG     = 0.4    # borne de |ln(f_pas)| ≈ facteur max 1.49 par itération



def _snap(x, y):
    """Clé hashable représentant un sommet arrondi."""
    return (round(x, SNAP_DIGITS), round(y, SNAP_DIGITS))



def _geom_to_lists(geom):
    if geom.wkbType() == QgsWkbTypes.MultiPolygon:
        return [[[QgsPointXY(pt) for pt in ann]
                 for ann in part]
                for part in geom.asMultiPolygon()], False
    return [[[QgsPointXY(pt) for pt in ann]
             for ann in geom.asPolygon()]], True


def _lists_to_geom(parties, mono):
    if mono:
        return QgsGeometry.fromPolygonXY(parties[0])
    return QgsGeometry.fromMultiPolygonXY(parties)


def _aire(parties):
    """Aire signée par Shoelace (trous en négatif)."""
    total = 0.0
    for partie in parties:
        for ia, anneau in enumerate(partie):
            s = 0.0
            n = len(anneau)
            for i in range(n):
                j = (i + 1) % n
                s += anneau[i].x() * anneau[j].y()
                s -= anneau[j].x() * anneau[i].y()
            a = abs(s) / 2.0
            total += a if ia == 0 else -a
    return total


def _barycentre(parties):
    """Barycentre simple de tous les sommets."""
    pts = [pt for part in parties for ann in part for pt in ann]
    return (sum(p.x() for p in pts) / len(pts),
            sum(p.y() for p in pts) / len(pts))


class GravityAnamorphosis(QgsProcessingAlgorithm):
    """
    Anamorphose thématique par déplacement radial itératif des sommets.
    La surface des polygones en sortie est proportionnelle au champ attributaire choisi.
    """

    # ── Identifiants des paramètres ───────────────────────────────────────────
    INPUT          = 'INPUT'
    FIELD          = 'FIELD'
    OUTPUT_GPKG    = 'OUTPUT_GPKG'
    N_ITERATIONS   = 'N_ITERATIONS'
    SEUIL          = 'THRESHOLD'
    DAMPING        = 'DAMPING'
    OUTPUT_LAYER   = 'OUTPUT_LAYER'
    RESISTANCE  ='RESISTANCE'
    LIMIT='LIMIT'

    def tr(self, string,context=''):
        if context == '':
            context = self.__class__.__name__
        return QCoreApplication.translate(context, string)

    def createInstance(self):
        return GravityAnamorphosis()

    def name(self):
        return 'gravity_anamorphosis'

    def displayName(self):
        return self.tr('Gravity anamorphosis')

    def group(self):
        return self.tr('Analysis')

    def groupId(self):
        return 'Analysis'

    def shortHelpString(self):
        return self.tr(
            """Generate a anamorphosis by iteratively moving each vertex
            along the radius from other polygons centroid weighted by distance
            The weighting uses √(D_target / D_mean) so that the final areas are proportional to the chosen attribute value.
            
            Parameters:
            • Source layer — polygons to be deformed
            • Value field — numeric field to represent
            • Output file — path to a .gpkg
            • Max iterations — 20 is often sufficient
            • Convergence threshold — maximum acceptable displacement (in layer units)
            • Damping factor — damping 0 < d ≤ 1 (0.4 recommended))
            • Resistance parameter - Calibrate the weighting (higher value give more inlfuence of further polygons)
            • Power parameter — Power of the weighting function""")
            
        

    # ── Définition des paramètres ─────────────────────────────────────────────
    def initAlgorithm(self, config=None):

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr('Source layer'),
                [QgsProcessing.TypeVectorPolygon]
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.FIELD,
                self.tr('Value field'),
                parentLayerParameterName=self.INPUT,
                type=QgsProcessingParameterField.Numeric
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.OUTPUT_GPKG,
                self.tr('Output File'),
                fileFilter='GeoPackage (*.gpkg)'
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.N_ITERATIONS,
                self.tr('Max number of iterations'),
                type=QgsProcessingParameterNumber.Integer,
                minValue=1,
                maxValue=5000,
                defaultValue=300
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.SEUIL,
                self.tr('Convergence threshold (layer units)'),
                type=QgsProcessingParameterNumber.Double,
                minValue=0.0001,
                defaultValue=1.0
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.DAMPING,
                self.tr('Damping factor (0 < d ≤ 1)'),
                type=QgsProcessingParameterNumber.Double,
                minValue=0.01,
                maxValue=1.0,
                defaultValue=0.4
            )
        )
        
        self.addParameter(
            QgsProcessingParameterNumber(
                self.RESISTANCE,
                self.tr('Resistance parameter'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=10
            )
        )
        
        self.addParameter(
            QgsProcessingParameterNumber(
                self.LIMIT,
                self.tr('Power parameter'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=2
            )
        )

        self.addOutput(
            QgsProcessingOutputVectorLayer(
                self.OUTPUT_LAYER,
                self.tr('Output layer')
            )
        )

    # ── Algorithme principal ──────────────────────────────────────────────────
    def processAlgorithm(self, parameters, context, feedback):

        # — Récupération des paramètres —
        source       = self.parameterAsSource(parameters, self.INPUT, context)
        field_name   = self.parameterAsString(parameters, self.FIELD, context)
        gpkg_path    = self.parameterAsFileOutput(parameters, self.OUTPUT_GPKG, context)
        n_iter       = self.parameterAsInt(parameters, self.N_ITERATIONS, context)
        seuil        = self.parameterAsDouble(parameters, self.SEUIL, context)
        damping      = self.parameterAsDouble(parameters, self.DAMPING, context)
        resistance   = self.parameterAsDouble(parameters, self.RESISTANCE, context)
        limite   = self.parameterAsDouble(parameters, self.LIMIT, context)
        if source is None:
            raise QgsProcessingException(self.tr('Invalid layer source'))

        features = list(source.getFeatures())
        total    = len(features)
        if total == 0:
            raise QgsProcessingException(self.tr('The layer is empty'))

        feedback.pushInfo(f"{total} {self.tr('loaded polygons.')}")

        # — Centroïde global —
        extent         = source.sourceExtent()
        centre_global  = QgsPointXY(extent.center())
        feedback.pushInfo(
            f"{self.tr('Global centroid')} : ({centre_global.x():.2f}, {centre_global.y():.2f})"
        )

        # — Valeurs attributaires fixes (ne changent pas au fil des itérations) —
        valeurs = {}
        centres = {} 
        centres_save={}
        for feat in features:
            valeur = feat[field_name]
            valeurs[feat.id()] = float(valeur) if valeur is not None else 0.0

        # — Densité moyenne globale : ΣValeurs / ΣAires (invariante) —
        somme_valeurs = sum(valeurs.values())
        somme_aires   = sum(feat.geometry().area() for feat in features)
        if somme_aires <= 0 or somme_valeurs <= 0:
            raise QgsProcessingException(
                self.tr('the sum of the areas is null - check layer source.')
            )
        densite_moy = somme_valeurs / somme_aires
        feedback.pushInfo(f"{self.tr('Densité moyenne globale')} : {densite_moy:.6f}")

        # — Copie de travail des géométries —
        geoms = {feat.id(): QgsGeometry(feat.geometry()) for feat in features}

        
        # ── Boucle d'itérations ──────────────────────────────────────────────
        for iteration in range(n_iter):

            if feedback.isCanceled():
                break


            # — Recalcul des densités courantes sur les géométries déformées —
            # Valeur attributaire fixe / aire courante (qui évolue à chaque itération)
            densites_courantes = {}
            for fid, geom in geoms.items():
                aire_courante = geom.area()
                #feedback.pushInfo(f"{self.tr('Itération')} {valeurs[fid]:3f} — {self.tr('maximum offset')} = {aire_courante:.4f}")
                    
                densites_courantes[fid] = (valeurs[fid] / aire_courante )
                #centres[fid]=geom.poleOfInaccessibility(1.0)[0].asPoint()
                centres[fid]=geom.boundingBox().center()
                

            deplacement_max = 0.0
            nouveaux_geoms  = {}

            # Progression dans la barre QGIS
            feedback.setProgress(int(iteration / n_iter * 100))

            for fid, geom in geoms.items():

                if geom.wkbType() not in (
                    QgsWkbTypes.Polygon,
                    QgsWkbTypes.MultiPolygon
                ):
                    nouveaux_geoms[fid] = geom
                    continue

                # Gestion mono / multi-polygone
                if geom.wkbType() == QgsWkbTypes.MultiPolygon:
                    parties = geom.asMultiPolygon()
                    mono    = False
                else:
                    parties = [geom.asPolygon()]
                    mono    = True

                nouvelles_parties = []

                for partie in parties:
                    nouveaux_anneaux = []

                    for anneau in partie:
                        nouveaux_pts = []

                        for pt in anneau:
                            sommet = QgsPointXY(pt)
                            nouveau_x=0
                            nouveau_y=0

 

                            # Rayon centre → sommet

                            facteur=0
                            dist_poids=0
                            '''geoms_copy=cloned = {
                                k: v
                                for k, v in geoms.items()}'''

                            for cid, gomme in geoms.items():
                                #if math.hypot(centres_save[fid].x()-centres_save[cid].x(),centres_save[fid].y()-centres_save[cid].y())<limite*1000:
            
                            
                                
                                
                                dx     = sommet.x() - centres[cid].x()
                                dy     = sommet.y() - centres[cid].y()
                                dist_cv = math.hypot(dx, dy)*0.001
                                
                                

                                
                                d_c   = densites_courantes[cid]
                                scale = ((float(d_c)/float(densite_moy)) **0.5)-1
                                dist_poids+= scale * (2**(-(float(dist_cv)/float(resistance))**limite))


                                # Déplacement radial amorti + garde-fou anti-retournement
                                facteur += (2**(-(float(dist_cv)/float(resistance))**limite))
                                #feedback.pushInfo(f"{dist_poids:.4f} {facteur:.4f}")

                    
                            nouveau_x +=  float(dx) * dist_poids*damping
                            nouveau_y +=  float(dy) * dist_poids*damping

                            deplacement_max=max(deplacement_max, math.hypot(nouveau_x,nouveau_y))

                            if facteur==0:
                                feedback.pushInfo(self.tr("Increase the parameter value"))
                                break
                            nouveaux_pts.append(QgsPointXY(sommet.x()+nouveau_x/facteur, sommet.y()+nouveau_y/facteur))

                        nouveaux_anneaux.append(nouveaux_pts)
                    nouvelles_parties.append(nouveaux_anneaux)

                if mono:
                    nouveaux_geoms[fid] = QgsGeometry.fromPolygonXY(
                        nouvelles_parties[0]
                    ).makeValid()
                else:
                    nouveaux_geoms[fid] = QgsGeometry.fromMultiPolygonXY(
                        nouvelles_parties
                    ).makeValid()

            geoms = nouveaux_geoms
            #deplacement_max=1-(somme_valeurs/somme_aires)

            feedback.pushInfo(
                f"{self.tr('Itération')} {iteration+1:3d} — {self.tr('maximum offset')} = {deplacement_max:.4f}"
            )

            if deplacement_max < seuil:
                feedback.pushInfo(
                    f"✓ {self.tr('Convergence atteinte à l\'itération')} {iteration+1}."
                )
                break

        # ── Construction de la couche résultat ──────────────────────────────
        result_layer = QgsVectorLayer(
            f"Polygon?crs={source.sourceCrs().authid()}",
            "anamorphosis",
            "memory"
        )
        pr = result_layer.dataProvider()
        pr.addAttributes(source.fields())
        result_layer.updateFields()

        feat_map = {f.id(): f for f in features}
        result_feats = []
        for fid, new_geom in geoms.items():
            new_feat = QgsFeature(feat_map[fid])
            new_feat.setGeometry(new_geom)
            result_feats.append(new_feat)
        pr.addFeatures(result_feats)
        result_layer.updateExtents()

        # ── Écriture du GeoPackage ───────────────────────────────────────────
        save_options              = QgsVectorFileWriter.SaveVectorOptions()
        save_options.driverName   = 'GPKG'
        save_options.fileEncoding = 'UTF-8'
        save_options.layerName    = 'anamorphosis'

        # writeAsVectorFormatV3 retourne un tuple de 2 ou 3 éléments selon la
        # version de QGIS — on indexe directement pour être compatible partout.
        retval       = QgsVectorFileWriter.writeAsVectorFormatV3(
            result_layer,
            gpkg_path,
            QgsCoordinateTransformContext(),
            save_options
        )
        write_result = retval[0]
        error_msg    = retval[1]

        if write_result != QgsVectorFileWriter.NoError:
            raise QgsProcessingException(
                self.tr(f'Errors writing GPKG : {error_msg}')
            )

        feedback.pushInfo(f"GeoPackage  : {gpkg_path}")

        # Chargement de la couche depuis le fichier pour le retour
        output_layer = QgsVectorLayer(
            f"{gpkg_path}|layername=anamorphosis",
            "anamorphosis",
            "ogr"
        )
        QgsProject.instance().addMapLayer(output_layer)

        return {self.OUTPUT_LAYER: output_layer.id()}
