# Copyright (C) 2026 Wilhelmus Dingenus de Pooter
# SPDX-License-Identifier: GPL-3.0-or-later
"""
sourceval/reconcile.py — the reconciliation engine.
================================================================================

Given the ENC point features (``ChartFeature``), the source features
(``SourceFeature``) and a ``Recipe``, run the four checks and return a ``Report``
(the Auditor's report type, so HTML/CSV/flags-layer all work unchanged).

Spatial match (v1): two points are the "same feature" if within ``recipe.buffer_m``
metres. Distance is the haversine great-circle metre distance — at a 0–5 m buffer
its error versus a full ellipsoidal computation is far below the tolerance, and
it needs no QGIS. (Inside QGIS the GUI can substitute ``QgsDistanceArea`` for an
exact ellipsoidal figure; the engine interface is the same.)

Matching is brute force O(chart × source). Point layers are small enough that
this is fine for v1; a spatial index is a later optimisation if needed.

Comparison rule (deliberately low-noise for v1):
  * A mismatch is only raised when *both* sides carry a value and they disagree
    after the recipe's value mapping. A blank chart attribute is the Auditor's
    business (missing mandatory data), not the Source Validator's.
  * If the recipe's value map does not cover the chart value, the pair is skipped
    (an incomplete recipe is not a data finding).
  * The per-mapping ``unmapped`` policy decides whether a source value the recipe
    never maps to is flagged or ignored.
"""

import math
from datetime import datetime, timezone

try:                                   # QGIS package import
    from ..validation.engine import Report, Finding
except ImportError:                    # headless harness
    from validation.engine import Report, Finding
from . import checks

_EARTH_R = 6371008.8     # IUGG mean radius, metres


def distance_m(lon1, lat1, lon2, lat2):
    """Great-circle distance between two lon/lat points, in metres (haversine)."""
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dphi = math.radians(lat2 - lat1)
    dlmb = math.radians(lon2 - lon1)
    a = (math.sin(dphi / 2) ** 2
         + math.cos(p1) * math.cos(p2) * math.sin(dlmb / 2) ** 2)
    return 2 * _EARTH_R * math.asin(min(1.0, math.sqrt(a)))


def _norm(v):
    """Normalise a value for comparison: numbers compare numerically, everything
    else as a case-folded, stripped string."""
    if v is None:
        return None
    if isinstance(v, bool):
        return str(v)
    if isinstance(v, (int, float)):
        return float(v)
    s = str(v).strip()
    try:
        return float(s)
    except ValueError:
        return s.casefold()


def _nearest(chart_f, sources, buffer_m):
    """The nearest source feature to ``chart_f`` within ``buffer_m`` (or None)."""
    best, best_d = None, None
    for s in sources:
        d = distance_m(chart_f.lon, chart_f.lat, s.lon, s.lat)
        if d <= buffer_m and (best_d is None or d < best_d):
            best, best_d = s, d
    return best


def _any_within(lon, lat, others, buffer_m):
    for o in others:
        if distance_m(lon, lat, o.lon, o.lat) <= buffer_m:
            return True
    return False


def _expected_source_value(mapping, chart_raw, catalogue):
    """The source value the recipe expects for a given ENC value, plus a label
    for messages. Returns (expected, enc_label, covered)."""
    if mapping.enumerated:
        enc_label = catalogue.label_for(mapping.enc_attr, chart_raw)
        if enc_label not in mapping.value_map:
            return None, enc_label, False          # recipe doesn't cover it
        return mapping.value_map[enc_label], enc_label, True
    # free-text / numeric
    key = str(chart_raw)
    if mapping.value_map:
        if key not in mapping.value_map:
            return None, key, False
        return mapping.value_map[key], key, True
    return chart_raw, key, True                     # direct comparison


def _source_value_known(mapping, source_val):
    """Is this source value one the recipe maps to from some ENC value?"""
    target = _norm(source_val)
    return any(_norm(v) == target for v in mapping.value_map.values())


