#!/usr/bin/env python3
"""Re-runs the evidence check for the post organize-digital-photos.

Lab: builds a synthetic 500-photo, 15-day, five-city trip (fixed seed), clusters it
with DBSCAN at fixed radii and with the parameters Viallo's own organizer picks,
and compares every number with dataset.json.
Pages: fetches every cited URL and reports whether the quoted sentence is still on it.

Exit 0 when every number matches and every quote is still on its page, 1 otherwise.
--recheck re-reads the pages and skips a lab step that took over a minute when recorded.
Needs Python 3 and numpy.
"""
import json
import os
import re
import sys
import time
import unicodedata
import urllib.request
from html.parser import HTMLParser

import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
UA = "Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/131.0.0.0 Safari/537.36"

# --- lab inputs -------------------------------------------------------------

SEED = 20261026
N_PHOTOS = 500
DAYS = 15
SIGMA_M = 250.0  # per-axis Gaussian spread of each photo around its district centre
CITY_NAMES = ["Prague", "Vienna", "Budapest", "Krakow", "Berlin"]
# Hand-placed approximate centres; only the spacing of districts affects the result.
CITY_CENTRES = [(50.0755, 14.4378), (48.2082, 16.3738), (47.4979, 19.0402),
                (50.0647, 19.9450), (52.5200, 13.4050)]
DISTRICT_OFFSETS_M = [0.0, 2200.0, 4500.0]  # one district per day, three days per city
DISTRICT_BEARING_DEG = [0.0, 90.0, 180.0]
FIXED_RADII_M = [500, 1000, 2000, 5000]
FIXED_MIN_POINTS = 4
BURST_SIZE = 5       # bursty timeline: 5 photos 30 s apart, bursts spread over 09:00-19:00
BURST_STEP_S = 30


def offset(lat, lon, metres, bearing_deg):
    b = np.radians(bearing_deg)
    return (lat + metres * np.cos(b) / 111320.0,
            lon + metres * np.sin(b) / (111320.0 * np.cos(np.radians(lat))))


