"""Where are the points wrong beyond a band? For a cached descriptor: the far-wrong points' layers and parts, the answer
in every query mode and at the two neighbouring layer centres, so a query fault (modes or centres disagree with the
blend) is told apart from an encoding fault (the rings themselves are off).

    env -u PYTHONPATH .venv/bin/python scripts/diag_far.py <object> <cache .fcp> [band % diag]
"""
import sys
from collections import Counter
from pathlib import Path

import numpy as np
import open3d as o3d

sys.path.insert(0, str(Path(__file__).parent))
from common import CACHE, fcp, load_mesh  # noqa: E402
from fcpq import FCPQuery  # noqa: E402

LPAD = 1.1  # fcp._slice.L

name, f = sys.argv[1], sys.argv[2]
band = float(sys.argv[3]) if len(sys.argv) > 3 else 0.5
m = load_mesh(name)
sc = o3d.t.geometry.RaycastingScene()
sc.add_triangles(o3d.core.Tensor(np.asarray(m.vertices, np.float32)), o3d.core.Tensor(np.asarray(m.faces, np.uint32)))
np.random.seed(0)
rng = np.random.default_rng(0)
diag = np.linalg.norm(m.extents)
S, fi = m.sample(1_000_000, return_index=True)
P = S + m.face_normals[fi] * rng.uniform(-0.02, 0.02, (len(S), 1)) * diag
T = o3d.core.Tensor(P.astype(np.float32))
gt = sc.compute_occupancy(T, nsamples=3).numpy() > 0
dist = np.abs(sc.compute_distance(T).numpy()) / diag * 100
d = fcp.FCP.from_bytes((CACHE / f).read_bytes())
q = FCPQuery(d, nt=max(512, 8 * d.N))
Pc = d.pose.to_canonical(P).astype(np.float32)
ins = q.inside_cpu(Pc)
bad = np.flatnonzero((ins != gt) & (dist > band))
print(f"{f}: parts {q.P}, wrong beyond {band} % diag: {len(bad)} (false inside {int(ins[bad].sum())})")
X = Pc[bad]
lay = (X[:, 2] + LPAD) / (2 * LPAD) * d.layers - 0.5
for mode in ("scan",):
    print(f"  {mode} agrees with table on {int((q.inside_cpu(X, mode) == ins[bad]).sum())} / {len(bad)}")
c0, c1 = X.copy(), X.copy()
c0[:, 2] = -LPAD + (np.floor(lay) + 0.5) * 2 * LPAD / d.layers
c1[:, 2] = -LPAD + (np.floor(lay) + 1.5) * 2 * LPAD / d.layers
a0, a1 = q.inside_cpu(c0), q.inside_cpu(c1)
print(f"  at the lower centre right: {int((a0 == gt[bad]).sum())}, upper centre right: {int((a1 == gt[bad]).sum())}, "
      f"both right: {int(((a0 == gt[bad]) & (a1 == gt[bad])).sum())}")
L = np.floor(lay).astype(int)
print("  layers p0/p50/p100", np.percentile(lay, [0, 50, 100]).round(1), " top layers", Counter(L).most_common(8))
top = Counter(L).most_common(1)[0][0]
print(f"  parts spanning layer {top}:", [(int(q.sign[p]), int(q.h0[p]), int(q.h1[p])) for p in range(q.P) if q.h0[p] - 1 <= top <= q.h1[p] + 1])
np.save(CACHE / f"diag_{Path(f).stem}.npy", np.column_stack([P[bad], lay, ins[bad], gt[bad], dist[bad]]))