#!/usr/bin/env python3
# equivalence-test --- the two-tier equivalence test for the BSPPFit
# mono-file model.
#
# BSPPFit converged (iteration 0.98.0) on a mono-file model: a single
# bsppfit.scad that produces all eleven trade categories, the
# dedicated per-category models (Pipe.scad, Cap.scad, ...) being kept as
# the reference ORACLE. This harness checks that the mono-file reproduces
# the oracle. It implements the two-tier method agreed with the project
# author:
#
#   Tier 1 --- functional-dimension equality. The two implementations are
#     NOT byte-identical meshes: the dedicated models are purpose-built
#     solids, while bsppfit builds a hub with arms, so wall and hub
#     geometry differ legitimately. What must match is the FUNCTIONAL
#     geometry --- the dimensions that decide whether parts mate: the male
#     thread major diameter, the bore diameter, and (where measurable) the
#     female socket minor diameter. These are extracted from rendered STLs
#     by axis-radius analysis and compared within a tolerance.
#
#   Tier 2 --- visual confirmation. A side-by-side still of the oracle and
#     the mono-file equivalent, for the eye to confirm topology and form.
#
# A CASE maps one oracle invocation to its mono-file equivalent and lists
# the functional dimensions to compare. The harness renders both to STL,
# measures, compares, and (optionally) renders the stills.
#
# Usage:
#   equivalence-test [--case NAME ...] [--tol MM] [--no-visual]
#                    [--outdir DIR] [--models-dir DIR]
#
# With no --case, runs every registered case. Exit status is non-zero if
# any case fails Tier 1.
#
# This file is part of BSPPFit. GNU GPL v3 or later.

import argparse
import math
import os
import re
import struct
import subprocess
import sys
import tempfile

# --- Environment -----------------------------------------------------

MODELS_DIR_DEFAULT = os.environ.get(
    "BSPPFIT_MODELS_DIR",
    os.path.dirname(os.path.abspath(__file__)))
OPENSCAD = os.environ.get("OPENSCAD", "openscad")
FN = int(os.environ.get("BSPPFIT_TEST_FN", "64"))
RENDER_TIMEOUT = int(os.environ.get("BSPPFIT_TEST_TIMEOUT", "400"))

# --- A test case -----------------------------------------------------
#
# Each case names the oracle module call and the mono-file call, plus the
# functional dimensions expected to agree. `dims` lists checks; each is a
# (label, kind, args) tuple where kind selects a measurement:
#   "thread_major"  --- max in-plane radius in a +Z or -Z tip band -> Ø
#   "bore"          --- global min in-plane radius -> Ø
# Both implementations are measured the same way and compared.

class Case:
    def __init__(self, name, oracle_call, mono_call, dims, header=""):
        self.name = name
        self.oracle_call = oracle_call
        self.mono_call = mono_call
        self.dims = dims          # list of (label, kind, args)
        self.header = header      # extra SCAD (e.g. use<> beyond defaults)


# The harness wraps each call in a standard preamble.
PREAMBLE = """include <BOSL2/std.scad>
include <BOSL2/threading.scad>
use <{model}>
$fn={fn};
"""


def scad_source(model, fn, call, header=""):
    return PREAMBLE.format(model=model, fn=fn) + header + "\n" + call + ";\n"


def render_stl(models_dir, model, call, header, outpath):
    """Render a single module call to an ASCII STL."""
    src = scad_source(model, FN, call, header)
    with tempfile.NamedTemporaryFile("w", suffix=".scad", delete=False,
                                     dir=models_dir) as f:
        f.write(src)
        scadfile = f.name
    try:
        env = dict(os.environ)
        env["OPENSCADPATH"] = (os.environ.get("OPENSCADPATH", "")
                               + os.pathsep + models_dir)
        r = subprocess.run(
            [OPENSCAD, "-o", outpath, "--export-format=asciistl", scadfile],
            capture_output=True, text=True, timeout=RENDER_TIMEOUT, env=env)
        ok = (r.returncode == 0 and os.path.exists(outpath)
              and os.path.getsize(outpath) > 0)
        return ok, r.stderr
    finally:
        os.unlink(scadfile)


def read_vertices(stl):
    """Return list of (x,y,z) from an ASCII STL."""
    v = []
    rx = re.compile(r"\s*vertex\s+(\S+)\s+(\S+)\s+(\S+)")
    with open(stl) as fh:
        for line in fh:
            m = rx.match(line)
            if m:
                v.append(tuple(map(float, m.groups())))
    return v