def build_inputs():
    """lat, lon, district id, city id, and two timelines in ms: evenly spaced and bursty."""
    rng = np.random.default_rng(SEED)
    per_day = np.full(DAYS, N_PHOTOS // DAYS)
    per_day[: N_PHOTOS - per_day.sum()] += 1
    lat, lon, district, city, even, bursty = [], [], [], [], [], []
    for day in range(DAYS):
        c, d = day // 3, day % 3
        dlat, dlon = offset(*CITY_CENTRES[c], DISTRICT_OFFSETS_M[d], DISTRICT_BEARING_DEG[d])
        n = int(per_day[day])
        n_bursts = -(-n // BURST_SIZE)
        for i in range(n):
            lat.append(dlat + rng.normal(0, SIGMA_M) / 111320.0)
            lon.append(dlon + rng.normal(0, SIGMA_M) / (111320.0 * np.cos(np.radians(dlat))))
            start = day * 86400 + 9 * 3600
            even.append(start + int(i * (10 * 3600 / max(n - 1, 1))))
            b, k = divmod(i, BURST_SIZE)
            bursty.append(start + int(b * (10 * 3600 / max(n_bursts - 1, 1))) + k * BURST_STEP_S)
            district.append(day)
            city.append(c)
    return (np.array(lat), np.array(lon), np.array(district), np.array(city),
            np.array(even, dtype=float) * 1000, np.array(bursty, dtype=float) * 1000)


# --- ports of apps/api/src/albums/organizer/*.ts ------------------------------

def haversine_matrix(lat, lon):
    # geo.service.ts haversineMeters, R = 6371000
    la, lo = np.radians(lat)[:, None], np.radians(lon)[:, None]
    h = np.sin((la - la.T) / 2) ** 2 + np.cos(la) * np.cos(la.T) * np.sin((lo - lo.T) / 2) ** 2
    return 2 * 6371000 * np.arcsin(np.sqrt(np.clip(h, 0.0, 1.0)))


def dbscan(dist, t_ms, eps_m, eps_t_ms, min_points):
    """clustering.service.ts dbscan: a point counts itself among its neighbours."""
    nb = dist <= eps_m
    if eps_t_ms is not None:
        nb &= np.abs(t_ms[:, None] - t_ms[None, :]) <= eps_t_ms
    n = dist.shape[0]
    visited = np.zeros(n, bool)
    assigned = np.zeros(n, bool)
    clusters = []
    for i in range(n):
        if visited[i]:
            continue
        visited[i] = True
        neighbours = np.flatnonzero(nb[i])
        if neighbours.size < min_points:
            continue
        cluster = [i]
        assigned[i] = True
        queue = [int(x) for x in neighbours]
        in_queue = set(queue)
        qi = 0
        while qi < len(queue):
            m = queue[qi]
            qi += 1
            if not visited[m]:
                visited[m] = True
                mn = np.flatnonzero(nb[m])
                if mn.size >= min_points:
                    for x in mn:
                        x = int(x)
                        if x not in in_queue:
                            in_queue.add(x)
                            queue.append(x)
            if not assigned[m]:
                assigned[m] = True
                cluster.append(m)
        clusters.append(cluster)
    return clusters


def find_knee_eps(s):
    n = len(s)
    if n == 0:
        return 2000.0
    if n == 1 or s[-1] == s[0]:
        return float(s[0])
    best, knee = -1.0, n * 3 // 4
    for i in range(n):
        d = abs(i / (n - 1) - (s[i] - s[0]) / (s[-1] - s[0])) / 2 ** 0.5
        if d > best:
            best, knee = d, i
    return float(s[knee])


def adaptive_params(dist, t_ms):
    """clustering.service.ts computeAdaptiveParams."""
    n = dist.shape[0]
    hour = 3600_000
    if n < 5:
        return 2000, 4 * hour, max(2, n - 1)
    min_points = max(2, min(8, int(np.floor(np.log2(n)))))
    d = dist.copy()
    np.fill_diagonal(d, np.inf)
    d.sort(axis=1)
    kd = np.sort(d[:, min(min_points - 1, n - 2)])
    eps = max(100.0, min(10_000.0, find_knee_eps(kd)))
    gaps = np.diff(t_ms)
    gaps = np.sort(gaps[gaps > 0])
    eps_t = 4 * hour
    if gaps.size:
        eps_t = max(hour, min(6 * hour, gaps[gaps.size // 2] * 3))
    # Math.round in JS rounds .5 up
    return int(np.floor(eps + 0.5)), int(np.floor(eps_t + 0.5)), min_points


def merge_count(clusters, lat, lon, radius):
    """albums-organizer.service.ts: centroid per cluster, union-find within radius."""
    if not clusters:
        return 0
    c_lat = np.array([lat[c].mean() for c in clusters])
    c_lon = np.array([lon[c].mean() for c in clusters])
    d = haversine_matrix(c_lat, c_lon)
    parent = list(range(len(clusters)))

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

    for i in range(len(clusters)):
        for j in range(i + 1, len(clusters)):
            if d[i, j] <= radius:
                ri, rj = find(i), find(j)
                if ri != rj:
                    parent[rj] = ri
    return len({find(i) for i in range(len(clusters))})


def product_pipeline(lat, lon, t_ms):
    """Primary ST-DBSCAN, merge into places, then the 500 m secondary pass on what is left.
    Pass A (noise to the nearest visit) and name-based merging after geocoding are not
    modelled; the first needs the visits, the second a geocoder."""
    order = np.argsort(t_ms, kind="stable")
    lat, lon, t_ms = lat[order], lon[order], t_ms[order]
    dist = haversine_matrix(lat, lon)
    eps, eps_t, min_pts = adaptive_params(dist, t_ms)
    clusters = dbscan(dist, t_ms, eps, eps_t, min_pts)
    merge_radius = max(eps * 1.5, 200)
    places = merge_count(clusters, lat, lon, merge_radius)
    left = np.setdiff1d(np.arange(len(lat)), [i for c in clusters for i in c])
    secondary = dbscan(dist[np.ix_(left, left)], None, 500, None, 2) if left.size else []
    singletons = int(left.size - sum(len(c) for c in secondary))
    return {
        "eps_m": eps, "eps_time_min": eps_t / 60000, "min_points": min_pts,
        "visits": len(clusters), "merge_radius_m": merge_radius, "places_from_visits": places,
        "photos_outside_visits": int(left.size), "secondary_places_500m": len(secondary),
        "singletons": singletons,
    }


def measure():
    lat, lon, district, city, even, bursty = build_inputs()
    dist = haversine_matrix(lat, lon)
    out = {"photos": int(len(lat)), "districts": int(len(set(district.tolist()))),
           "cities": len(CITY_NAMES)}
    for r in FIXED_RADII_M:
        cl = dbscan(dist, None, r, None, FIXED_MIN_POINTS)
        out["fixed_%dm_groups" % r] = len(cl)
        out["fixed_%dm_unclustered" % r] = int(len(lat) - sum(len(c) for c in cl))
        out["fixed_%dm_groups_spanning_districts" % r] = sum(len(set(district[c].tolist())) > 1 for c in cl)
    for name, t in (("even", even), ("bursty", bursty)):
        for k, v in product_pipeline(lat, lon, t).items():
            out["viallo_%s_%s" % (name, k)] = v
    return out


# --- pages ------------------------------------------------------------------

class _Text(HTMLParser):
    def __init__(self):
        super().__init__(convert_charrefs=True)
        self.skip, self.out = 0, []

    def handle_starttag(self, tag, attrs):
        if tag in ("script", "style", "noscript", "svg"):
            self.skip += 1
        if tag in ("p", "div", "li", "br", "h1", "h2", "h3", "h4", "tr", "td", "section"):
            self.out.append("\n")

    def handle_endtag(self, tag):
        if tag in ("script", "style", "noscript", "svg") and self.skip:
            self.skip -= 1

    def handle_data(self, data):
        if not self.skip:
            self.out.append(data)


def norm(s):
    s = unicodedata.normalize("NFKC", s)
    for a, b in [("\u2018", "'"), ("\u2019", "'"), ("\u02bc", "'"), ("\u201c", '"'), ("\u201d", '"'),
                 ("\u2013", "-"), ("\u2014", "-"), ("\u00a0", " "), ("\u00ae", ""), ("\u2122", "")]:
        s = s.replace(a, b)
    return re.sub(r"\s+", " ", s).strip()


def fetch_text(url, tries=3):
    last = None
    for attempt in range(tries):
        try:
            req = urllib.request.Request(url, headers={
                "User-Agent": UA, "Accept": "text/html,application/xhtml+xml",
                "Accept-Language": "en-US,en;q=0.9"})
            with urllib.request.urlopen(req, timeout=60) as r:
                p = _Text()
                p.feed(r.read().decode("utf-8", errors="replace"))
                return norm("".join(p.out))
        except Exception as exc:
            last = exc
            time.sleep(2 + 2 * attempt)
    raise last


def main():
    recheck = "--recheck" in sys.argv
    with open(os.path.join(HERE, "dataset.json")) as fh:
        data = json.load(fh)
    failures = 0

    print("lab")
    expected = data["results"]["lab"]
    if recheck and data.get("labSeconds", 0) > 60:
        print("  skipped: recorded run took %ss" % data["labSeconds"])
    else:
        t0 = time.time()
        got = measure()
        for key, want in expected.items():
            ok = got.get(key) == want
            failures += not ok
            print("  %s %-38s want %-8s got %s" % ("ok  " if ok else "DRIFT", key, want, got.get(key)))
        print("  lab rebuilt in %.1fs" % (time.time() - t0))

    print("pages")
    cache = {}
    for row in data["results"]["pages"]:
        url = row["url"]
        if url not in cache:
            try:
                cache[url] = fetch_text(url)
            except Exception as exc:
                cache[url] = exc
        if isinstance(cache[url], Exception):
            failures += 1
            print("  NOT CHECKED %-14s %-28s %s (%s)" % (row["option"], row["property"], url, cache[url]))
            continue
        found = norm(row["quote"]) in cache[url]
        failures += not found
        print("  %s %-14s %-28s %s" % ("ok  " if found else "GONE", row["option"], row["property"], url))
        if not found:
            print("       quote: %s" % row["quote"])

    print("\n%d failure(s)" % failures)
    return 0 if failures == 0 else 1


if __name__ == "__main__":
    if "--measure" in sys.argv:
        print(json.dumps(measure(), indent=1))
        sys.exit(0)
    sys.exit(main())