import math
from datetime import datetime, timedelta
from collections import defaultdict

from qgis.PyQt.QtCore import QVariant
from qgis.core import (
    QgsFeature, QgsField, QgsFields, QgsGeometry, QgsRectangle,
    QgsProcessing, QgsProcessingAlgorithm,
    QgsProcessingParameterFeatureSource, QgsProcessingParameterString,
    QgsProcessingParameterFeatureSink, QgsProcessingParameterExtent, 
    QgsWkbTypes, QgsFeatureRequest, QgsFeatureSink
)

from .common import DEFAULT_REF_DATE, parse_ref_date

class HovmollerStatsExtractor(QgsProcessingAlgorithm):
    INPUT_GRID = 'INPUT_GRID'
    AREA_OF_INTEREST = 'AREA_OF_INTEREST'
    INPUT_RAW = 'INPUT_RAW'
    INPUT_SHORELINES = 'INPUT_SHORELINES'
    DATE_FIELD = 'DATE_FIELD'
    REF_DATE = 'REF_DATE'

    OUTPUT_PTS = 'OUTPUT_PTS'
    OUTPUT_LINES = 'OUTPUT_LINES'

    def name(self): return 'hovmoller_step3_stats_extractor'
    def displayName(self): return '2. Interactive: Extract Stats & Zoom'
    def group(self): return 'Shoreline Analysis'
    def groupId(self): return 'shoreline'
    def createInstance(self): return HovmollerStatsExtractor()

    def shortHelpString(self):
        return (
            "<h3>Instructions</h3>"
            "<ol>"
            "<li><b>MANDATORY:</b> In the Layers panel on the left, click and highlight (make active) the <b>Grid Points</b> layer.</li>"
            "<li>For the second parameter, click <b>[...] -> 'Select Extent on Canvas'</b> and draw a rectangle over your area of interest on the Hovmöller diagram.</li>"
            "<li>Specify the other required layers (Raw Intersection Points and Clipped Shorelines).</li>"
            "<li>Click Run. Upon completion, QGIS will automatically calculate all statistics and zoom directly to the real beach area.</li>"
            "</ol>"
        )

    def initAlgorithm(self, config=None):
        # 1. Grid layer selection with warning
        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT_GRID, '1. Grid Points Layer (WARNING: Make this layer active in the panel!)', [QgsProcessing.TypeVectorPoint]))
        
        # 2. Extent selection on canvas
        self.addParameter(QgsProcessingParameterExtent(self.AREA_OF_INTEREST, '2. Select area on the diagram (Click [...] -> Select Extent on Canvas)'))
        
        # Other parameters
        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT_RAW, '3. All Raw Intersection Points', [QgsProcessing.TypeVectorPoint]))
        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT_SHORELINES, '4. All Clipped Shorelines', [QgsProcessing.TypeVectorLine]))
        self.addParameter(QgsProcessingParameterString(self.DATE_FIELD, 'Date Field (e.g., date_clean)', defaultValue='date_clean'))
        self.addParameter(QgsProcessingParameterString(self.REF_DATE, 'Reference Date (epoch for day index, YYYY-MM-DD)', defaultValue=DEFAULT_REF_DATE))

        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT_PTS, 'Extracted Points with Stats'))
        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT_LINES, 'Extracted Shorelines'))

    def processAlgorithm(self, parameters, context, feedback):
        bbox = self.parameterAsExtent(parameters, self.AREA_OF_INTEREST, context)
        bbox_geom = QgsGeometry.fromRect(bbox)

        grid_source = self.parameterAsSource(parameters, self.INPUT_GRID, context)
        raw_source = self.parameterAsSource(parameters, self.INPUT_RAW, context)
        lines_source = self.parameterAsSource(parameters, self.INPUT_SHORELINES, context)
        date_f = self.parameterAsString(parameters, self.DATE_FIELD, context)
        ref_date = parse_ref_date(self.parameterAsString(parameters, self.REF_DATE, context))

        # 1. Collect unique longi and days from the selected extent
        selected_longi = set()
        selected_days = set()
        for f in grid_source.getFeatures():
            if f.geometry().intersects(bbox_geom):
                selected_longi.add(f['longi'])
                selected_days.add(f['days'])

        if not selected_longi:
            raise Exception("No grid points found in the selected area. Please check your selection or ensure the correct layer is active.")

        # 2. Filter raw points
        raw_pts_by_longi = defaultdict(list)
        original_raw_fields = raw_source.fields()
        
        for f in raw_source.getFeatures():
            l, d = f['longi'], f['days']
            if l in selected_longi and d in selected_days:
                raw_pts_by_longi[l].append({
                    'days': d,
                    'dist': f['dist_m'],
                    'geom': QgsGeometry(f.geometry()),
                    'attrs': f.attributes()
                })

        # 3. Prepare output fields (Adding 'real_date' and Stats)
        out_fields = QgsFields(original_raw_fields)
        out_fields.append(QgsField('real_date', QVariant.String)) # New Field for human-readable date
        
        stats_columns = ['SCE', 'NSM', 'EPR', 'MSP', 'SD', 'CV', 'LRR', 'LSE', 'LCI', 'R2']
        for stat in stats_columns:
            out_fields.append(QgsField(stat, QVariant.Double))

        (pts_sink, pts_id) = self.parameterAsSink(parameters, self.OUTPUT_PTS, context, out_fields, QgsWkbTypes.Point, raw_source.sourceCrs())

        self.zoom_extent = QgsRectangle()
        first_point_for_extent = True

        # 4. Calculate statistics
        for l, pts in raw_pts_by_longi.items():
            pts.sort(key=lambda x: x['days'])
            
            dists = [p['dist'] for p in pts]
            days = [p['days'] for p in pts]
            n = len(pts)

            sce = nsm = epr = msp = sd = cv = lrr = lse_rate = lci = r2 = 0.0

            if n > 0:
                sce = max(dists) - min(dists)
                nsm = dists[-1] - dists[0]
                msp = sum(dists) / n
                
                if n > 1:
                    variance = sum((d - msp)**2 for d in dists) / (n - 1)
                    sd = math.sqrt(variance)
                    if msp != 0:
                        cv = (sd / abs(msp)) * 100

                    years = [(d - days[0]) / 365.25 for d in days]
                    total_years = years[-1]
                    if total_years > 0:
                        epr = nsm / total_years

                    mean_x = sum(years) / n
                    ss_xx = sum((x - mean_x)**2 for x in years)
                    ss_yy = sum((y - msp)**2 for y in dists)
                    ss_xy = sum((x - mean_x)*(y - msp) for x, y in zip(years, dists))

                    if ss_xx != 0:
                        lrr = ss_xy / ss_xx
                        intercept = msp - lrr * mean_x
                        
                        if n > 2:
                            sse = sum((y - (lrr * x + intercept))**2 for x, y in zip(years, dists))
                            lse_estimate = math.sqrt(sse / (n - 2))
                            lse_rate = lse_estimate / math.sqrt(ss_xx)
                            
                            if (ss_xx * ss_yy) != 0:
                                r2 = (ss_xy**2) / (ss_xx * ss_yy)
                            
                            df = n - 2
                            t_val = self.get_t_value(df)
                            lci = t_val * lse_rate

            for p in pts:
                new_f = QgsFeature(out_fields)
                new_f.setGeometry(p['geom'])
                
                # Convert 'days' integer back to 'YYYY-MM-DD' date string
                real_date_str = (ref_date + timedelta(days=p['days'])).strftime('%Y-%m-%d')
                
                attrs = p['attrs'] + [real_date_str, sce, nsm, epr, msp, sd, cv, lrr, lse_rate, lci, r2]
                new_f.setAttributes(attrs)
                pts_sink.addFeature(new_f, QgsFeatureSink.FastInsert)
                
                # Expand Bounding Box for automatic zoom
                if first_point_for_extent:
                    self.zoom_extent = p['geom'].boundingBox()
                    first_point_for_extent = False
                else:
                    self.zoom_extent.combineExtentWith(p['geom'].boundingBox())

        # 5. Extract Shorelines
        (lines_sink, lines_id) = self.parameterAsSink(parameters, self.OUTPUT_LINES, context, lines_source.fields(), QgsWkbTypes.LineString, lines_source.sourceCrs())
        
        skipped = 0
        for f in lines_source.getFeatures():
            try:
                dt = datetime.strptime(str(f[date_f]).strip(), '%Y%m%d')
            except (ValueError, TypeError, KeyError):  # unparseable date value
                skipped += 1
                continue
            day_idx = (dt - ref_date).days
            if day_idx in selected_days:
                lines_sink.addFeature(f, QgsFeatureSink.FastInsert)
        if skipped:
            feedback.pushWarning(
                f"Skipped {skipped} shoreline feature(s) with unparseable '{date_f}' dates.")

        return {
            self.OUTPUT_PTS: pts_id,
            self.OUTPUT_LINES: lines_id
        }

    def postProcessAlgorithm(self, context, feedback):
        """ Automatically zooms map canvas to the processed area upon completion. """
        try:
            import qgis.utils
            if qgis.utils.iface and hasattr(self, 'zoom_extent') and not self.zoom_extent.isEmpty():
                self.zoom_extent.scale(1.2)
                qgis.utils.iface.mapCanvas().setExtent(self.zoom_extent)
                qgis.utils.iface.mapCanvas().refresh()
        except Exception as e:
            feedback.pushInfo(f"Automatic zoom failed: {str(e)}")
            
        return super().postProcessAlgorithm(context, feedback)

    def get_t_value(self, df):
        # Student's t-distribution table for 95% confidence interval (two-tailed)
        t_table = {
            1: 12.706, 2: 4.303, 3: 3.182, 4: 2.776, 5: 2.571,
            6: 2.447, 7: 2.365, 8: 2.306, 9: 2.262, 10: 2.228,
            11: 2.201, 12: 2.179, 13: 2.160, 14: 2.145, 15: 2.131,
            16: 2.120, 17: 2.110, 18: 2.101, 19: 2.093, 20: 2.086,
            21: 2.080, 22: 2.074, 23: 2.069, 24: 2.064, 25: 2.060,
            26: 2.056, 27: 2.052, 28: 2.048, 29: 2.045, 30: 2.042
        }
        if df in t_table: return t_table[df]
        if df > 30: return 1.96 # Approximation for large samples
        return 0.0