def measure(kind, args, verts):
    """Measure a functional dimension from vertices. Returns a diameter."""
    if not verts:
        return None
    zs = [z for _, _, z in verts]
    zmin, zmax = min(zs), max(zs)
    if kind == "thread_major":
        # band: a slab of `depth` mm at the +Z (end="top") or -Z tip.
        end = args.get("end", "top")
        depth = args.get("depth", 7.0)
        if end == "top":
            band = [math.hypot(x, y) for x, y, z in verts if z > zmax - depth]
        else:
            band = [math.hypot(x, y) for x, y, z in verts if z < zmin + depth]
        return 2 * max(band) if band else None
    if kind == "thread_major_axis":
        # Thread major Ø at the tip of a port pointing along unit vector
        # `axis` (any direction). Project every vertex onto the axis, take
        # the slab `depth` mm below the furthest projection (the tip),
        # and measure the maximum radius PERPENDICULAR to the axis there.
        # This compares a thread tip regardless of where the port points,
        # so it works for elbows and other non-axial fittings.
        ax = args["axis"]
        n = math.sqrt(sum(c * c for c in ax))
        u = [c / n for c in ax]
        depth = args.get("depth", 6.0)
        proj = [(x * u[0] + y * u[1] + z * u[2], (x, y, z)) for x, y, z in verts]
        pmax = max(p for p, _ in proj)
        radii = []
        for p, (x, y, z) in proj:
            if p > pmax - depth:
                # perpendicular distance from the axis line through origin
                dot = x * u[0] + y * u[1] + z * u[2]
                px, py, pz = x - dot * u[0], y - dot * u[1], z - dot * u[2]
                radii.append(math.sqrt(px * px + py * py + pz * pz))
        return 2 * max(radii) if radii else None
    if kind == "bore":
        # the through bore: smallest in-plane radius over the whole body.
        r = [math.hypot(x, y) for x, y, z in verts]
        return 2 * min(r)
    if kind == "socket_minor":
        # Female socket minor diameter: the smallest in-plane radius in a
        # band at the socket-mouth end, EXCLUDING the deep interior (a
        # blind socket's closed floor, or a dome) which would read as a
        # near-zero radius. We take the band [zc-depth, zc] near the mouth
        # and report the minor radius of the thread crest there, found as
        # the smallest radius that still lies on the socket wall (above a
        # floor fraction of the nominal bore).
        end = args.get("end", "bottom")
        depth = args.get("depth", 6.0)
        floor = args.get("floor", 3.0)   # ignore radii below this (axis/dome)
        if end == "top":
            band = [math.hypot(x, y) for x, y, z in verts if z > zmax - depth]
        else:
            band = [math.hypot(x, y) for x, y, z in verts if z < zmin + depth]
        wall = [r for r in band if r > floor]
        return 2 * min(wall) if wall else None
    raise ValueError("unknown measurement kind: " + kind)


# --- Case registry ---------------------------------------------------