def _compare_attrs(chart_f, source_f, recipe, catalogue):
    """Yield (enc_attr, detail) for each mapped attribute that disagrees."""
    for m in recipe.mappings:
        chart_raw = chart_f.attrs.get(m.enc_attr)
        if chart_raw is None:
            continue                                # blank chart attr → not ours
        source_val = source_f.attrs.get(m.source_field)
        if m.is_null(source_val):
            continue                                # source unspecified / NULL
            # sentinel (e.g. "L") → skip
        expected, enc_label, covered = _expected_source_value(
            m, chart_raw, catalogue)
        if not covered:
            continue                                # recipe gap, not a finding
        if _norm(expected) == _norm(source_val):
            continue                                # match
        if (m.unmapped == "ignore"
                and not _source_value_known(m, source_val)):
            continue                                # user chose to ignore these
        shown_expected = "blank" if expected is None else expected
        detail = (f"{m.enc_attr}={enc_label!r} expects "
                  f"{m.source_field}={shown_expected!r}, source has "
                  f"{source_val!r}")
        yield m.enc_attr, detail


def _chart_area(chart_features, margin_m):
    """Bounding box of the chart features, padded by margin_m (converted to
    degrees at the box's mean latitude). Returns (minlon, minlat, maxlon, maxlat)
    or None if there are no features."""
    if not chart_features:
        return None
    lons = [c.lon for c in chart_features]
    lats = [c.lat for c in chart_features]
    minlon, maxlon = min(lons), max(lons)
    minlat, maxlat = min(lats), max(lats)
    mean_lat = (minlat + maxlat) / 2.0
    dlat = margin_m / 111320.0
    dlon = margin_m / (111320.0 * max(0.1, math.cos(math.radians(mean_lat))))
    return (minlon - dlon, minlat - dlat, maxlon + dlon, maxlat + dlat)


def _in_ring(x, y, ring):
    """Ray-casting point-in-polygon for a single ring of (lon, lat) vertices."""
    inside = False
    n = len(ring)
    if n < 3:
        return False
    j = n - 1
    for i in range(n):
        xi, yi = ring[i]
        xj, yj = ring[j]
        if ((yi > y) != (yj > y)) and \
                (x < (xj - xi) * (y - yi) / (yj - yi + 1e-30) + xi):
            inside = not inside
        j = i
    return inside


def _ring_bbox(ring):
    xs = [p[0] for p in ring]
    ys = [p[1] for p in ring]
    return (min(xs), min(ys), max(xs), max(ys))


def point_in_coverage(x, y, coverage):
    """Is (lon, lat) inside the coverage? ``coverage`` is a list of polygons,
    each ``(exterior_ring, [hole_ring, ...], bbox)`` with rings as lon/lat
    vertex lists. A point is covered if inside an exterior ring and outside that
    polygon's holes. The bbox is a cheap pre-reject."""
    for ext, holes, bbox in coverage:
        if not (bbox[0] <= x <= bbox[2] and bbox[1] <= y <= bbox[3]):
            continue
        if _in_ring(x, y, ext) and not any(_in_ring(x, y, h) for h in holes):
            return True
    return False


def prepare_coverage(polygons):
    """Attach a bbox to each (exterior, holes) polygon for fast rejection."""
    return [(ext, holes, _ring_bbox(ext)) for ext, holes in polygons if ext]


