"""Fairway network engine — bind live AIS positions to the Dutch fairway network.

Pure Python (stdlib + networkx only): no shapely / GDAL / pyproj, so it runs and
is testable anywhere — the QGIS plugin, an offline test, or a future field node /
central server. QGIS is used only to *render* results, never to compute them.

Pipeline:
  load_network(gpkg)          read 'vaarwegen' (lines) + 'junctions' (points)
  FairwayNetwork.build(...)   project-and-split noding -> routable graph (networkx)
                              + a segment grid index carrying linear-referencing
  net.match(lon, lat, cog)    -> Match(routeid, route_name, km, offset_m, ...) | None

Geometry is done in a local equirectangular metric frame (deg -> m via cos(lat0)),
accurate to ~0.1-0.2 % over the Netherlands — ample for snapping/offset work where
distances are metres. Junctions lie ON the lines (not on vertices), so noding splits
each line at the junction's *projected* position.
"""
import math
import re
import sqlite3
import struct

try:
    import networkx as nx
except Exception:                       # networkx optional for pure matching
    nx = None


# Major NL fairway entrances from the North Sea / Wadden / estuaries. Used to
# orient "toward the sea". Editable — adding/refining gateways tweaks the few
# routes that come out wrong. (name, lon, lat)
DEFAULT_SEA_GATEWAYS = [
    ("Westerschelde mouth", 3.36, 51.45),
    ("Maasmond (Hoek van Holland)", 4.02, 51.99),
    ("IJmuiden", 4.55, 52.46),
    ("Den Helder (Marsdiep)", 4.78, 52.96),
    ("Vliestroom (Harlingen)", 5.20, 53.30),
    ("Lauwersoog", 6.20, 53.42),
    ("Eems (Delfzijl)", 6.83, 53.46),
    ("Den Oever (Afsluitdijk W)", 5.04, 52.93),
    ("Kornwerderzand (Afsluitdijk E)", 5.34, 53.07),
]


# Optional route-name overrides, keyed by routeid. The binding network is built
# from the RWS ``route`` layer, where each routeid is already ONE full-stretch
# route with a proper name (e.g. 20940 = 'Westerschelde'), so this is normally
# empty — add entries only to relabel a specific route in the route manager.
ROUTE_NAME_OVERRIDES = {}

# Custom (user-drawn) routes carry NEGATIVE routeids and are ranked ahead of every
# native route in candidates() by subtracting this large bonus from their cost —
# so when a vessel is within tolerance of both a custom and a native route, it
# binds to the custom one. Native binding is unchanged when no custom is in range.
_CUSTOM_PRIORITY = 1.0e7


# --------------------------------------------------------------------------- #
#  GeoPackage reading (sqlite3 + manual WKB decode)
# --------------------------------------------------------------------------- #
def _parse_gpkg_geom(blob):
    """Decode a GeoPackage geometry blob -> ('Point'|'LineString'|..., coords)."""
    if blob[:2] != b"GP":
        raise ValueError("not a GeoPackage geometry blob")
    flags = blob[3]
    env = (flags >> 1) & 0x07
    envsz = {0: 0, 1: 32, 2: 48, 3: 48, 4: 64}[env]
    off = 8 + envsz
    bo = "<" if blob[off] == 1 else ">"
    off += 1
    (typ,) = struct.unpack(bo + "I", blob[off:off + 4]); off += 4
    base = typ % 1000                   # strip ISO Z/M thousands
    if base == 1:                       # Point
        x, y = struct.unpack(bo + "dd", blob[off:off + 16])
        return "Point", (x, y)
    if base == 2:                       # LineString
        (n,) = struct.unpack(bo + "I", blob[off:off + 4]); off += 4
        v = struct.unpack(bo + "d" * (2 * n), blob[off:off + 16 * n])
        return "LineString", list(zip(v[0::2], v[1::2]))
    if base == 5:                       # MultiLineString
        (ng,) = struct.unpack(bo + "I", blob[off:off + 4]); off += 4
        parts = []
        for _ in range(ng):
            bo2 = "<" if blob[off] == 1 else ">"; off += 1
            off += 4                    # inner type
            (n,) = struct.unpack(bo2 + "I", blob[off:off + 4]); off += 4
            v = struct.unpack(bo2 + "d" * (2 * n), blob[off:off + 16 * n])
            off += 16 * n
            parts.append(list(zip(v[0::2], v[1::2])))
        return "MultiLineString", parts
    raise ValueError("unsupported geometry type %d" % base)


