#!/usr/bin/env python3
"""Finite Lorentz transformations and the positive-shell measure, as an experiment.

Run: python lorentz-transformations.py --output-dir a-new-directory
Requires Python >=3.10 and NumPy; verified with Python 3.12.14 / NumPy 2.3.5.
No network access, installation, plotting package, or local helper is needed.

Conventions: c=hbar=1, eta=diag(1,-1,-1,-1), column four-vectors (t,x,y,z).
B(xi,n) is PASSIVE: the new frame moves with beta=tanh(xi)*n; its mixed
entries are -sinh(xi)*n. Rodrigues R(theta,n) is defined by positive component
rotation n cross v; pass -theta to describe axes turned by +theta. A@B acts
with B first. Signed rapidity is distinct from the unit direction.

Read the vector CSV as a finite experiment, not a proof of the Lorentz group.
Rotations and boosts generate SO+(1,3), not all four components of O(1,3).
P, T and PT are separately tabulated; classical T is not a quantum operator.
The measure experiment checks d^3p/(2E), not a scattering cross section or
the invariance of an independently imposed spherical momentum cutoff.

Exercises: reverse the boost signs and inspect the rest-particle check;
reverse two noncollinear boosts and compare their output; vary --fd-step
and compare all THREE finite-difference spacings; increase --max-rapidity
and compare absolute versus scaled metric errors and the Doppler reference.
Float64 matrix tests are restricted to |xi|<=12. A separate Decimal check
at xi=40 illustrates why cancellation can destroy the small photon energy.
Fraction arithmetic independently checks rational boosts, a rotation and
compositions, including the exact massive/massless shell derivative Jacobian.

References: Rindler, Introduction to Special Relativity (2nd ed., 1991);
Weinberg, The Quantum Theory of Fields, Vol. I (1995), sections 2.3 and 2.5.
Local physics owners: /relativistic-qm/lorentz-transformations/,
/relativistic-qm/four-vectors/, /relativistic-qm/relativistic-phase-space/.
"""

from __future__ import annotations

import argparse
import csv
from decimal import Decimal, localcontext
from fractions import Fraction
import hashlib
from itertools import permutations
import json
import math
from pathlib import Path
import platform
import subprocess
import sys

import numpy as np

ETA = np.diag([1.0, -1.0, -1.0, -1.0])
EPS = np.finfo(float).eps
ALGEBRA_TOL = 5e-13
FD_TOL = 2e-7
DECIMAL_TOL = 1e-40


def direction(values):
    """Normalize a nonzero finite direction; boosts are not parametrized by beta."""
    n = np.asarray(values, dtype=float)
    if n.shape != (3,) or not np.all(np.isfinite(n)):
        raise ValueError('direction must have three finite entries')
    norm = float(np.linalg.norm(n))
    if norm == 0 or not math.isfinite(norm):
        raise ValueError('direction must have a finite nonzero norm')
    return n / norm


def boost(rapidity, axis):
    """Return B for the passive frame velocity tanh(rapidity)*unit(axis)."""
    xi = float(rapidity)
    if not math.isfinite(xi) or abs(xi) > 12:
        raise ValueError('float64 boost rapidity must be finite and |xi|<=12')
    n = direction(axis)
    sh, ch = math.sinh(xi), math.cosh(xi)
    # 2*sinh(xi/2)^2 avoids subtracting 1 from cosh(xi) near xi=0.
    cm1 = 2 * math.sinh(xi / 2) ** 2
    result = np.eye(4)
    result[0, 0] = ch
    result[0, 1:] = result[1:, 0] = -sh * n
    result[1:, 1:] += cm1 * np.outer(n, n)
    return result


def rotation(angle, axis):
    """Positive right-hand component rotation; passive axes use minus angle."""
    if not math.isfinite(angle):
        raise ValueError('rotation angle must be finite')
    x, y, z = direction(axis)
    cross = np.array([[0, -z, y], [z, 0, -x], [-y, x, 0]])
    result = np.eye(4)
    result[1:, 1:] += math.sin(angle) * cross + 2 * math.sin(angle / 2) ** 2 * (cross @ cross)
    return result


def mdot(a, b):
    return float(np.asarray(a) @ ETA @ np.asarray(b))


