# -*- 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
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 ThematicAnamorphosis(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'

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

    def createInstance(self):
        return ThematicAnamorphosis()

    def name(self):
        return 'thematic_anamorphosis'

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

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

    def groupId(self):
        return 'Analysis'

    def shortHelpString(self):
        return self.tr(
            """Generate a cartogram by iteratively moving each vertex
            along the radius that connects it to the global centroid of the layer.
            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)""")
        

    # ── 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=200,
                defaultValue=30
            )
        )

        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.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)

        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 = {}
        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

            # Reconstruction de l'index spatial sur les géométries courantes
            temp_feats = []
            for fid, geom in geoms.items():
                f = QgsFeature(fid)
                f.setGeometry(geom)
                temp_feats.append(f)
            spatial_index  = QgsSpatialIndex()
            geom_dict_iter = {}
            for f in temp_feats:
                spatial_index.insertFeature(f)
                geom_dict_iter[f.id()] = f.geometry()

            # — 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()
                densites_courantes[fid] = (
                    valeurs[fid] / aire_courante 
                )

            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)
                            dx     = sommet.x() - centre_global.x()
                            dy     = sommet.y() - centre_global.y()
                            dist_cv = math.hypot(dx, dy)

 

                            # Rayon centre → sommet
                            rayon = QgsGeometry.fromPolylineXY(
                                [centre_global, sommet]
                            )

                            # Intersections avec polygones traversés.
                            # Pour chaque segment lon_seg traversant le polygone p,
                            # la distance cible est lon_seg * √(D_moy / D_courante(p)) :
                            #   D_courante > D_moy  → scale < 1 → polygone à comprimer
                            #   D_courante < D_moy  → scale > 1 → polygone à dilater
                            candidats       = spatial_index.intersects(
                                rayon.boundingBox()
                            )
                            dist_cible      = 0.0
                            longueur_totale = 0.0

                            for cid in candidats:
                                poly_geom = geom_dict_iter[cid]
                                inter     = rayon.intersection(poly_geom)
                                if inter.isEmpty():
                                    continue
                                lon_seg = float(inter.length())

                                d_c   = densites_courantes[cid]
                                scale = (float(d_c)/float(densite_moy)) **0.5
                                dist_cible      += lon_seg * scale
                                longueur_totale += lon_seg



                            # Déplacement radial amorti + garde-fou anti-retournement
                            facteur   = 1.0+damping * ((dist_cible / longueur_totale)-1.0)
                            
                            nouveau_x = centre_global.x() + dx * facteur
                            nouveau_y = centre_global.y() + dy * facteur

                            deplacement = math.hypot(
                                nouveau_x - sommet.x(),
                                nouveau_y - sommet.y()
                            )
                            deplacement_max = max(deplacement_max, deplacement)

                            nouveaux_pts.append(QgsPointXY(nouveau_x, nouveau_y))

                        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

            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()}
