"""Routing graph noding (docs/offroad/PLAN.md R10.1a, Appendix E.7; platform since docs/backcountry/PLAN.md B13.2). PC
only, never on the VM. Mode-independent: each app's graph build (offroad-graph/official_graph.py, osm_graph.py) reads
its rows, calls read_parts and build here, then adds its own per-edge attributes.

read_parts takes rows as GeoJSON features with source, source_id, id (and number, name) properties; with require_clean
any tier other than clean fails the build (RJ-18). It explodes them into parts (parts under 1 m are dropped) and build
nodes them, byte-identically on every run, into:

  edges        one LineString per edge: a, b (node ids, nodeid.py), src, sid, rid (the row id), len (m), sorted by
               (a, b, src, sid);
  gaps         each dead end 3-50 m from another route, never bridged (for track RT);
  crossings    each crossing without a shared vertex: noded between sources, not within one (QA, R10.8), with runs of
               more than 3 per 100 m between two parts flagged as duplicate drawing;
  counts       for the build's report.json.

credits() writes credits.json: every source in the graph with its credit, plus the terrain's.

Noding rules (E.7): (a) a vertex shared by parts is a node where the network branches (other than two distinct
neighbouring vertices, or a part ends there), so parts drawn vertex for vertex on the same road bed keep parallel edges;
(b) every part endpoint is a node; (c) an endpoint within GRAPH_SNAP_V1 (3 m) of another part's interior splits it there
(onto that part's node instead when one is within 3 m); (d) endpoints within 3 m merge (the first in sort order keeps its
place); (e) larger gaps are never bridged; (f) bare crossings are noded only between different sources.

Distances are metres in a local equirectangular frame around each tested point (k = cos of its latitude), and candidates
come from a k-d tree over 3D points (chord = ground distance at these scales) of 5 m pieces of every segment; so there
are no chunk seams anywhere. A crossing on a piece boundary is found from both pieces and kept once. Python stdlib +
NumPy/SciPy only (no pgRouting, no GPL code).
"""
import json
import math
import os
import sys

import numpy as np
from scipy.spatial import cKDTree

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path[:0] = [HERE, os.path.join(HERE, '..', '..', 'lib')]  # route-graph/lib, then tools/geodata/lib (terrarium)
from nodeid import node_id  # noqa: E402

R = 6371008.8  # WGS84 mean radius (m): haversine lengths and the 3D index
DEG = math.pi / 180 * R  # metres per degree of latitude
GRAPH_SNAP_V1 = 3.0
GAP_MAX_M = 50.0
SAMPLE_M = 5.0
MIN_PART_M = 1.0
VERTEX_LAND_M = 0.5  # a T-snap this close to a vertex lands on the vertex
SAME_POINT_M = 0.05  # crossings this close to a node, or to each other, are that node / one crossing
DUP_COUNT, DUP_WINDOW_M = 3, 100.0


class BuildError(Exception):
    pass