def shell_momentum(momentum, mass):
    p = np.asarray(momentum, dtype=float)
    if p.shape != (3,) or not np.all(np.isfinite(p)) or mass < 0 or not math.isfinite(mass):
        raise ValueError('finite three-momentum and nonnegative finite mass required')
    energy = math.hypot(mass, *p)
    if energy == 0:
        raise ValueError('massless p=0 is the nonregular cone apex; excluded')
    return np.r_[energy, p]


def shell_jacobian_fd(matrix, momentum, mass, h):
    """Five-point derivative of p -> spatial part of L(sqrt(p^2+m^2),p).

    Every displaced point is put on shell again. This function intentionally
    uses neither the derivative dE/dp nor the predicted determinant E'/E.
    """
    if not math.isfinite(h) or h <= 0:
        raise ValueError('finite positive differentiation step required')
    p = np.asarray(momentum, dtype=float)
    jacobian = np.empty((3, 3))
    for j in range(3):
        shift = np.eye(3)[j] * h
        fm2, fm1, fp1, fp2 = [
            (matrix @ shell_momentum(p + k * shift, mass))[1:]
            for k in (-2, -1, 1, 2)
        ]
        jacobian[:, j] = (fm2 - 8 * fm1 + 8 * fp1 - fp2) / (12 * h)
    return jacobian


class Checks:
    def __init__(self):
        self.items = []

    def check(self, name, case, residual, tolerance=ALGEBRA_TOL):
        residual = abs(float(residual))
        self.items.append(dict(name=name, case=case, residual=residual,
                               tolerance=tolerance,
                               passed=math.isfinite(residual) and residual <= tolerance))

    def require(self, name, case, condition):
        self.check(name, case, 0 if condition else 1, 0)

    def summary(self):
        names = sorted({item['name'] for item in self.items})
        return {name: dict(count=len(rows := [r for r in self.items if r['name'] == name]),
                           max_residual=max(r['residual'] for r in rows),
                           tolerance=rows[0]['tolerance'],
                           passed=all(r['passed'] for r in rows)) for name in names}


def matrix_checks(checks, name, matrix):
    metric_abs = float(np.max(np.abs(matrix.T @ ETA @ matrix - ETA)))
    inverse_abs = float(np.max(np.abs((ETA @ matrix.T @ ETA) @ matrix - np.eye(4))))
    # Product scales describe cancellation, not relative accuracy of a tiny invariant.
    scale = max(1.0, float(np.max(np.abs(matrix).T @ np.abs(matrix))))
    det_error = abs(float(np.linalg.det(matrix)) - 1)
    # For a Lorentz matrix inv(L)=eta L.T eta; this estimate avoids an SVD.
    condition_estimate = float(np.linalg.norm(matrix, np.inf) * np.linalg.norm(matrix.T, np.inf))
    checks.check('metric_product_scaled', name, metric_abs / scale)
    checks.check('metric_inverse_scaled', name, inverse_abs / scale)
    checks.check('determinant_condition_scaled', name, det_error / max(1.0, condition_estimate))
    checks.require('proper_time_orientation', name, matrix[0, 0] >= 1 - ALGEBRA_TOL)
    return dict(case=name, metric_absolute_error=metric_abs, metric_scale=scale,
                metric_scaled_error=metric_abs / scale, inverse_absolute_error=inverse_abs,
                determinant=float(np.linalg.det(matrix)), determinant_absolute_error=det_error,
                condition_estimate=condition_estimate,
                estimated_roundoff_floor=EPS * condition_estimate)


def decimal_doppler(rapidity, checks):
    """An 100-digit independent exp-based boost compared with exp(-xi)."""
    with localcontext() as context:
        context.prec = 100
        xi = Decimal(str(rapidity))
        ep, em = xi.exp(), (-xi).exp()
        ch, sh = (ep + em) / 2, (ep - em) / 2
        energy = ch - sh
        rel = abs((energy - em) / em)
        checks.check('decimal_photon_relative', 'decimal-high-rapidity', float(rel), DECIMAL_TOL)
        # This float64 demonstration is reported, not asserted to retain relative accuracy.
        float_energy = math.cosh(rapidity) - math.sinh(rapidity)
        return dict(rapidity=rapidity, decimal_precision_digits=context.prec,
                    decimal_energy=str(energy), expected_exp_minus_xi=str(em),
                    decimal_relative_error=str(rel), float64_energy=float_energy,
                    float64_relative_error=abs(float_energy / float(em) - 1),
                    warning='Float64 subtraction can return zero or a wrong sign at large xi; use stable light-cone variables or higher precision.')