def reconcile(chart_features, source_features, recipe, catalogue, coverage=None):
    """Run the four checks and return a populated ``Report``.

    ``coverage`` (optional) is a prepared coverage polygon set (see
    ``prepare_coverage``) — the real chart footprint (M_COVR). When given, the
    "in source, not in chart" check is clipped to it by point-in-polygon, which
    is exact (no false flags in the bounding-box corners). Without it, the check
    falls back to the padded chart bounding box."""
    report = Report()
    report.ran = list(checks.ALL)
    buf = recipe.buffer_m
    klass = recipe.enc_class

    # --- Known-fault exclusions ---------------------------------------------
    # Remove excluded features from BOTH populations up front (before any test),
    # so an excluded beacon can't anti-match a nearby source buoy either, and
    # record the counts so the report can disclose exactly what was withheld —
    # an exclusion you can't see is indistinguishable from a dropped-feature bug.
    ex_classes = recipe.excluded_class_set()
    ex_lnam = recipe.excluded_lnam_set()
    ex_source = recipe.excluded_source_set()
    report.excluded_classes = sorted(ex_classes)
    report.enc_excluded_classes_count = 0
    report.enc_excluded_features_count = 0
    report.source_excluded_count = 0
    # The flagged known faults themselves (LNAM/OBJNAM, source id), so the
    # exported report can list exactly which features were withheld.
    report.excluded_enc_list = [dict(e) for e in recipe.enc_exclude]
    report.excluded_source_list = [dict(e) for e in recipe.source_exclude]
    if ex_classes or ex_lnam:
        kept = []
        for c in chart_features:
            if (c.klass or "").upper() in ex_classes:
                report.enc_excluded_classes_count += 1
                continue
            ln = c.lnam()
            if ln and ln.upper() in ex_lnam:
                report.enc_excluded_features_count += 1
                continue
            kept.append(c)
        chart_features = kept
    if ex_source:
        kept = []
        for s in source_features:
            if str(getattr(s, "ident", "")) in ex_source:
                report.source_excluded_count += 1
                continue
            kept.append(s)
        source_features = kept

    # Clip the "in source, not in chart" check to the chart's real area, so a
    # national register doesn't report out-of-area features as omissions.
    # Coverage polygon (exact) if supplied, else the padded bounding box.
    use_coverage = bool(coverage)
    area = None
    out_of_area = 0
    if not use_coverage and getattr(recipe, "restrict_to_area", True):
        area = _chart_area(chart_features, getattr(recipe, "area_margin_m", 250.0))
    report.clip_mode = ("coverage" if use_coverage
                        else "bbox" if area is not None else "none")

    def in_area(s):
        if use_coverage:
            return point_in_coverage(s.lon, s.lat, coverage)
        if area is None:
            return True
        return area[0] <= s.lon <= area[2] and area[1] <= s.lat <= area[3]

    # --- Check 1: in source, not in chart -----------------------------------
    for s in source_features:
        if not in_area(s):
            out_of_area += 1
            continue
        if not _any_within(s.lon, s.lat, chart_features, buf):
            f = Finding(
                checks.IN_SOURCE_NOT_CHART, "(source)", klass, s.ident,
                detail=f"source id {s.ident} @ {s.lat:.5f}, {s.lon:.5f}",
                geom_wkt=f"POINT ({s.lon:.7f} {s.lat:.7f})")
            f.exclude = ("source", str(s.ident), "")
            report.add(f)
    report.source_out_of_area = out_of_area

    # --- Checks 2 & 3: chart features vs their nearest source ---------------
    for c in chart_features:
        enc_excl = ("enc", c.lnam() or "", c.attrs.get("OBJNAM") or "",
                    c.klass or "")
        match = _nearest(c, source_features, buf)
        if match is None:
            name = c.attrs.get("OBJNAM") or ""
            where = f"@ {c.lat:.5f}, {c.lon:.5f}"
            det = f"{name} {where}".strip() if name else f"fid {c.fid} {where}"
            f = Finding(
                checks.IN_CHART_NOT_SOURCE, c.cell, klass, c.fid,
                detail=det, geom_wkt=c.wkt())
            f.exclude = enc_excl
            report.add(f)
            continue
        for _attr, detail in _compare_attrs(c, match, recipe, catalogue):
            f = Finding(
                checks.VALUE_MISMATCH, c.cell, klass, c.fid,
                detail=detail, geom_wkt=c.wkt())
            f.exclude = enc_excl
            report.add(f)

    # --- Check 4: same identifier, different location -----------------------
    if recipe.identity:
        enc_key = recipe.identity.get("enc")
        src_key = recipe.identity.get("source")
        if enc_key and src_key:
            src_by_id = {}
            for s in source_features:
                ident = s.attrs.get(src_key)
                if ident is not None and str(ident).strip() != "":
                    src_by_id.setdefault(str(ident), []).append(s)
            for c in chart_features:
                ident = c.identity(enc_key)
                if ident is None:
                    continue
                peers = src_by_id.get(str(ident))
                if not peers:
                    continue
                if not _any_within(c.lon, c.lat, peers, buf):
                    near = min(peers, key=lambda s: distance_m(
                        c.lon, c.lat, s.lon, s.lat))
                    d = distance_m(c.lon, c.lat, near.lon, near.lat)
                    f = Finding(
                        checks.IDENTITY_DISPLACED, c.cell, klass, c.fid,
                        detail=f"identity {ident!r} differs by {d:.0f} m",
                        geom_wkt=c.wkt())
                    f.exclude = ("enc", c.lnam() or "",
                                 c.attrs.get("OBJNAM") or "", c.klass or "")
                    report.add(f)

    report.finished = datetime.now(timezone.utc)
    return report