def registered_cases():
    cases = []

    # --- Straight family (G-01) ---
    # Pipe: oracle Pipe(size=1/2) (default M x M), mono bsppfit_pipe_ports.
    cases.append(Case(
        name="pipe-g12-mm",
        oracle_call='Pipe(size = 1/2)',
        mono_call='bsppfit_pipe(size = 1/2)',
        dims=[
            ("male thread major Ø (+Z tip)", "thread_major", {"end": "top"}),
            ("male thread major Ø (-Z tip)", "thread_major", {"end": "bottom"}),
            ("bore Ø", "bore", {}),
        ],
    ))
    # Nipple: a short M x M connector; same functional thread/bore as the
    # pipe, so the mono equivalent is the same pipe port list.
    cases.append(Case(
        name="nipple-g12-mm",
        oracle_call='Nipple(size = 1/2)',
        mono_call='bsppfit_nipple(size = 1/2)',
        dims=[
            ("male thread major Ø (+Z tip)", "thread_major", {"end": "top"}),
            ("male thread major Ø (-Z tip)", "thread_major", {"end": "bottom"}),
            ("bore Ø", "bore", {}),
        ],
    ))

    # --- Cap (G-06) ---
    # The dedicated Cap has a wrench HEAD (hex by default); bsppfit_cap_ports
    # is the cap's mate feature alone --- one OPEN female port, a tapped
    # socket with a solid back. (A "closed" port is now a thread-less
    # integral seal, not a cap, so the reproduction uses an open port.) The
    # FUNCTIONAL dimension --- the female socket minor Ø that decides whether
    # the cap threads onto a pipe --- must match; the head is an un-migrated
    # decoration, checked separately below, not a mate dimension.
    cases.append(Case(
        name="cap-g12",
        oracle_call='Cap(size = 1/2)',
        mono_call='bsppfit(ports = bsppfit_cap_ports(size = 1/2, gender = "F"))',
        dims=[
            ("socket minor Ø (mouth)", "socket_minor", {"end": "bottom"}),
        ],
    ))

    # --- Plug (G-06) ---
    # Dedicated Plug also has a wrench head; the male thread is the mate
    # dimension. bsppfit_plug_ports is the plug's mate feature alone --- one
    # OPEN male port (its solid body and head are the bsppfit_plug
    # decoration).
    cases.append(Case(
        name="plug-g12",
        oracle_call='Plug(size = 1/2)',
        mono_call='bsppfit(ports = bsppfit_plug_ports(size = 1/2))',
        dims=[
            ("male thread major Ø (tip)", "thread_major", {"end": "bottom"}),
        ],
    ))

    # --- Head-decorated cap and plug (G-06) ---
    # With the head decoration completed (iteration 0.102.0), the mono-file
    # bsppfit_cap / bsppfit_plug reproduce the dedicated Cap / Plug
    # including the hex wrench head, not only the mate dimension. The hex
    # head reads as an across-corners span; we check the mate dimension
    # (socket minor / thread major) and the head across-corners Ø.
    cases.append(Case(
        name="caphead-g12",
        oracle_call='Cap(size = 1/2)',
        mono_call='bsppfit_cap(size = 1/2)',
        dims=[
            ("socket minor Ø (mouth)", "socket_minor", {"end": "bottom"}),
            ("hex head across-corners Ø", "thread_major", {"end": "top"}),
        ],
    ))
    cases.append(Case(
        name="plughead-g12",
        oracle_call='Plug(size = 1/2)',
        mono_call='bsppfit_plug(size = 1/2)',
        dims=[
            ("male thread major Ø (tip)", "thread_major", {"end": "bottom"}),
            ("hex head across-corners Ø", "thread_major", {"end": "top"}),
        ],
    ))

    # --- Elbow (G-02), 90 degrees, M x M ---
    # The swept mono elbow now poses exactly like the oracle (both arms run
    # +X and +Z), so each port's thread Ø is measured along the shared axis.
    # The functional dimension --- the male thread major Ø at each mouth ---
    # is what must agree.
    cases.append(Case(
        name="elbow-g12-mm",
        oracle_call='Elbow(size = 1/2)',
        mono_call='bsppfit_elbow(size = 1/2)',
        dims=[
            ("male thread major Ø (arm A)", "thread_major_axis",
             {"axis": [1, 0, 0]}),                    # +X arm
            ("male thread major Ø (arm B)", "thread_major_axis",
             {"axis": [0, 0, 1]}),                    # +Z arm
        ],
    ))

    # --- Reducer (G-05), 1/2 -> 3/8, M x M ---
    cases.append(Case(
        name="reducer-g12-g38-mm",
        oracle_call='Reducer(size_large = 1/2, size_small = 3/8)',
        mono_call='bsppfit_reducer(size_large = 1/2, size_small = 3/8)',
        dims=[
            ("large male thread major Ø", "thread_major", {"end": "bottom"}),
            ("small male thread major Ø", "thread_major", {"end": "top"}),
        ],
    ))

    # --- Tee (G-03), default female ---
    # The dedicated Tee is a one-branch Manifold; its functional ports are
    # the female sockets. Measure the bore (the cavity is open on all
    # axes) and a socket minor at the +Z run end.
    cases.append(Case(
        name="tee-g12-f",
        oracle_call='Manifold(size_run = 1/2, size_branch = 1/2, branch_count = 1, gender_run = "F", gender_branch = "F")',
        mono_call='bsppfit_tee(size_run = 1/2, size_branch = 1/2)',
        dims=[
            ("run socket minor Ø (+Z)", "socket_minor", {"end": "top"}),
            ("bore Ø", "bore", {}),
        ],
    ))

    # --- Cross (G-04), default female (two-sided one-branch Manifold) ---
    cases.append(Case(
        name="cross-g12-f",
        oracle_call='Manifold(size_run = 1/2, size_branch = 1/2, branch_count = 1, branch_sides = "both", gender_run = "F", gender_branch = "F")',
        mono_call='bsppfit_cross(size_run = 1/2, size_branch = 1/2)',
        dims=[
            ("run socket minor Ø (+Z)", "socket_minor", {"end": "top"}),
            ("bore Ø", "bore", {}),
        ],
    ))

    # --- Decoration sub-classes: oracle vs mono sub-class module ---
    # Hex bush (G-08): default oracle 1/2 male x 1/4 female; mono mirrors.
    cases.append(Case(
        name="hexbush-g12-g14",
        oracle_call='HexBush(size_male = 1/2, size_female = 1/4)',
        mono_call='bsppfit_hexbush(size_male = 1/2, size_female = 1/4)',
        dims=[
            ("male thread major Ø (tip)", "thread_major", {"end": "top"}),
            ("bore Ø", "bore", {}),
        ],
    ))
    # Hose tail (G-09): default oracle 1/2 male thread + barbs; mono mirrors.
    cases.append(Case(
        name="hosetail-g12",
        oracle_call='HoseTail(size = 1/2)',
        mono_call='bsppfit_hosetail(size = 1/2, gender = "M")',
        dims=[
            ("male thread major Ø (thread end)", "thread_major", {"end": "bottom"}),
            ("bore Ø", "bore", {}),
        ],
    ))
    # Union (G-07): 1/2 female-female; bore is gender-independent.
    cases.append(Case(
        name="union-g12-f",
        oracle_call='Union(size = 1/2, gender_a = "F", gender_b = "F")',
        mono_call='bsppfit_union(size = 1/2, gender_a = "F", gender_b = "F")',
        dims=[
            ("bore Ø", "bore", {}),
        ],
    ))

    return cases


