#!/usr/bin/env python3
"""Finite probes for the Green compact inverse, using only Python's standard library.

No finite run proves infinite-dimensional convergence. The displayed error ranges
combine a finite coefficient sum with the analytic p-series tail bound in the
companion article. Decimal endpoints use binary64 arithmetic, not outward-rounded
interval arithmetic; they are numerical evaluations of analytic bounds.
"""
from __future__ import annotations
import argparse
import json
import math
import sys

MAX_TERMS = 200_000
PI2 = math.pi**2
PI4 = math.pi**4


def validate(delta: float, mode: int, terms: int) -> None:
    if not math.isfinite(delta) or not 0.0 < delta <= 1.0:
        raise ValueError("delta must be finite and satisfy 0 < delta <= 1 in this normalized example")
    if isinstance(mode, bool) or not isinstance(mode, int) or not 1 <= mode <= MAX_TERMS:
        raise ValueError(f"mode must be an integer between 1 and {MAX_TERMS}")
    if isinstance(terms, bool) or not isinstance(terms, int) or not 1 <= terms <= MAX_TERMS:
        raise ValueError(f"terms must be an integer between 1 and {MAX_TERMS}")
    if terms < mode:
        raise ValueError("terms must be at least mode, so the noisy coefficient is included exactly")


def sigma(n: int) -> float:
    return 1.0 / (PI2 * n * n)


def run_probe(delta: float = 1e-6, mode: int = 100, terms: int = 2000) -> dict:
    validate(delta, mode, terms)
    tau = math.sqrt(delta)
    alpha = delta
    cutoff = math.floor(1.0 / (math.pi * math.sqrt(tau)))
    if cutoff > 2**50:
        raise ValueError("delta is too small to resolve adjacent retained-mode indices in binary64")
    # Correct any binary64 floor discrepancy against the actual threshold test.
    while cutoff > 0 and sigma(cutoff) < tau:
        cutoff -= 1
    while sigma(cutoff + 1) >= tau:
        cutoff += 1
    coefficients = []
    tsvd_squares = []
    tikh_squares = []
    for n in range(1, terms + 1):
        s = sigma(n)
        truth = PI2 / n
        noise_data = delta if n == mode else 0.0
        direct_noise = noise_data / s
        kept = n <= cutoff
        tsvd = truth + direct_noise if kept else 0.0
        # Stable expression for the error avoids subtracting two nearby truths.
        ratio = alpha / (s*s + alpha)
        tikh_error = -ratio * truth + s * noise_data / (s*s + alpha)
        tsvd_error = direct_noise if kept else -truth
        tsvd_squares.append(tsvd_error**2)
        tikh_squares.append(tikh_error**2)
        if n <= 6 or n == mode:
            coefficients.append({"n": n, "sigma": s, "truth": truth,
                                 "tsvd": tsvd, "tikhonov": truth + tikh_error})
    tsvd_head = math.fsum(tsvd_squares)
    tikh_head = math.fsum(tikh_squares)
    # All noise lies among these finite coefficients because terms >= mode.
    # TSVD tail begins beyond both the finite sum and the cutoff.
    tail_start = max(terms, cutoff)
    tsvd_range = [math.sqrt(tsvd_head + PI4/(tail_start + 1)),
                  math.sqrt(tsvd_head + PI4/tail_start)]
    tikh_range = [math.sqrt(tikh_head), math.sqrt(tikh_head + PI4/terms)]
    # A separate source-compatible truth: w=e1, rho=1, x*=T*T w=pi^-4 e1.
    source_alpha = delta**(2.0/3.0)
    s1 = sigma(1)
    first_bias = -source_alpha * s1*s1 / (s1*s1 + source_alpha)
    sm = sigma(mode)
    source_noise = sm * delta / (sm*sm + source_alpha)
    source_actual = abs(first_bias + source_noise) if mode == 1 else math.hypot(first_bias, source_noise)
    result = {
        "scope": "finite coefficients plus analytic infinite-tail bounds; no simulation proves convergence",
        "arithmetic": "binary64 evaluations, not outward-rounded interval certificates",
        "input": {"delta": delta, "noise_mode": mode, "terms": terms, "truth_data": "g_3"},
        "parameters": {"tau": tau, "alpha": alpha, "retained_modes": cutoff,
                       "cutoff_arithmetic": "adjacent indices checked with binary64 sigma >= tau; exact equality follows the article, not rounded floats"},
        "direct_inverse_error": PI2 * mode * mode * delta,
        "tsvd_error_range_from_tail_bounds": tsvd_range,
        "tikhonov_error_range_from_tail_bounds": tikh_range,
        "uniform_noise_bounds": {"tsvd": delta/tau, "tikhonov": delta/(2*math.sqrt(alpha))},
        "selected_solution_coefficients": coefficients,
        "separate_source_example": {
            "truth": "pi^-4 e_1; w=e_1 and rho=1", "alpha": source_alpha,
            "actual_error_for_same_noise": source_actual,
            "bound_alpha_plus_noise": source_alpha + delta/(2*math.sqrt(source_alpha)),
            "bound_three_halves_delta_two_thirds": 1.5*delta**(2.0/3.0),
        },
    }
    if not all(math.isfinite(x) for x in [*tsvd_range, *tikh_range, source_actual]):
        raise ArithmeticError("a floating-point output is non-finite; choose less extreme parameters")
    return result