def haversine(a, b):
    lon1, lat1, lon2, lat2 = (np.radians(v) for v in (a[..., 0], a[..., 1], b[..., 0], b[..., 1]))
    h = np.sin((lat2 - lat1) / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2
    return 2 * R * np.arcsin(np.sqrt(np.minimum(h, 1.0)))


def local(origin, pts):
    """Metres east and north of `origin` (lon, lat; broadcasts) for `pts`."""
    k = np.cos(np.radians(origin[..., 1]))
    return np.stack(((pts[..., 0] - origin[..., 0]) * k * DEG, (pts[..., 1] - origin[..., 1]) * DEG), -1)


def unit3(ll):
    lon, lat = np.radians(ll[..., 0]), np.radians(ll[..., 1])
    return np.stack((np.cos(lat) * np.cos(lon), np.cos(lat) * np.sin(lon), np.sin(lat)), -1) * R


def cross2(u, v):
    return u[..., 0] * v[..., 1] - u[..., 1] * v[..., 0]


def row_key(p):
    rid = p.get('id')
    return (p.get('source') or '', p.get('source_id') or '', (0, rid, '') if isinstance(rid, int) else (1, 0, str(rid)))


def read_parts(lines, require_clean=True):
    """Clean rows in (source, source_id, id) order, exploded into parts. Returns (parts, meta, counts). The OSM graph
    (osm_graph.py, R10.3) has no tiers: OSM is never an official rule, so it passes require_clean=False."""
    rows = []
    for n, line in enumerate(lines, 1):
        if not line.strip():
            continue
        f = json.loads(line)
        p = f.get('properties') or {}
        if require_clean and p.get('tier') != 'clean':
            raise BuildError(f"row {p.get('id')} ({p.get('source')}) on input line {n} has tier {p.get('tier')!r}: "
                             'only clean rows enter the routing graph (RJ-18)')
        g = f.get('geometry')
        rows.append((row_key(p) + (json.dumps(g, sort_keys=True),), p, g))
    rows.sort(key=lambda r: r[0])
    parts, meta, short = [], [], 0
    for _, p, g in rows:
        if not g or g.get('type') not in ('LineString', 'MultiLineString'):
            continue
        for coords in [g['coordinates']] if g['type'] == 'LineString' else g['coordinates']:
            ll = np.round(np.array([c[:2] for c in coords], dtype=np.float64), 7) if coords else np.zeros((0, 2))
            if len(ll):
                ll = ll[np.concatenate([[True], np.any(np.diff(ll, axis=0) != 0, axis=1)])]
            if len(ll) < 2 or haversine(ll[:-1], ll[1:]).sum() < MIN_PART_M:
                short += 1
                continue
            parts.append(ll)
            meta.append({'src': p.get('source'), 'sid': p.get('source_id') or str(p.get('id')), 'rid': p.get('id'),
                         'num': p.get('number'), 'name': p.get('name'), 'row': p})
    return parts, meta, {'rows': len(rows), 'short_parts_dropped': short}


def key(xy):
    return (round(float(xy[0]), 7), round(float(xy[1]), 7))


class Pieces:
    """Every segment cut into pieces of at most SAMPLE_M, indexed by their midpoints."""

    def __init__(self, parts):
        part, seg, k, n, a, b = [], [], [], [], [], []
        for i, p in enumerate(parts):
            pa, pb = p[:-1], p[1:]
            cnt = np.maximum(1, np.ceil(haversine(pa, pb) / SAMPLE_M)).astype(np.int64)
            rep = np.repeat(np.arange(len(pa)), cnt)
            kk = np.concatenate([np.arange(c) for c in cnt])
            nn = np.repeat(cnt, cnt)
            a.append(pa[rep] + (pb[rep] - pa[rep]) * (kk / nn)[:, None])
            b.append(pa[rep] + (pb[rep] - pa[rep]) * ((kk + 1) / nn)[:, None])
            part.append(np.full(len(rep), i))
            seg.append(rep)
            k.append(kk)
            n.append(nn)
        self.part, self.seg, self.k, self.n, self.a, self.b = map(np.concatenate, (part, seg, k, n, a, b))
        self.tree = cKDTree(unit3((self.a + self.b) / 2))

    def position(self, c, t):
        """(segment, fraction along it) for piece c at its parameter t; the next vertex when it lands on it."""
        seg, frac = int(self.seg[c]), float((self.k[c] + t) / self.n[c])
        return (seg + 1, 0.0) if frac >= 1.0 else (seg, max(frac, 0.0))


def seg_point(p, a, b):
    """Distance (m) from p to each segment a-b and the parameter t in [0, 1], in a frame around p."""
    A, B = local(p, a), local(p, b)
    ab = B - A
    l2 = (ab ** 2).sum(-1)
    t = np.where(l2 > 0, -(A * ab).sum(-1) / np.where(l2 > 0, l2, 1), 0).clip(0, 1)
    return np.hypot(*(A + ab * t[..., None]).T), t


class Graph:
    def __init__(self, parts, meta, snap=GRAPH_SNAP_V1):
        self.parts, self.meta, self.snap = parts, meta, snap
        self.cuts = [dict() for _ in parts]  # part -> {(segment, fraction): (lon, lat)}
        self.moves = {}  # endpoint index (2 * part + 0|1) -> new (lon, lat)
        self.gaps, self.crossings = [], []
        self.qa = {'t_snaps': 0, 'endpoints_merged': 0, 'crossings_noded': 0, 'crossings_unnoded': 0}
        self.ends = np.array([p[k] for p in parts for k in (0, -1)])
        self.connected = set()  # endpoint indices joined to another part already

    def end_pos(self, e):
        return self.moves.get(e, self.ends[e])

    def part_nodes(self, i):
        return [self.end_pos(2 * i), self.end_pos(2 * i + 1), *self.cuts[i].values()]

    def node(self):
        self.branch_vertices()
        self.merge_endpoints()
        self.pieces = Pieces(self.parts)
        self.snap_t_junctions()
        self.node_crossings()

    def branch_vertices(self):
        """(a) and (b): shared vertices where the network branches, or where a part ends, are nodes."""
        owners, neighbours = {}, {}
        keys = [[key(v) for v in p] for p in self.parts]
        for i, ks in enumerate(keys):
            for j, k in enumerate(ks):
                owners.setdefault(k, set()).add(i)
                nb = neighbours.setdefault(k, set())
                if j:
                    nb.add(ks[j - 1])
                if j < len(ks) - 1:
                    nb.add(ks[j + 1])
        end_keys = {ks[0] for ks in keys} | {ks[-1] for ks in keys}
        for i, ks in enumerate(keys):
            for j in range(1, len(ks) - 1):
                if len(neighbours[ks[j]]) != 2 or ks[j] in end_keys:
                    self.cuts[i][(j, 0.0)] = tuple(self.parts[i][j])
            for e, k in ((2 * i, ks[0]), (2 * i + 1, ks[-1])):
                if len(owners[k]) > 1:
                    self.connected.add(e)

    def merge_endpoints(self):
        """(d): endpoints of different parts within the snap distance become one node, at the first one's place."""
        parent = list(range(len(self.ends)))

        def find(x):
            while parent[x] != x:
                parent[x] = parent[parent[x]]
                x = parent[x]
            return x

        for e1, e2 in cKDTree(unit3(self.ends)).query_pairs(self.snap, output_type='ndarray'):
            if e1 // 2 != e2 // 2 and key(self.ends[e1]) != key(self.ends[e2]):
                r1, r2 = find(int(e1)), find(int(e2))
                if r1 != r2:
                    parent[max(r1, r2)] = min(r1, r2)
        for e in range(len(self.ends)):
            r = find(e)
            if r != e:
                self.connected.update((e, r))
                if key(self.ends[r]) != key(self.ends[e]):
                    self.moves[e] = tuple(self.ends[r])
                    self.qa['endpoints_merged'] += 1

    def snap_t_junctions(self):
        """(c) and (e): a dead end within the snap distance of another part joins it there; 3-50 m is only reported."""
        pc = self.pieces
        for e in range(len(self.ends)):
            if e in self.connected:
                continue
            i, pe = e // 2, np.array(self.end_pos(e))
            cand = np.array([c for c in pc.tree.query_ball_point(unit3(pe), GAP_MAX_M + SAMPLE_M) if pc.part[c] != i],
                            dtype=np.int64)
            if not len(cand):
                continue
            d, t = seg_point(pe, pc.a[cand], pc.b[cand])
            # nearest first; ties go to (source, source_id), then sort order
            order = sorted(range(len(cand)), key=lambda x: (round(float(d[x]), 6), self.meta[pc.part[cand[x]]]['src'] or '',
                                                            self.meta[pc.part[cand[x]]]['sid'] or '', int(pc.part[cand[x]]), int(cand[x])))
            best = order[0]
            c, dist = int(cand[best]), float(d[best])
            other = int(pc.part[c])
            if dist > GAP_MAX_M:
                continue
            if dist > self.snap:
                self.gaps.append((e, dist, other))
                continue
            nodes = np.array(self.part_nodes(other))
            nd = np.hypot(*local(pe, nodes).T)
            if nd.min() <= self.snap:
                q = tuple(nodes[int(np.argmin(nd))])  # onto that part's node, never a new cut beside it
            else:
                p = self.parts[other]
                seg, frac = pc.position(c, float(t[best]))
                q = tuple(pc.a[c] + (pc.b[c] - pc.a[c]) * t[best])
                for j in (seg, seg + 1):
                    if j < len(p) and np.hypot(*local(np.array(q), p[j])) < VERTEX_LAND_M:
                        q, seg, frac = tuple(p[j]), j, 0.0
                        break
                if not (frac == 0.0 and seg in (0, len(p) - 1)):
                    self.cuts[other][(seg, frac)] = q
            self.moves[e] = q
            self.connected.add(e)
            self.qa['t_snaps'] += 1

    def node_crossings(self):
        """(f): bare crossings are noded between different sources and only counted within one source."""
        pc = self.pieces
        pairs = pc.tree.query_pairs(SAMPLE_M * (1 + 1e-9), output_type='ndarray')
        pairs = pairs[pc.part[pairs[:, 0]] != pc.part[pairs[:, 1]]]
        if not len(pairs):
            return
        c1, c2 = pairs[:, 0], pairs[:, 1]
        o = pc.a[c1]
        r, a2, b2 = local(o, pc.b[c1]), local(o, pc.a[c2]), local(o, pc.b[c2])
        s = b2 - a2
        den = cross2(r, s)
        with np.errstate(divide='ignore', invalid='ignore'):
            t = cross2(a2, s) / den
            u = cross2(a2, r) / den
        eps = 1e-9
        # parallel (overlapping runs drawn vertex for vertex) never cross
        ok = (np.abs(den) > 1e-9 * np.hypot(*r.T) * np.hypot(*s.T)) & (t >= -eps) & (t <= 1 + eps) & (u >= -eps) & (u <= 1 + eps)
        found = {}
        for k in np.nonzero(ok)[0]:
            x = pc.a[c1[k]] + (pc.b[c1[k]] - pc.a[c1[k]]) * np.clip(t[k], 0, 1)
            first = int(pc.part[c1[k]]) < int(pc.part[c2[k]])
            ci, ti, cj, tj = (c1[k], t[k], c2[k], u[k]) if first else (c2[k], u[k], c1[k], t[k])
            i, j = int(pc.part[ci]), int(pc.part[cj])
            pos_i, pos_j = pc.position(ci, float(np.clip(ti, 0, 1))), pc.position(cj, float(np.clip(tj, 0, 1)))
            for part, (seg, f) in ((i, pos_i), (j, pos_j)):
                if f == 0.0:
                    x = self.parts[part][seg]  # on a vertex: both lines take its exact coordinate
            found.setdefault((i, j), []).append((pos_i, pos_j, tuple(x)))
        for (i, j), hits in sorted(found.items()):
            nodes = np.array(self.part_nodes(i) + self.part_nodes(j))
            kept = []
            for pos_i, pos_j, x in sorted(hits):
                xa = np.array(x)
                if np.hypot(*local(xa, nodes).T).min() <= SAME_POINT_M:
                    continue  # already a node both lines share (a shared vertex or a snap)
                if any(np.hypot(*local(xa, np.array(y[2]))) <= SAME_POINT_M for y in kept):
                    continue  # the same crossing found from the next piece
                kept.append((pos_i, pos_j, x))
            if not kept:
                continue
            noded = self.meta[i]['src'] != self.meta[j]['src']
            along = self.along(i, [k[0] for k in kept])
            dup = any(sum(1 for b in along if a <= b < a + DUP_WINDOW_M) > DUP_COUNT for a in along)
            for pos_i, pos_j, x in kept:
                if noded:
                    self.cuts[i][pos_i] = x
                    self.cuts[j][pos_j] = x
                self.crossings.append((i, j, x, noded, dup))
                self.qa['crossings_noded' if noded else 'crossings_unnoded'] += 1

    def along(self, i, positions):
        p = self.parts[i]
        cum = np.concatenate([[0], np.cumsum(haversine(p[:-1], p[1:]))])
        return [float(cum[s] + (cum[min(s + 1, len(p) - 1)] - cum[s]) * f) for s, f in positions]

    def edges(self):
        """Every part cut at its nodes: (part, coordinates) in sort order."""
        for i, p in enumerate(self.parts):
            p = p.copy()
            p[0], p[-1] = self.end_pos(2 * i), self.end_pos(2 * i + 1)
            pts = [(j, 0.0, tuple(p[j])) for j in range(len(p))]
            pts += [(seg, f, xy) for (seg, f), xy in self.cuts[i].items() if f > 0]
            pts.sort(key=lambda x: (x[0], x[1]))
            cut = set(self.cuts[i])
            cur = []
            for n, (seg, f, xy) in enumerate(pts):
                cur.append(xy)
                if (n == len(pts) - 1 or (seg, f) in cut) and len(cur) > 1:
                    yield i, cur
                    cur = [xy]


def build(parts, meta, snap=GRAPH_SNAP_V1):
    """Nodes the parts. Returns (edge features, gap features, crossing features, counts)."""
    g = Graph(parts, meta, snap)
    if parts:
        g.node()
    canon, raw, dropped = {}, [], 0
    for i, coords in g.edges():
        ll = [key(c) for c in coords]
        ll = [c for n, c in enumerate(ll) if n == 0 or c != ll[n - 1]]
        if len(ll) < 2:
            dropped += 1
            continue
        a, b = node_id(*ll[0]), node_id(*ll[-1])
        canon.setdefault(a, ll[0])  # a node's coordinate: the first member's, in sort order
        canon.setdefault(b, ll[-1])
        raw.append((i, a, b, ll))
    edges, seen = [], set()
    for i, a, b, ll in raw:
        ll = [canon[a], *ll[1:-1], canon[b]]
        ll = [c for n, c in enumerate(ll) if n == 0 or c != ll[n - 1]]
        arr = np.array(ll)
        length = round(float(haversine(arr[:-1], arr[1:]).sum()), 1) if len(ll) > 1 else 0.0
        m = meta[i]
        dedupe = (a, b, m['src'], m['sid'], tuple(ll))
        if len(ll) < 2 or length == 0 or dedupe in seen:
            dropped += 1
            continue
        seen.add(dedupe)
        edges.append({'type': 'Feature', 'properties': {'a': a, 'b': b, 'src': m['src'], 'sid': m['sid'], 'rid': m['rid'], 'len': length},
                      'geometry': {'type': 'LineString', 'coordinates': [list(c) for c in ll]}})
    edges.sort(key=lambda f: (f['properties']['a'], f['properties']['b'], f['properties']['src'] or '', f['properties']['sid'] or '',
                              str(f['properties']['rid']), json.dumps(f['geometry']['coordinates'])))

    def who(i):
        m = meta[i]
        return {'src': m['src'], 'sid': m['sid'], 'num': m['num'], 'name': m['name']}

    def label(i):
        m = meta[i]
        return m['num'] or m['name'] or m['sid']

    gaps = [{'type': 'Feature',
             'properties': {'gap_m': round(d, 1), 'from': who(e // 2), 'to': who(o),
                            'note': f'gap {round(d)} m between {label(e // 2)} and {label(o)}'},
             'geometry': {'type': 'Point', 'coordinates': list(key(g.end_pos(e)))}} for e, d, o in g.gaps]
    crossings = [{'type': 'Feature',
                  'properties': {'noded': noded, 'duplicate_drawing': dup, 'a': who(i), 'b': who(j)},
                  'geometry': {'type': 'Point', 'coordinates': list(key(x))}} for i, j, x, noded, dup in g.crossings]
    nodes = {n for f in edges for n in (f['properties']['a'], f['properties']['b'])}
    counts = {'parts': len(parts), 'nodes': len(nodes), 'edges': len(edges), 'gaps': len(gaps), 'edges_dropped': dropped,
              'duplicate_drawing_crossings': sum(1 for c in g.crossings if c[4]), **g.qa}
    return edges, gaps, crossings, counts


def write_seq(path, features):
    with open(path, 'w', encoding='utf-8', newline='\n') as out:
        for f in features:
            out.write(json.dumps(f, separators=(',', ':'), ensure_ascii=False) + '\n')


def credits(edges, sources_path, terrain):
    """credits.json: the credit of every route source the graph's edges come from, then the terrain's."""
    import terrarium
    with open(sources_path, encoding='utf-8') as f:
        known = {c['src']: c for c in json.load(f)}
    used = sorted({e['properties']['src'] for e in edges})
    missing = [s for s in used if not (known.get(s) or {}).get('attribution')]
    if missing:
        raise BuildError(f"no credit recorded in geo.dataset_imports for {', '.join(missing)}: re-import before building")
    out = {'sources': [{'src': s, 'attribution': known[s]['attribution'], 'license': known[s].get('license')} for s in used]}
    if terrain:
        out['terrain'] = terrarium.LAND_CREDIT
    out['attribution'] = '; '.join(dict.fromkeys([c['attribution'] for c in out['sources']] + [out.get('terrain')] * terrain))
    return out