def reconcile_lines(chart_features, source_features, recipe, catalogue,
                    coverage=None):
    """Line-geometry reconciliation (quay walls, shoreline construction, …).

    Pairs charted lines with source lines (``geomatch.match_lines``) and judges
    each pair's geometry with the generalization DEAD-BAND (``geomatch.
    compare_lines`` → SRC.5): within-tolerance and almost-similar shapes pass
    silently; only an obvious difference — a contiguous stretch outside the
    tolerance corridor, or a worst deviation beyond the cap — is flagged. The
    presence checks (SRC.1/SRC.2) and the attribute comparison on matched pairs
    (SRC.3) run exactly as for points.

    For a line recipe ``recipe.buffer_m`` is the geometry tolerance in metres
    (the corridor half-width is 2× it — see geomatch's constants). Features on
    both sides must carry ``coords`` = [(lon, lat), …]; their lon/lat is the
    representative midpoint used for coverage clipping and report locations.
    """
    from . import geomatch
    report = Report()
    report.ran = list(checks.LINE_CHECKS)
    tol = recipe.buffer_m
    klass = recipe.enc_class

    # Known-fault exclusions — same discipline (and disclosure) as points.
    ex_classes = recipe.excluded_class_set()
    ex_lnam = recipe.excluded_lnam_set()
    ex_source = recipe.excluded_source_set()
    report.excluded_classes = sorted(ex_classes)
    report.enc_excluded_classes_count = 0
    report.enc_excluded_features_count = 0
    report.source_excluded_count = 0
    report.excluded_enc_list = [dict(e) for e in recipe.enc_exclude]
    report.excluded_source_list = [dict(e) for e in recipe.source_exclude]
    if ex_classes or ex_lnam:
        kept = []
        for c in chart_features:
            if (c.klass or "").upper() in ex_classes:
                report.enc_excluded_classes_count += 1
                continue
            ln = c.lnam()
            if ln and ln.upper() in ex_lnam:
                report.enc_excluded_features_count += 1
                continue
            kept.append(c)
        chart_features = kept
    if ex_source:
        kept = []
        for s in source_features:
            if str(getattr(s, "ident", "")) in ex_source:
                report.source_excluded_count += 1
                continue
            kept.append(s)
        source_features = kept

    # Coverage clip for the "in source, not in chart" side (midpoint test).
    use_coverage = bool(coverage)
    area = None
    if not use_coverage and getattr(recipe, "restrict_to_area", True):
        area = _chart_area(chart_features, getattr(recipe, "area_margin_m",
                                                   250.0))
    report.clip_mode = ("coverage" if use_coverage
                        else "bbox" if area is not None else "none")

    def in_area(s):
        if use_coverage:
            return point_in_coverage(s.lon, s.lat, coverage)
        if area is None:
            return True
        return area[0] <= s.lon <= area[2] and area[1] <= s.lat <= area[3]

    chart = [c for c in chart_features
             if getattr(c, "coords", None) and len(c.coords) >= 2]
    src = [s for s in source_features
           if getattr(s, "coords", None) and len(s.coords) >= 2]

    pairs, un_chart, un_src = geomatch.match_lines(
        [c.coords for c in chart], [s.coords for s in src], tol)

    def _enc_excl(c):
        return ("enc", c.lnam() or "", c.attrs.get("OBJNAM") or "",
                c.klass or "")

    # --- matched pairs: SRC.5 geometry dead-band + SRC.3 attributes ---------
    for i, j in pairs:
        c, s = chart[i], src[j]
        r = geomatch.compare_lines(c.coords, s.coords, tol)
        if r["flag"]:
            wr = r["worst_run"]
            det = (f"{int(round(100 * r['fraction']))}% of "
                   f"{r['length_m']:.0f} m within "
                   f"{geomatch.CORRIDOR_X * tol:.1f} m of the source; "
                   f"worst {r['hausdorff_m']:.1f} m")
            if wr:
                det += f" (a {wr[0]:.0f} m stretch leaves the corridor)"
            pt = wr[2] if wr else (c.lon, c.lat)
            f = Finding(checks.GEOMETRY_DEVIATION, c.cell, klass, c.fid,
                        detail=det,
                        geom_wkt=f"POINT ({pt[0]:.7f} {pt[1]:.7f})")
            f.exclude = _enc_excl(c)
            report.add(f)
        for _attr, detail in _compare_attrs(c, s, recipe, catalogue):
            f = Finding(checks.VALUE_MISMATCH, c.cell, klass, c.fid,
                        detail=detail, geom_wkt=c.wkt())
            f.exclude = _enc_excl(c)
            report.add(f)

    # --- unmatched chart lines: nothing plausible in the source -------------
    for i in un_chart:
        c = chart[i]
        f = Finding(checks.IN_CHART_NOT_SOURCE, c.cell, klass, c.fid,
                    detail=f"no source line within "
                           f"{geomatch.MATCH_GATE_X * tol:.0f} m "
                           f"@ {c.lat:.5f}, {c.lon:.5f}",
                    geom_wkt=c.wkt())
        f.exclude = _enc_excl(c)
        report.add(f)

    # --- unmatched source lines inside the chart area -----------------------
    out_of_area = 0
    for j in un_src:
        s = src[j]
        if not in_area(s):
            out_of_area += 1
            continue
        f = Finding(checks.IN_SOURCE_NOT_CHART, "(source)", klass, s.ident,
                    detail=f"source id {s.ident} @ {s.lat:.5f}, {s.lon:.5f}",
                    geom_wkt=f"POINT ({s.lon:.7f} {s.lat:.7f})")
        f.exclude = ("source", str(s.ident), "")
        report.add(f)
    report.source_out_of_area = out_of_area

    report.finished = datetime.now(timezone.utc)
    return report