CASE_MODEL = {
    # case-name-prefix -> oracle model file
    "pipe": "Pipe.scad",
    "nipple": "Nipple.scad",
    "cap": "Cap.scad",
    "caphead": "Cap.scad",
    "plug": "Plug.scad",
    "plughead": "Plug.scad",
    "elbow": "Elbow.scad",
    "reducer": "Reducer.scad",
    "tee": "Manifold.scad",
    "cross": "Manifold.scad",
    "hexbush": "HexBush.scad",
    "hosetail": "HoseTail.scad",
    "union": "Union.scad",
}


def oracle_model_for(case):
    for prefix, model in CASE_MODEL.items():
        if case.name.startswith(prefix):
            return model
    raise ValueError("no oracle model registered for case " + case.name)


# --- Runner ----------------------------------------------------------

def render_still(models_dir, model, call, header, outpath):
    """Render a single module call to a PNG still (3/4 view)."""
    src = scad_source(model, FN, call, header)
    with tempfile.NamedTemporaryFile("w", suffix=".scad", delete=False,
                                     dir=models_dir) as f:
        f.write(src)
        scadfile = f.name
    try:
        env = dict(os.environ)
        env["OPENSCADPATH"] = (os.environ.get("OPENSCADPATH", "")
                               + os.pathsep + models_dir)
        r = subprocess.run(
            [OPENSCAD, "-o", outpath, "--imgsize=600,600",
             "--camera=0,0,0,62,0,25,180", "--autocenter", "--viewall",
             scadfile],
            capture_output=True, text=True, timeout=RENDER_TIMEOUT, env=env)
        return (r.returncode == 0 and os.path.exists(outpath)), r.stderr
    finally:
        os.unlink(scadfile)