def exact_fixtures(checks):
    """Small rational examples independent of NumPy and the rapidity constructor.

    The shell derivative is analytic here: J_ij=L_ij+L_i0*p_j/E.
    The separate finite-difference experiment never uses this formula.
    Exact finite fixtures do not constitute a symbolic proof for arbitrary L.
    """
    F = Fraction

    def transpose(matrix):
        return [list(column) for column in zip(*matrix)]

    def multiply(a, b):
        return [[sum((x * y for x, y in zip(row, column)), F(0))
                 for column in zip(*b)] for row in a]

    def determinant(matrix):
        size = len(matrix)
        answer = F(0)
        for order in permutations(range(size)):
            inversions = sum(order[i] > order[j] for i in range(size) for j in range(i + 1, size))
            term = F((-1) ** inversions)
            for i in range(size):
                term *= matrix[i][order[i]]
            answer += term
        return answer

    def scalar(a, b):
        return sum((F(1 if i == 0 else -1) * a[i] * b[i] for i in range(4)), F(0))

    def transform(matrix, vector):
        return [row[0] for row in multiply(matrix, [[value] for value in vector])]

    identity = [[F(int(i == j)) for j in range(4)] for i in range(4)]
    eta = [[F((1 if i == 0 else -1) * int(i == j)) for j in range(4)] for i in range(4)]
    # Passive beta_x=3/5: gamma=5/4 and gamma*beta=3/4.
    bx = [[F(5, 4), F(-3, 4), F(0), F(0)],
          [F(-3, 4), F(5, 4), F(0), F(0)],
          [F(0), F(0), F(1), F(0)], [F(0), F(0), F(0), F(1)]]
    # Passive beta_y=-5/13: gamma=13/12 and -gamma*beta=5/12.
    by = [[F(13, 12), F(0), F(5, 12), F(0)],
          [F(0), F(1), F(0), F(0)],
          [F(5, 12), F(0), F(13, 12), F(0)],
          [F(0), F(0), F(0), F(1)]]
    # Positive component rotation about z, cos(theta)=3/5, sin(theta)=4/5.
    rz = [[F(1), F(0), F(0), F(0)],
          [F(0), F(3, 5), F(-4, 5), F(0)],
          [F(0), F(4, 5), F(3, 5), F(0)],
          [F(0), F(0), F(0), F(1)]]
    massive = [F(173, 115), F(72, 115), F(96, 115), F(48, 115)]
    massless = [F(1), F(2, 3), F(2, 3), F(1, 3)]
    event = [F(2), F(1, 3), F(-2, 7), F(5, 9)]
    checks.require('exact_massive_shell', 'rational-p', scalar(massive, massive) == 1)
    checks.require('exact_massless_shell', 'rational-p', scalar(massless, massless) == 0)
    fixture_rows = []
    for name, matrix in [('Bx', bx), ('Rz', rz), ('By-after-Bx', multiply(by, bx)),
                         ('Rz-after-By-after-Bx', multiply(rz, multiply(by, bx)))]:
        inverse = multiply(eta, multiply(transpose(matrix), eta))
        checks.require('exact_metric', name, multiply(transpose(matrix), multiply(eta, matrix)) == eta)
        checks.require('exact_inverse', name, multiply(inverse, matrix) == identity)
        checks.require('exact_determinant', name, determinant(matrix) == 1)
        checks.require('exact_scalar_product', name,
                       scalar(transform(matrix, massive), transform(matrix, event)) == scalar(massive, event))
        for shell_name, momentum in [('massive-m=1', massive), ('massless-nonzero', massless)]:
            mapped = transform(matrix, momentum)
            jacobian = [[matrix[i + 1][j + 1] + matrix[i + 1][0] * momentum[j + 1] / momentum[0]
                         for j in range(3)] for i in range(3)]
            det_j, energy_ratio = determinant(jacobian), mapped[0] / momentum[0]
            checks.require('exact_future_energy', f'{name}/{shell_name}', mapped[0] > 0)
            checks.require('exact_shell_derivative_jacobian', f'{name}/{shell_name}', det_j == energy_ratio)
            fixture_rows.append(dict(transformation=name, shell=shell_name,
                                     input_momentum=[str(value) for value in momentum],
                                     transformed_energy=str(mapped[0]),
                                     determinant_derivative=str(det_j), energy_ratio=str(energy_ratio),
                                     exact_equality=det_j == energy_ratio))
    twice = multiply(bx, bx)
    checks.require('exact_collinear_composition', 'two-beta-3/5-boosts',
                   twice[0][0] == F(17, 8) and twice[0][1] == F(-15, 8))
    checks.require('exact_passive_rest_sign', 'beta=3/5',
                   transform(bx, [F(1), F(0), F(0), F(0)]) == [F(5, 4), F(-3, 4), F(0), F(0)])
    return dict(arithmetic='Python Fraction, exact rational operations; no NumPy operations',
                determinant_method='Leibniz permutation sum, independent of numpy.linalg.det',
                scope='Four rational matrices and two shell points; finite exact fixtures, not a universal proof',
                shell_derivative='J_ij=L_(i+1)(j+1)+L_(i+1)0*p_(j+1)/E', rows=fixture_rows)