def reconcile_areas(chart_features, source_features, recipe, catalogue,
                    coverage=None):
    """Area-geometry reconciliation (basins, land areas vs source polygons).

    The source polygons are FIRST CLIPPED to the chart's coverage (M_COVR) —
    both sides then describe the same footprint, so a source polygon running
    past the cell border never reads as a disagreement (see areamatch's module
    docstring for the full design). The verdict is the coverage-clipped
    symmetric difference + IoU inside a dead-band; only OBVIOUS disagreement
    is flagged (SRC.5), with the surface disagreement reported in m²/ha and a
    zoom-to point on the largest lobe. Presence (SRC.1/SRC.2) and attribute
    checks (SRC.3) run on matched pairs as for points and lines.

    ``recipe.buffer_m`` is the geometry tolerance in metres (drives the
    sliver-erosion width and the lobe floor). ``coverage`` is the prepared
    coverage set from :func:`prepare_coverage`. No-op findings when OGR/GEOS
    is unavailable (mirrors the geometry checks)."""
    from . import areamatch
    report = Report()
    report.ran = list(checks.LINE_CHECKS)
    tol = recipe.buffer_m
    klass = recipe.enc_class

    # Known-fault exclusions — same discipline (and disclosure) as points.
    ex_classes = recipe.excluded_class_set()
    ex_lnam = recipe.excluded_lnam_set()
    ex_source = recipe.excluded_source_set()
    report.excluded_classes = sorted(ex_classes)
    report.enc_excluded_classes_count = 0
    report.enc_excluded_features_count = 0
    report.source_excluded_count = 0
    report.excluded_enc_list = [dict(e) for e in recipe.enc_exclude]
    report.excluded_source_list = [dict(e) for e in recipe.source_exclude]
    if ex_classes or ex_lnam:
        kept = []
        for c in chart_features:
            if (c.klass or "").upper() in ex_classes:
                report.enc_excluded_classes_count += 1
                continue
            ln = c.lnam()
            if ln and ln.upper() in ex_lnam:
                report.enc_excluded_features_count += 1
                continue
            kept.append(c)
        chart_features = kept
    if ex_source:
        kept = []
        for s in source_features:
            if str(getattr(s, "ident", "")) in ex_source:
                report.source_excluded_count += 1
                continue
            kept.append(s)
        source_features = kept

    use_coverage = bool(coverage)
    # areamatch wants raw ring sets; prepared coverage is (ext, holes, bbox).
    cov_rings = ([[c[0]] + list(c[1] or []) for c in coverage]
                 if use_coverage else None)
    report.clip_mode = "coverage" if use_coverage else "none"

    def in_area(s):
        if use_coverage:
            return point_in_coverage(s.lon, s.lat, coverage)
        return True

    chart = [c for c in chart_features if getattr(c, "areas", None)]
    src = [s for s in source_features if getattr(s, "areas", None)]

    matched = areamatch.match_areas([c.areas for c in chart],
                                    [s.areas for s in src], tol)
    if matched is None:                       # no OGR/GEOS: judge nothing
        report.finished = datetime.now(timezone.utc)
        return report
    pairs, un_chart, un_src = matched

    def _enc_excl(c):
        return ("enc", c.lnam() or "", c.attrs.get("OBJNAM") or "",
                c.klass or "")

    # --- matched pairs: SRC.5 dead-band + SRC.3 attributes ------------------
    for i, j in pairs:
        c, s = chart[i], src[j]
        r = areamatch.compare_areas(c.areas, s.areas, tol,
                                    coverage_ll=cov_rings)
        if r and r["flag"]:
            ha = r["symdiff_m2"] / 10000.0
            det = (f"IoU {r['iou']:.2f}; surface disagreement "
                   f"{ha:.2f} ha ({r['symdiff_m2']:.0f} m²)")
            lobe = r["worst_lobe"]
            if lobe:
                det += f"; largest lobe {lobe[0] / 10000.0:.2f} ha"
                pt = lobe[1]
            else:
                pt = (c.lon, c.lat)
            f = Finding(checks.GEOMETRY_DEVIATION, c.cell, klass, c.fid,
                        detail=det,
                        geom_wkt=f"POINT ({pt[0]:.7f} {pt[1]:.7f})")
            f.exclude = _enc_excl(c)
            report.add(f)
        for _attr, detail in _compare_attrs(c, s, recipe, catalogue):
            f = Finding(checks.VALUE_MISMATCH, c.cell, klass, c.fid,
                        detail=detail, geom_wkt=c.wkt())
            f.exclude = _enc_excl(c)
            report.add(f)

    # --- unmatched chart polygons --------------------------------------------
    for i in un_chart:
        c = chart[i]
        f = Finding(checks.IN_CHART_NOT_SOURCE, c.cell, klass, c.fid,
                    detail=f"no source polygon overlaps "
                           f"@ {c.lat:.5f}, {c.lon:.5f}",
                    geom_wkt=c.wkt())
        f.exclude = _enc_excl(c)
        report.add(f)

    # --- unmatched source polygons inside the chart area ---------------------
    out_of_area = 0
    for j in un_src:
        s = src[j]
        if not in_area(s):
            out_of_area += 1
            continue
        f = Finding(checks.IN_SOURCE_NOT_CHART, "(source)", klass, s.ident,
                    detail=f"source id {s.ident} @ {s.lat:.5f}, {s.lon:.5f}",
                    geom_wkt=f"POINT ({s.lon:.7f} {s.lat:.7f})")
        f.exclude = ("source", str(s.ident), "")
        report.add(f)
    report.source_out_of_area = out_of_area

    report.finished = datetime.now(timezone.utc)
    return report