def self_test() -> dict:
    count = 0
    def check(condition: bool, message: str) -> None:
        nonlocal count
        if not condition:
            raise RuntimeError("self-test failed: " + message)
        count += 1
    for delta, mode, terms in [(1e-6, 100, 2000), (1e-3, 1, 20), (1.0, 3, 3), (1e-12, 2, 2)]:
        r = run_probe(delta, mode, terms)
        lo, hi = r["tsvd_error_range_from_tail_bounds"]
        check(0 <= lo <= hi, "TSVD tail range is ordered")
        lo, hi = r["tikhonov_error_range_from_tail_bounds"]
        check(0 <= lo <= hi, "Tikhonov tail range is ordered")
        s = r["separate_source_example"]
        check(s["actual_error_for_same_noise"] <= s["bound_alpha_plus_noise"]*(1+1e-12), "source-condition bound")
        check(math.isclose(s["bound_alpha_plus_noise"], s["bound_three_halves_delta_two_thirds"], rel_tol=1e-12), "three-halves constant")
    for delta, mode, terms in [(0.,1,1), (-1.,1,1), (math.nan,1,1), (math.inf,1,1), (2.,1,1), (1e-6,0,1), (1e-6,2,1), (1e-6,1,0), (1e-6,1,MAX_TERMS+1), (1e-200,1,1)]:
        try:
            run_probe(delta, mode, terms)
        except ValueError:
            count += 1
        else:
            raise RuntimeError("self-test failed: invalid input was accepted")
    for n in range(1, 9):
        r = run_probe(sigma(n)**2, 1, 10)
        tau = r["parameters"]["tau"]
        k = r["parameters"]["retained_modes"]
        check(k == sum(sigma(j) >= tau for j in range(1, 20)), "threshold-adjacent binary64 comparison")
    for n in range(1, 5):
        omega = (n-0.5)*math.pi
        check(abs(math.cos(omega)) < 1e-14, "Volterra mixed boundary")
        # Tests finite antiderivative identities at fixed points, not completeness.
        for t in (0., 0.25, 1.):
            left = math.sqrt(2)*(math.cos(omega*t)-math.cos(omega))/omega
            right = math.sqrt(2)*math.cos(omega*t)/omega
            check(math.isclose(left, right, abs_tol=1e-14), "Volterra adjoint integral")
    return {"finite_checks_passed": count, "scope": "finite arithmetic and input errors only; analytic proofs remain in the article"}


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument("--delta", type=float, default=1e-6, help="normalized noise norm, 0 < delta <= 1")
    parser.add_argument("--mode", type=int, default=100, help="positive Green noise mode")
    parser.add_argument("--terms", type=int, default=2000, help="finite sum length, must be >= mode")
    parser.add_argument("--self-test", action="store_true", help="run finite checks and invalid-input checks")
    args = parser.parse_args()
    try:
        result = self_test() if args.self_test else run_probe(args.delta, args.mode, args.terms)
    except (ValueError, ArithmeticError, RuntimeError) as error:
        print("error: " + str(error), file=sys.stderr)
        return 2
    print(json.dumps(result, ensure_ascii=False, indent=2, allow_nan=False))
    return 0


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