def run_experiment(args):
    checks = Checks()
    exact_report = exact_fixtures(checks)
    rng = np.random.default_rng(args.seed)
    x, y, oblique = np.eye(3)[0], np.eye(3)[1], direction([1, -2, 3])
    cases = [('identity', boost(0, x)), ('tiny-positive', boost(1e-10, oblique)),
             ('tiny-negative', boost(-1e-10, oblique)), ('passive-3over5', boost(math.log(2), x)),
             ('negative-oblique', boost(-1.1, oblique)), ('oblique', boost(0.9, oblique)),
             ('high-positive', boost(args.max_rapidity, oblique)),
             ('high-negative', boost(-args.max_rapidity, oblique)),
             ('rotation', rotation(0.7, oblique))]
    for i in range(args.samples):
        axis = direction(rng.normal(size=3))
        xi = rng.uniform(-min(args.max_rapidity, 3), min(args.max_rapidity, 3))
        cases.append((f'seeded-{i}', rotation(rng.uniform(-math.pi, math.pi), rng.normal(size=3)) @ boost(xi, axis)))

    b1, b2 = boost(0.9, x), boost(-0.7, y)
    r = rotation(0.4, oblique)
    cases.extend([('B2-after-B1', b2 @ b1), ('B1-after-B2', b1 @ b2), ('R-after-B2-after-B1', r @ b2 @ b1)])
    vector_rows, matrix_rows, composition_rows = [], [], []
    scalar_absolute_max = 0.0
    vectors = [('timelike-rest', np.array([1., 0, 0, 0])),
               ('timelike-moving', shell_momentum([.3, -.7, 1.1], 0.8)),
               ('null-axis', np.array([1., 1, 0, 0])),
               ('null-oblique', np.r_[1., oblique]),
               ('spacelike', np.array([.3, 1., -2, .5])),
               ('past-timelike', np.array([-1., 0, 0, 0]))]
    for name, matrix in cases:
        matrix_rows.append(matrix_checks(checks, name, matrix))
        for label, v in vectors:
            w = matrix @ v
            original, transformed = mdot(v, v), mdot(w, w)
            invariant_scale = max(1.0, float(w @ w), float(v @ v))
            error = abs(transformed - original)
            checks.check('vector_norm_product_scaled', f'{name}/{label}', error / invariant_scale)
            cone_violation = max(0., float(np.linalg.norm(w[1:])) - w[0])
            if label.startswith(('timelike', 'null')):
                checks.require('future_energy_positive', f'{name}/{label}', w[0] > 0)
                # Include cancellation in the transformation, as for scalar products.
                cone_scale = max(1., float(np.max(np.abs(matrix) @ np.abs(v))))
                checks.check('future_cone_scaled', f'{name}/{label}', cone_violation / cone_scale)
            vector_rows.append(dict(case=name, vector=label, **{f'input_{i}': float(v[i]) for i in range(4)},
                                    **{f'output_{i}': float(w[i]) for i in range(4)},
                                    norm_before=original, norm_after=transformed,
                                    invariant_absolute_error=error, invariant_scale=invariant_scale,
                                    invariant_scaled_error=error / invariant_scale,
                                    future_cone_absolute_violation=cone_violation if label.startswith(('timelike', 'null')) else 'not-applicable'))
        for ia, (_, a) in enumerate(vectors):
            for _, b in vectors[ia:]:
                ap, bp = matrix @ a, matrix @ b
                # Include cancellation INSIDE each matrix-vector multiplication.
                # For a redshifted null vector, |L a| is tiny but |L| |a| is large.
                # Scaling only by the already-computed vectors misses that error.
                scale = max(1., float((np.abs(matrix) @ np.abs(a)) @ (np.abs(matrix) @ np.abs(b))),
                            float(np.abs(a) @ np.abs(b)))
                error = abs(mdot(ap, bp) - mdot(a, b))
                scalar_absolute_max = max(scalar_absolute_max, error)
                checks.check('scalar_product_scaled', name, error / scale)

    # These sign and order tests distinguish a correct Lorentz matrix from the wrong convention.
    checks.check('rest_passive_sign', 'beta=3/5', np.max(np.abs(boost(math.log(2), x) @ vectors[0][1] - [1.25, -.75, 0, 0])))
    checks.check('rotation_orientation', 'z-pi/2', np.max(np.abs(rotation(math.pi / 2, [0, 0, 1]) @ [0, 1, 0, 0] - [0, 0, 1, 0])))
    for xi in [0, 1e-10, -1.3, args.max_rapidity]:
        b = boost(xi, oblique)
        scale = max(1., float(np.max(np.abs(b).T @ np.abs(b))))
        checks.check('opposite_rapidity_inverse_scaled', str(xi), np.max(np.abs(boost(-xi, oblique) @ b - np.eye(4))) / scale)
        origin_worldline = np.r_[1., math.tanh(xi) * oblique]
        checks.check('moving_origin_at_rest_scaled', str(xi), np.max(np.abs((b @ origin_worldline)[1:])) / max(1., float(np.max(np.abs(b)))))
    for first, second in [(.3, .8), (-.7, .2), (0., 0.), (1e-10, -1e-10)]:
        lhs, rhs = boost(second, oblique) @ boost(first, oblique), boost(first + second, oblique)
        error = float(np.max(np.abs(lhs - rhs)))
        beta_formula = (math.tanh(first) + math.tanh(second)) / (1 + math.tanh(first) * math.tanh(second))
        checks.check('collinear_rapidity_addition', f'{first}+{second}', error)
        checks.check('collinear_velocity_addition', f'{first}+{second}', beta_formula - math.tanh(first + second))
        composition_rows.append(dict(test='collinear', first_rapidity=first, second_rapidity=second,
                                     residual=error, expected='B(second) @ B(first) = B(first+second)'))
    v = np.array([2., .2, -.4, .7])
    sequential_error = float(np.max(np.abs((b2 @ b1) @ v - b2 @ (b1 @ v))))
    commutator_gap = float(np.max(np.abs(b2 @ b1 - b1 @ b2)))
    checks.check('rightmost_acts_first', 'noncollinear', sequential_error)
    checks.require('noncollinear_order_matters', 'noncollinear', commutator_gap > .1)
    composition_rows.extend([
        dict(test='rightmost-first', first_rapidity=.9, second_rapidity=-.7, residual=sequential_error, expected='B2 @ B1 @ v = B2 @ (B1 @ v)'),
        dict(test='noncommutation-gap', first_rapidity=.9, second_rapidity=-.7, residual=commutator_gap, expected='A nonzero gap is expected, not an error')])
    wrong_sign = float(np.max(np.abs(boost(-math.log(2), x) @ vectors[0][1] - [1.25, -.75, 0, 0])))
    wrong_inverse = float(np.max(np.abs(b1.T @ b1 - np.eye(4))))
    checks.require('negative_control_wrong_boost_sign_detected', 'passive-rest', wrong_sign > .5)
    checks.require('negative_control_transpose_inverse_detected', 'boost', wrong_inverse > .5)

    component_rows = []
    for name, diagonal in [('I', [1, 1, 1, 1]), ('P', [1, -1, -1, -1]),
                           ('T', [-1, 1, 1, 1]), ('PT', [-1, -1, -1, -1])]:
        matrix = np.diag(diagonal)
        det = int(round(np.linalg.det(matrix)))
        future = bool((matrix @ vectors[0][1])[0] > 0)
        checks.check('component_metric', name, np.max(np.abs(matrix.T @ ETA @ matrix - ETA)))
        component_rows.append(dict(component=name, determinant=det, time_time=int(matrix[0, 0]),
                                   proper=det == 1, orthochronous=future,
                                   positive_shell_preserved=future))

    # Isolate differentiation from high-boost cancellation: FD uses |xi|<=3.
    # This axis is independent of the much larger rapidity matrix stress test.
    measure_cases = [('massive-rest-identity', 1., [0, 0, 0], np.eye(4)),
                     ('massive-oblique', .8, [.3, -.7, 1.1], boost(.9, oblique)),
                     ('massless-parallel', 0., [1, 0, 0], boost(3., x)),
                     ('massless-antiparallel', 0., [1, 0, 0], boost(-3., x)),
                     ('massless-soft', 0., [1e-8, -2e-8, 3e-8], boost(.6, oblique)),
                     ('massless-rotation-composition', 0., [.3, -.7, 1.1], r @ b2 @ b1)]
    for i in range(args.samples):
        mass = 0. if i % 2 == 0 else rng.uniform(.1, 2)
        measure_cases.append((f'seeded-{i}', mass, rng.normal(size=3), boost(rng.uniform(-3, 3), rng.normal(size=3))))
    measure_rows = []
    for name, mass, p, matrix in measure_cases:
        q = shell_momentum(p, mass)
        qp = matrix @ q
        predicted = float(qp[0] / q[0])
        relative_errors = []
        determinants = []
        for divisor in (1, 2, 4):
            h = args.fd_step * q[0] / divisor
            determinant = float(np.linalg.det(shell_jacobian_fd(matrix, p, mass, h)))
            relative_error = abs(determinant / predicted - 1)
            checks.check('shell_measure_fd_relative', f'{name}/h-div-{divisor}', relative_error, FD_TOL)
            checks.require('shell_jacobian_positive', name, determinant > 0)
            relative_errors.append(relative_error)
            determinants.append(determinant)
            measure_rows.append(dict(case=name, mass=mass, px=float(p[0]), py=float(p[1]), pz=float(p[2]),
                                     energy=float(q[0]), transformed_energy=float(qp[0]),
                                     step=float(h), step_divisor=divisor, determinant_fd=determinant,
                                     expected_energy_ratio=predicted, relative_error=relative_error,
                                     invariant_measure_ratio=determinant / predicted))
        # A refinement difference is a diagnostic, not a guaranteed truncation bound.
        for row in measure_rows[-3:]:
            row['finest_pair_relative_difference'] = abs(determinants[-1] - determinants[-2]) / predicted
    try:
        shell_momentum([0, 0, 0], 0)
    except ValueError:
        checks.require('massless_apex_rejected', 'p=0', True)
    else:
        checks.require('massless_apex_rejected', 'p=0', False)

    rapidity_rows = []
    for xi in sorted(set([-args.max_rapidity, -3, -1, -1e-10, 0, 1e-10, 1, 3, args.max_rapidity])):
        expected = math.exp(-xi)
        measured = float((boost(xi, x) @ [1., 1., 0, 0])[0])
        relative_error = abs(measured / expected - 1)
        # Loss of significance grows as eps*exp(2*xi) for the redshifted photon.
        checks.check('photon_doppler_condition_scaled', str(xi), relative_error / max(1., math.exp(2 * xi)))
        rapidity_rows.append(dict(rapidity=xi, beta=math.tanh(xi), gamma=math.cosh(xi),
                                  gamma_beta=math.sinh(xi), expected_doppler=expected,
                                  matrix_doppler=measured, doppler_relative_error=relative_error,
                                  estimated_cancellation_floor=EPS * max(1., math.exp(2 * xi))))

    decimal_report = decimal_doppler(args.decimal_rapidity, checks)
    data = {'vectors': vector_rows, 'matrices': matrix_rows, 'rapidity': rapidity_rows,
            'composition': composition_rows, 'measure': measure_rows, 'components': component_rows}
    report = dict(passed=all(r['passed'] for r in checks.items), check_count=len(checks.items),
                  summary=checks.summary(), failed_checks=[r for r in checks.items if not r['passed']],
                  case_counts={name: len(rows) for name, rows in data.items()},
                  decimal_doppler=decimal_report,
                  exact_arithmetic=exact_report,
                  negative_controls={'wrong_passive_sign_residual': wrong_sign,
                                     'wrong_transpose_inverse_residual': wrong_inverse},
                  absolute_error_maxima={key: max(row[key] for row in matrix_rows)
                                         for key in ('metric_absolute_error', 'inverse_absolute_error', 'determinant_absolute_error')},
                  scalar_product_max_absolute_error=scalar_absolute_max,
                  measure_finest_pair_max_relative_difference=max(row['finest_pair_relative_difference'] for row in measure_rows),
                  measure_scope='Massive and nonzero massless shell; finite-difference rapidities |xi|<=3, plus noncollinear rotation/boost composition.',
                  limitations=[
                      'Finite deterministic cases and optional seeded samples do not prove group identities globally.',
                      'Scaled errors test floating-point consistency relative to product magnitudes; they do not certify relative accuracy of a small norm, determinant, or Doppler-shifted energy.',
                      'At high rapidity a tiny invariant is a difference of large terms. Inspect absolute errors and the Decimal comparison before interpreting digits.',
                      'Three finite-difference steps separate step sensitivity from matrix stress. Refinement need not be monotone after roundoff dominates and is not a rigorous error bound.',
                      'The massless cone apex p=0 is excluded; no derivative or measure density is asserted there.',
                      'The positive-shell measure also survives parity; time-reversing components do not preserve the selected positive-energy sheet.',
                      'No continuum integration, spacetime microcausality, Wigner-angle extraction, or quantum time-reversal implementation is tested.'])
    return data, report