def _safe_ident(name):
    """Validate a GeoPackage table name as a bare SQL identifier before it is
    interpolated into a query (table names cannot be bound as parameters). Our
    callers only ever pass fixed constants ('route'/'junctions'/'vaarwegen'); this
    guard makes that structurally injection-proof."""
    if not re.match(r"^[A-Za-z_][A-Za-z0-9_]*$", name or ""):
        raise ValueError("unsafe table name: %r" % (name,))
    return name


def load_network(gpkg_path, lines_table="vaarwegen", junc_table="junctions"):
    """Read fairway lines and junction points from a GeoPackage.

    Returns (lines, junctions):
      lines     : list of dicts {fid, routeid, name, kmbegin, kmend, coords[(lon,lat)]}
      junctions : list of (fid, (lon, lat))
    """
    lines_table = _safe_ident(lines_table)
    junc_table = _safe_ident(junc_table)
    con = sqlite3.connect(gpkg_path)
    cur = con.cursor()
    lines = []
    for fid, rid, name, kmb, kme, blob in cur.execute(
            f'SELECT fid, routeid, name, routekmbegin, routekmend, geom '
            f'FROM "{lines_table}"'):        # nosec B608 - table name is a validated identifier
        kind, coords = _parse_gpkg_geom(blob)
        if kind == "LineString":
            segs = [coords]
        elif kind == "MultiLineString":
            segs = coords
        else:
            continue
        for part in segs:
            if len(part) >= 2:
                lines.append({"fid": fid, "routeid": rid, "name": name,
                              "kmbegin": kmb, "kmend": kme, "coords": part})
    junctions = []
    for fid, blob in cur.execute(
            f'SELECT fid, geom FROM "{junc_table}"'):   # nosec B608 - validated identifier
        kind, coords = _parse_gpkg_geom(blob)
        if kind == "Point":
            junctions.append((fid, coords))
    con.close()
    return lines, junctions


# --------------------------------------------------------------------------- #
#  Geometry helpers (local equirectangular metric frame)
# --------------------------------------------------------------------------- #
def _project_point_on_seg(px, py, ax, ay, bx, by):
    """Return (distance, t, cx, cy): nearest point on segment a->b to p."""
    dx, dy = bx - ax, by - ay
    L2 = dx * dx + dy * dy
    if L2 == 0.0:
        return math.hypot(px - ax, py - ay), 0.0, ax, ay
    t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / L2))
    cx, cy = ax + t * dx, ay + t * dy
    return math.hypot(px - cx, py - cy), t, cx, cy


def _bearing_deg(ax, ay, bx, by):
    """Compass bearing (0=N, clockwise) of a->b in the metric frame."""
    return math.degrees(math.atan2(bx - ax, by - ay)) % 360.0


def _ang_diff_180(a, b):
    """Smallest angle between two bearings, folded to 0..90 (direction-agnostic)."""
    d = abs((a - b) % 180.0)
    return min(d, 180.0 - d)


class Match:
    """Result of binding a position to the fairway network."""
    __slots__ = ("routeid", "route_name", "km", "offset_m", "snap_lon",
                 "snap_lat", "bearing", "line_fid", "axis_side")

    def __init__(self, routeid, route_name, km, offset_m, snap_lon, snap_lat,
                 bearing, line_fid, axis_side=0):
        """Hold one position-to-fairway binding: which route, the linear-
        referenced km, the perpendicular offset, the snapped point, the axis
        bearing in the +km direction, and which side of the axis the fix lies."""
        self.routeid = routeid
        self.route_name = route_name
        self.km = km
        self.offset_m = offset_m
        self.snap_lon = snap_lon
        self.snap_lat = snap_lat
        self.bearing = bearing          # bearing of the segment in +km direction
        self.line_fid = line_fid
        self.axis_side = axis_side       # +1 left of +km dir, -1 right, 0 on-axis

    def __repr__(self):
        """Compact debug view: route id, name, km and offset."""
        return ("Match(route=%s '%s' km=%.3f off=%.1fm)"
                % (self.routeid, self.route_name, self.km, self.offset_m))


