"""MVP boundary-radial persistence loss.

The general routine computes capped H0 persistence on a polygonal boundary
 graph. Zero-length tie bars are discarded, as standard persistence diagrams
 discard diagonal points. The toy circle helper uses one canonical interval per
 connected circular boundary component, making the mechanism sweeps transparent.
Full XRPH relative/H1 computation and differentiable gradients are not included.
"""
from dataclasses import dataclass
import numpy as np
from scipy.optimize import linear_sum_assignment

@dataclass(frozen=True)
class Bar:
    birth: float
    death: float
    weight: float = 1.0

def boundary_graph_bars(points, components, tol=1e-10):
    points = np.asarray(points, dtype=float)
    if points.size == 0:
        return []
    points = points.reshape((-1, 2))
    r = np.linalg.norm(points, axis=1)
    edges = []
    for comp in components:
        comp = list(comp)
        for a, b in zip(comp, comp[1:] + comp[:1]):
            edges.append((max(r[a], r[b]), a, b))
    events = [(r[i], 0, i, i) for i in range(len(points))]
    events += [(f, 1, a, b) for f, a, b in edges]
    events.sort(key=lambda z: (z[0], z[1]))
    parent = list(range(len(points))); birth = r.copy()
    active = np.zeros(len(points), dtype=bool); bars = []
    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]; x = parent[x]
        return x
    for f, typ, a, b in events:
        if typ == 0:
            active[a] = True
        elif active[a] and active[b]:
            ra, rb = find(a), find(b)
            if ra != rb:
                keep, kill = (ra, rb) if birth[ra] <= birth[rb] else (rb, ra)
                if f - birth[kill] > tol:
                    bars.append(Bar(float(birth[kill]), float(f)))
                parent[kill] = keep
    maxf = float(max(r))
    roots = {find(i) for i in range(len(points)) if active[i]}
    bars.extend(Bar(float(birth[q]), maxf) for q in roots)
    return sorted(bars, key=lambda x: (x.birth, x.death))

def match_loss(pred, target, unmatched=1.0):
    """Hungarian L1 endpoint matching plus unmatched penalty."""
    n, m = len(pred), len(target)
    if n == 0 and m == 0: return 0.0, []
    k = n + m
    cost = np.full((k, k), float(unmatched))
    if n and m:
        cost[:n, :m] = [[abs(p.birth-t.birth)+abs(p.death-t.death)
                          for t in target] for p in pred]
    cost[n:, m:] = 0.0
    rows, cols = linear_sum_assignment(cost)
    return float(cost[rows, cols].sum()), list(zip(rows.tolist(), cols.tolist()))

def make_components(radii, n=24):
    pts, comps = [], []
    for rad in radii:
        start = len(pts)
        pts.extend([[rad*np.cos(2*np.pi*j/n), rad*np.sin(2*np.pi*j/n)]
                    for j in range(n)])
        comps.append(list(range(start, start+n)))
    return np.asarray(pts, dtype=float).reshape((-1, 2)), comps

def bars_from_radii(radii):
    """Canonical radial boundary-component bars used by the toy experiment."""
    return [Bar(float(rad), float(rad)) for rad in radii]

def radial_loss(pred_radii, target_radii, unmatched=1.0):
    return match_loss(bars_from_radii(pred_radii), bars_from_radii(target_radii), unmatched)[0]