def git_context(directory):
    """Revision metadata is qualified by dirty tracked state, never release evidence."""
    try:
        head = subprocess.run(['git', '-C', str(directory), 'rev-parse', 'HEAD'], capture_output=True, text=True, check=True).stdout.strip()
        status = subprocess.run(['git', '-C', str(directory), 'status', '--porcelain', '--untracked-files=no'], capture_output=True, text=True, check=True).stdout
        return {'head': head, 'tracked_worktree_dirty': bool(status),
                'qualification': 'HEAD alone does not identify the script or a dirty checkout; use script_sha256.'}
    except (OSError, subprocess.SubprocessError):
        return {'head': None, 'qualification': 'Git metadata unavailable; script SHA-256 remains recorded.'}


def reserve_paths(paths):
    """Refuse all existing outputs and absent-but-tracked destinations before writing."""
    if len(set(paths)) != len(paths):
        raise ValueError('output paths must be distinct')
    for path in paths:
        if path.exists() or path.is_symlink():
            raise ValueError(f'output already exists (never overwritten): {path}')
        if '.git' in [part.lower() for part in path.parts]:
            raise ValueError(f'output cannot be inside .git: {path}')
        ancestor = path.parent
        while not ancestor.exists():
            ancestor = ancestor.parent
        try:
            result = subprocess.run(['git', '-C', str(ancestor), 'rev-parse', '--show-toplevel'], capture_output=True, text=True)
            if result.returncode == 0:
                root = Path(result.stdout.strip()).resolve()
                relative = path.relative_to(root).as_posix()
                tracked = subprocess.run(['git', '-C', str(root), 'ls-files', '-z'], capture_output=True, text=True, check=True).stdout.split('\0')
                if relative.casefold() in {item.casefold() for item in tracked}:
                    raise ValueError(f'refusing tracked output destination: {path}')
        except FileNotFoundError:
            pass  # Exclusive file creation still prevents overwriting any existing artifact.


