# 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