# --------------------------------------------------------------------------- #
#  The network
# --------------------------------------------------------------------------- #
class FairwayNetwork:
    def __init__(self, lat0, cos0, seg_index, cell, graph=None, meta=None):
        """Hold a built network: the metric-frame constants (lat0/cos0), the
        grid-indexed segment records with linear-referencing, the optional
        networkx topology graph, and the meta dict (node coords, route names,
        per-route end nodes, sea-orientation flags)."""
        self.lat0 = lat0
        self.cos0 = cos0
        self.seg_index = seg_index      # cell -> list of segment records
        self.cell = cell                # grid cell size (metres)
        self.graph = graph              # networkx MultiGraph (topology)
        self.meta = meta or {}
        # user-drawn custom routes (negative routeids), indexed like seg_index;
        # populated at runtime by set_custom_routes(). They take binding priority.
        self.custom_index = {}          # cell -> list of segment records

    # ---- metric frame ---------------------------------------------------- #
    _MDEG = 111320.0

    def to_m(self, lon, lat):
        """Project lon/lat (deg) to local equirectangular metres (x east, y north)."""
        return lon * self._MDEG * self.cos0, lat * self._MDEG

    def to_lonlat(self, x, y):
        """Inverse of to_m: local metres back to lon/lat degrees."""
        return x / (self._MDEG * self.cos0), y / self._MDEG

    def _cell(self, x, y):
        """Grid-cell key (ix, iy) for a metric point, for the segment index."""
        return (int(x // self.cell), int(y // self.cell))

    # ---- build ----------------------------------------------------------- #
    @classmethod
    def build(cls, lines, junctions, cell=300.0, snap_m=1.0, build_graph=True):
        """Build the network from loaded lines + junctions.

        Picks a mean-latitude metric frame, indexes every line segment into a
        grid (carrying linear-referencing metadata), and — when networkx is
        present — nodes the lines by project-and-split at each on-line junction
        to produce a routable graph, then labels routes by sea orientation.
        """
        lat0 = sum(ln["coords"][0][1] for ln in lines) / max(len(lines), 1)
        cos0 = math.cos(math.radians(lat0))
        mdeg = cls._MDEG

        def to_m(lon, lat):
            """Local metric projection using this build's frame (closure)."""
            return lon * mdeg * cos0, lat * mdeg

        # junction spatial grid (metres) for project-and-split
        jgrid = {}
        for jfid, (jlon, jlat) in junctions:
            jx, jy = to_m(jlon, jlat)
            jgrid.setdefault((int(jx // cell), int(jy // cell)), []).append(
                (jfid, jx, jy))

        seg_index = {}
        graph = nx.MultiGraph() if (build_graph and nx is not None) else None
        jcoord = {jfid: (lo, la) for jfid, (lo, la) in junctions}
        node_xy = {}                    # graph node -> (lon, lat)
        route_ends = {}                 # routeid -> [(km, node), ...]

        for ln in lines:
            coords = ln["coords"]
            m = [to_m(lon, lat) for lon, lat in coords]
            cum = [0.0]
            for i in range(len(m) - 1):
                cum.append(cum[-1] + math.hypot(m[i + 1][0] - m[i][0],
                                                m[i + 1][1] - m[i][1]))
            total = cum[-1] or 1.0
            kmb, kme = ln["kmbegin"], ln["kmend"]
            if kmb is None:
                kmb = 0.0
            if kme is None:
                kme = total / 1000.0

            # index every segment with linear-referencing metadata
            for i in range(len(m) - 1):
                ax, ay = m[i]; bx, by = m[i + 1]
                rec = (ax, ay, bx, by, ln["fid"], ln["routeid"], ln["name"],
                       kmb, kme, total, cum[i],
                       coords[i][0], coords[i][1], coords[i + 1][0], coords[i + 1][1])
                c0 = (int(ax // cell), int(ay // cell))
                c1 = (int(bx // cell), int(by // cell))
                for cx in range(min(c0[0], c1[0]), max(c0[0], c1[0]) + 1):
                    for cy in range(min(c0[1], c1[1]), max(c0[1], c1[1]) + 1):
                        seg_index.setdefault((cx, cy), []).append(rec)

            # project-and-split: find junctions lying on this line -> graph edges
            if graph is not None:
                onsite, seen = [], set()
                for i in range(len(m) - 1):
                    ax, ay = m[i]; bx, by = m[i + 1]
                    c0 = (int(ax // cell), int(ay // cell))
                    c1 = (int(bx // cell), int(by // cell))
                    for cx in range(min(c0[0], c1[0]) - 1, max(c0[0], c1[0]) + 2):
                        for cy in range(min(c0[1], c1[1]) - 1, max(c0[1], c1[1]) + 2):
                            for (jfid, jx, jy) in jgrid.get((cx, cy), ()):
                                if jfid in seen:
                                    continue
                                d, t, _, _ = _project_point_on_seg(jx, jy, ax, ay, bx, by)
                                if d <= snap_m:
                                    along = cum[i] + t * math.hypot(bx - ax, by - ay)
                                    onsite.append((along, ("j", jfid)))
                                    seen.add(jfid)

                def endnode(idx):
                    """Node id for a line end: reuse a junction sitting on the
                    endpoint, else a coordinate-keyed endpoint node (closure)."""
                    target = 0.0 if idx == 0 else total
                    for along, nid in onsite:
                        if abs(along - target) <= snap_m:
                            return nid
                    p = coords[idx]
                    return ("e", round(p[0], 5), round(p[1], 5))

                seq = [(0.0, endnode(0))] + sorted(onsite) + [(total, endnode(-1))]
                dedup = []
                for along, nid in seq:
                    if not dedup or dedup[-1][1] != nid:
                        dedup.append((along, nid))
                for (a0, n0), (a1, n1) in zip(dedup, dedup[1:]):
                    if n0 != n1:
                        graph.add_edge(n0, n1, routeid=ln["routeid"],
                                       length_m=a1 - a0, fid=ln["fid"])
                # node coords + per-route end nodes, for sea orientation later
                for _a, nid in dedup:
                    if nid not in node_xy:
                        node_xy[nid] = (jcoord.get(nid[1]) if nid[0] == "j"
                                        else (nid[1], nid[2]))
                if dedup:
                    re = route_ends.setdefault(ln["routeid"], [])
                    re.append((kmb, dedup[0][1]))
                    re.append((kme, dedup[-1][1]))

        # route-manager name per routeid: a deliberate override where set, else
        # the (last-seen) segment name. Per-segment names stay on the segments.
        route_names = {ln["routeid"]: ln["name"] for ln in lines}
        for rid, nm in ROUTE_NAME_OVERRIDES.items():
            if rid in route_names:
                route_names[rid] = nm
        meta = {"n_lines": len(lines), "n_junctions": len(junctions),
                "node_xy": node_xy, "route_ends": route_ends, "route_seaward": {},
                "route_names": route_names}
        net = cls(lat0, cos0, seg_index, cell, graph, meta)
        if graph is not None:
            net.compute_sea_orientation()
        return net

    # ---- sea orientation (inbound / outbound) ---------------------------- #
    def compute_sea_orientation(self, gateways=None, snap_max_m=15000.0,
                                min_diff_m=500.0):
        """Label each route with whether *increasing km heads toward the sea*.

        Snaps the sea gateways to network nodes, computes each node's network
        distance to the nearest sea, then per route compares that distance at its
        low-km and high-km ends. Stored as meta['route_seaward'][routeid] =
        True (km rises seaward) / False (km rises inland) / None (ambiguous).
        """
        if self.graph is None or nx is None:
            return
        gateways = gateways or DEFAULT_SEA_GATEWAYS
        node_xy = self.meta.get("node_xy", {})
        route_ends = self.meta.get("route_ends", {})
        items = [(nid, xy) for nid, xy in node_xy.items() if xy]
        sea = set()
        for _nm, lo, la in gateways:
            px, py = self.to_m(lo, la)
            best, bd = None, float("inf")
            for nid, xy in items:
                qx, qy = self.to_m(*xy)
                dd = math.hypot(px - qx, py - qy)
                if dd < bd:
                    bd, best = dd, nid
            if best is not None and bd <= snap_max_m:
                sea.add(best)
        if not sea:
            return
        dist = nx.multi_source_dijkstra_path_length(self.graph, sea,
                                                    weight="length_m")
        flags = {}
        for rid, ends in route_ends.items():
            lo = min(ends, key=lambda e: e[0])[1]
            hi = max(ends, key=lambda e: e[0])[1]
            dl, dh = dist.get(lo), dist.get(hi)
            if dl is None or dh is None or abs(dh - dl) < min_diff_m:
                flags[rid] = None
            else:
                flags[rid] = dh < dl
        self.meta["route_seaward"] = flags
        self.meta["sea_nodes"] = list(sea)

    # ---- match ----------------------------------------------------------- #
    def candidates(self, lon, lat, cog=None, max_offset_m=150.0,
                   cog_weight=2.0, k=6, route_tol=None):
        """Ranked nearby routes (best per route). Returns list of (cost, Match).

        cost = offset_m + cog_weight * bearing-misalignment(deg) when cog given.
        Deduplicated to one entry per routeid (the closest segment of each), so
        the stateful matcher reasons about distinct routes, not raw segments.

        ``route_tol`` is an optional {routeid: max_offset_m} override — a route
        gets its own perpendicular tolerance (so a wide tidal fairway can bind
        vessels sailing well off the centreline) while every other route keeps
        ``max_offset_m``. The cell search widens to cover the largest tolerance
        in play, so the default (no overrides) is unchanged.

        Custom (user-drawn) routes in ``custom_index`` are scanned alongside the
        native network and, having negative routeids, are ranked ahead of any
        native route within tolerance (see ``_CUSTOM_PRIORITY``).
        """
        px, py = self.to_m(lon, lat)
        cx, cy = self._cell(px, py)
        search = max_offset_m
        if route_tol:
            search = max(search, max(route_tol.values()))
        rad = max(1, int(math.ceil(search / self.cell)))
        best = {}                       # routeid -> (cost, Match)

        def consider(rec):
            """Score one segment record and keep it if it's the best for its
            route (custom routes get the priority bonus)."""
            ax, ay, bx, by = rec[0], rec[1], rec[2], rec[3]
            d, t, qx, qy = _project_point_on_seg(px, py, ax, ay, bx, by)
            rid = rec[5]
            allowed = route_tol.get(rid, max_offset_m) if route_tol \
                else max_offset_m
            if d > allowed:
                return
            seg_bear = _bearing_deg(ax, ay, bx, by)
            cost = d
            if cog is not None:
                cost += cog_weight * _ang_diff_180(seg_bear, cog)
            if rid < 0:                 # custom route: outrank every native route
                cost -= _CUSTOM_PRIORITY
            prev = best.get(rid)
            if prev is not None and prev[0] <= cost:
                return
            (_, _, _, _, fid, _, name, kmb, kme, total, cum_i,
             alon, alat, blon, blat) = rec
            seg_len = math.hypot(bx - ax, by - ay)
            along = cum_i + t * seg_len
            km = kmb + (along / total) * (kme - kmb)
            snap_lon = alon + t * (blon - alon)
            snap_lat = alat + t * (blat - alat)
            cross = (bx - ax) * (py - ay) - (by - ay) * (px - ax)
            side = 1 if cross > 0 else (-1 if cross < 0 else 0)
            best[rid] = (cost, Match(rid, name, km, d, snap_lon,
                                     snap_lat, seg_bear, fid, side))

        for ddx in range(-rad, rad + 1):
            for ddy in range(-rad, rad + 1):
                cell = (cx + ddx, cy + ddy)
                for rec in self.seg_index.get(cell, ()):
                    consider(rec)
                for rec in self.custom_index.get(cell, ()):
                    consider(rec)
        return sorted(best.values(), key=lambda cm: cm[0])[:k]

    def match(self, lon, lat, cog=None, max_offset_m=150.0, cog_weight=2.0):
        """Bind a position to the single nearest fairway. Returns Match or None."""
        c = self.candidates(lon, lat, cog, max_offset_m, cog_weight, k=1)
        return c[0][1] if c else None


def build_or_load(gpkg_path, cache_path=None, lines_table="vaarwegen",
                  junc_table="junctions", **kw):
    """Build the network from the GeoPackage.

    ``lines_table`` / ``junc_table`` select the GeoPackage tables (the RWS
    ``route`` resource uses 'route' + 'junctions'); remaining kwargs go to
    ``FairwayNetwork.build`` (e.g. ``snap_m``). The build is fast (~1-2 s) and is
    done fresh each session — there is no on-disk cache (a pickled cache would be
    an untrusted-deserialization sink). ``cache_path`` is accepted and ignored
    for backward compatibility."""
    lines, junctions = load_network(gpkg_path, lines_table, junc_table)
    return FairwayNetwork.build(lines, junctions, **kw)