def parse_args(argv=None):
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('--output-dir', type=Path, default=Path(__file__).resolve().parent / 'lorentz-results' / 'default',
                        help='directory for six CSV files and JSON (default: lorentz-results/default beside script; choose a NEW directory for reruns)')
    parser.add_argument('--json', type=Path, help='override JSON report path; CSV files still use --output-dir')
    parser.add_argument('--samples', type=int, default=16, help='additional seeded matrix and measure cases, 0..10000 (default: 16)')
    parser.add_argument('--seed', type=int, default=20261002, help='NumPy PCG64 seed, nonnegative integer (default: 20261002)')
    parser.add_argument('--max-rapidity', type=float, default=8., help='magnitude of float64 stress boost, 0..12 (default: 8); fixed teaching cases still run')
    parser.add_argument('--decimal-rapidity', type=float, default=40., help='separate 100-digit Doppler reference, 0..60 (default: 40)')
    parser.add_argument('--fd-step', type=float, default=1e-4, help='dimensionless base h/E; also test h/2 and h/4, 1e-7..1e-2 (default: 1e-4)')
    args = parser.parse_args(argv)
    if not 0 <= args.samples <= 10000 or args.seed < 0:
        parser.error('--samples must be 0..10000 and --seed nonnegative')
    for value, low, high, name in [(args.max_rapidity, 0., 12., '--max-rapidity'),
                                   (args.decimal_rapidity, 0., 60., '--decimal-rapidity'),
                                   (args.fd_step, 1e-7, 1e-2, '--fd-step')]:
        if not math.isfinite(value) or not low <= value <= high:
            parser.error(f'{name} must be finite and in [{low}, {high}]')
    return args