def run_case(case, models_dir, tol, outdir, do_visual):
    oracle_model = oracle_model_for(case)
    mono_model = "bsppfit.scad"
    o_stl = os.path.join(outdir, case.name + "-oracle.stl")
    m_stl = os.path.join(outdir, case.name + "-mono.stl")

    ok_o, err_o = render_stl(models_dir, oracle_model, case.oracle_call,
                             case.header, o_stl)
    ok_m, err_m = render_stl(models_dir, mono_model, case.mono_call,
                             case.header, m_stl)
    if not ok_o:
        return False, ["oracle render failed: " + err_o.strip()[:200]]
    if not ok_m:
        return False, ["mono render failed: " + err_m.strip()[:200]]

    vo = read_vertices(o_stl)
    vm = read_vertices(m_stl)
    lines = []
    passed = True
    for label, kind, args in case.dims:
        # Allow a dim to specify different axes for oracle and mono (used
        # when the two implementations pose a fitting differently).
        args_o = dict(args)
        args_m = dict(args)
        if "axis_oracle" in args:
            args_o["axis"] = args.pop("axis_oracle") if False else args["axis_oracle"]
            args_m["axis"] = args["axis_mono"]
        do = measure(kind, args_o, vo)
        dm = measure(kind, args_m, vm)
        if do is None or dm is None:
            lines.append("  %-32s MEASURE FAILED" % label)
            passed = False
            continue
        diff = abs(do - dm)
        ok = diff <= tol
        passed = passed and ok
        lines.append("  %-32s oracle=%7.3f  mono=%7.3f  Δ=%6.3f  %s"
                     % (label, do, dm, diff, "OK" if ok else "FAIL"))

    if do_visual:
        o_png = os.path.join(outdir, case.name + "-oracle.png")
        m_png = os.path.join(outdir, case.name + "-mono.png")
        vo_ok, _ = render_still(models_dir, oracle_model, case.oracle_call,
                                case.header, o_png)
        vm_ok, _ = render_still(models_dir, mono_model, case.mono_call,
                                case.header, m_png)
        lines.append("  visual stills: oracle=%s mono=%s"
                     % ("ok" if vo_ok else "FAILED",
                        "ok" if vm_ok else "FAILED"))

    return passed, lines


def main():
    ap = argparse.ArgumentParser(description="BSPPFit two-tier equivalence test")
    ap.add_argument("--case", action="append", default=None,
                    help="run only the named case(s)")
    ap.add_argument("--tol", type=float, default=0.05,
                    help="functional-dimension tolerance in mm (default 0.05)")
    ap.add_argument("--no-visual", action="store_true",
                    help="skip the Tier-2 visual stills")
    ap.add_argument("--outdir", default=None,
                    help="output directory (default: a temp dir)")
    ap.add_argument("--models-dir", default=MODELS_DIR_DEFAULT)
    args = ap.parse_args()

    outdir = args.outdir or tempfile.mkdtemp(prefix="bsppfit-equiv-")
    os.makedirs(outdir, exist_ok=True)

    # The per-fitting oracle .scad files this test compares bsppfit.scad's
    # mono-file wrappers against were retired at iteration 0.143.0, once
    # equivalence was confirmed passing for all ten topologies (see the
    # ResearchJournal entry of that date). Detect that state directly
    # (one representative oracle file is enough, since all ten retire
    # together) and say so plainly, rather than let every case fail with
    # a confusing "file not found" from deep inside the OpenSCAD run.
    if not os.path.isfile(os.path.join(args.models_dir, "Nipple.scad")):
        print("equivalence-test: the oracle .scad files this test compares",
              file=sys.stderr)
        print("  against were retired at iteration 0.143.0 (all ten",
              file=sys.stderr)
        print("  topologies confirmed equivalent first). This test has",
              file=sys.stderr)
        print("  nothing left to compare bsppfit.scad's wrappers against",
              file=sys.stderr)
        print("  and is not runnable until/unless repurposed --- see the",
              file=sys.stderr)
        print("  2026-09-21 ResearchJournal entry for the full account.",
              file=sys.stderr)
        return 2

    cases = registered_cases()
    if args.case:
        cases = [c for c in cases if c.name in args.case]
        if not cases:
            print("no matching case", file=sys.stderr)
            return 2

    print("BSPPFit equivalence test --- tolerance %.3f mm, $fn=%d" % (args.tol, FN))
    print("models: %s" % args.models_dir)
    print("output: %s" % outdir)
    print()

    all_ok = True
    for c in cases:
        print("CASE %s" % c.name)
        passed, lines = run_case(c, args.models_dir, args.tol, outdir,
                                 not args.no_visual)
        for l in lines:
            print(l)
        print("  => %s" % ("PASS" if passed else "FAIL"))
        print()
        all_ok = all_ok and passed

    print("RESULT: %s" % ("all cases passed" if all_ok else "FAILURES present"))
    return 0 if all_ok else 1


if __name__ == "__main__":
    sys.exit(main())