def main(argv=None):
    args = parse_args(argv)
    script = Path(__file__).resolve()
    output = args.output_dir.expanduser().resolve()
    filenames = {name: output / f'lorentz-transformations-{name}.csv'
                 for name in ('vectors', 'matrices', 'rapidity', 'composition', 'measure', 'components')}
    json_path = args.json.expanduser().resolve() if args.json else output / 'lorentz-transformations-report.json'
    try:
        reserve_paths([*filenames.values(), json_path])
        data, report = run_experiment(args)
        report.update(schema_version=1, program='lorentz-transformations',
                      script_sha256=hashlib.sha256(script.read_bytes()).hexdigest(),
                      environment={'python': platform.python_version(), 'numpy': np.__version__,
                                   'platform': platform.platform(), 'float': 'IEEE-754 binary64',
                                   'rng': 'numpy.random.default_rng / PCG64'},
                      git=git_context(script.parent),
                      parameters={'samples': args.samples, 'seed': args.seed, 'max_rapidity': args.max_rapidity,
                                  'decimal_rapidity': args.decimal_rapidity, 'fd_step': args.fd_step},
                      units='c=hbar=1; momenta and masses in any one consistent energy unit',
                      conventions='eta=(+---); passive boosts with -sinh(xi) mixed entries; column vectors; rightmost matrix acts first; Rodrigues positive component rotation',
                      method='Exact Fraction fixtures; float64 Lorentz matrices and vector products; independent five-point on-shell spatial Jacobian at three steps; 100-digit Decimal exponential Doppler reference.',
                      outputs={**{name: str(path) for name, path in filenames.items()}, 'report': str(json_path)})
        for name, path in filenames.items():
            path.parent.mkdir(parents=True, exist_ok=True)
            with path.open('x', encoding='utf-8', newline='') as handle:
                writer = csv.DictWriter(handle, fieldnames=list(data[name][0]))
                writer.writeheader()
                writer.writerows(data[name])
        report['output_sha256'] = {name: hashlib.sha256(path.read_bytes()).hexdigest() for name, path in filenames.items()}
        json_path.parent.mkdir(parents=True, exist_ok=True)
        with json_path.open('x', encoding='utf-8') as handle:
            json.dump(report, handle, indent=2, allow_nan=False)
            handle.write('\n')
    except (OSError, ValueError, OverflowError, subprocess.SubprocessError) as exc:
        print(f'error: {exc}', file=sys.stderr)
        return 2
    print(f"{'PASS' if report['passed'] else 'FAIL'}: {report['check_count']} checks; Python {platform.python_version()}, NumPy {np.__version__}")
    print(f"Metric absolute max={report['absolute_error_maxima']['metric_absolute_error']:.6g}; scaled max={report['summary']['metric_product_scaled']['max_residual']:.6g}")
    print(f"Mass-shell finite-difference relative max={report['summary']['shell_measure_fd_relative']['max_residual']:.6g}; tolerance={FD_TOL:g}")
    print(f"Report: {json_path}")
    return 0 if report['passed'] else 1


if __name__ == '__main__':
    raise SystemExit(main())