#!/usr/bin/env python3
"""Deterministic, RMS-matched diagnostic test; requires Python 3.9+ and NumPy.

Run: python3 generate_matched_test.py [--output-dir DIRECTORY]
Writes matched-test.json only. No network or other project files required.

Protocol fixed before these calculations: beta=0, active Majoranas 1..7,
lambda=0,.1,.2,.3; inclusive time grid 12..24 in steps .1. D is the arithmetic
mean of |G_i(t)| across those operators and grid points. Directional hypothesis:
D(lambda)<D(0) for EVERY tested nonzero lambda. No window/operator/lambda tuning.
This newly stated hypothesis is not preregistered, blind, or a novelty claim;
earlier unscaled calculations were inspected. Refinement to .05 is a numerical
sensitivity check, not a statistical confidence interval. The primary estimator
is a discrete mean, not trapezoidal quadrature; we also report normalized
trapezoidal integrals at both spacings as a separate time-integration check.
"""
import argparse
import hashlib
import json
import math
from functools import reduce
from pathlib import Path
import numpy as np

DIM=16
TOL=1e-11
LAMBDAS=(0., .1, .2, .3)
TERMS=((-0.36,(1,2,4,5)),(0.19,(1,3,4,7)),(-0.71,(1,3,5,6)),
       (0.22,(2,3,4,6)),(0.49,(2,3,5,7)))
PERTURBATION=(1,2,3,5)


def maxabs(a):
    return float(np.max(np.abs(a)))


def checked(error, tolerance=TOL):
    error=float(error)
    if not math.isfinite(error) or error>tolerance:
        raise AssertionError(f'Numerical check failed: {error} > {tolerance}')
    return {'passed':True,'max_absolute_error':error,'tolerance':tolerance}


def matrices():
    I=np.eye(2,dtype=complex)
    X=np.array([[0,1],[1,0]],dtype=complex)
    Y=np.array([[0,-1j],[1j,0]],dtype=complex)
    Z=np.diag([1,-1]).astype(complex)
    psi=[reduce(np.kron,[Z]*j+[P]+[I]*(3-j))/math.sqrt(2)
         for j in range(4) for P in (X,Y)]
    term=lambda indices:reduce(np.matmul,[psi[i-1] for i in indices])
    products=[term(idx) for _,idx in TERMS]
    return psi,products,sum(c*M for (c,_),M in zip(TERMS,products)),term(PERTURBATION)


def center_and_rms(H):
    centered=H-np.trace(H)/DIM*np.eye(DIM)
    return centered,float(np.sqrt(np.trace(centered@centered).real/DIM))


def independent_exp(A):
    norm=float(np.linalg.norm(A,ord=np.inf))
    power=max(0,math.ceil(math.log2(norm))) if norm else 0
    B=A/2**power
    out=np.eye(DIM,dtype=complex)
    term=out.copy()
    for k in range(1,150):
        term=term@B/k
        out+=term
        if np.linalg.norm(term,ord=np.inf)<1e-17:
            break
    else:
        raise AssertionError('Exponential series did not converge')
    for _ in range(power):
        out=out@out
    return out


def correlations(H,psi,times):
    E,V=np.linalg.eigh(H)
    gaps=E[:,None]-E[None,:]
    weights=[np.abs(V.conj().T@p@V)**2 for p in psi]
    return np.array([[2/DIM*np.sum(w*np.cos(gaps*t)) for t in times]
                     for w in weights]),E,V


def statistics(curves,step):
    per_time=np.mean(np.abs(curves[:7]),axis=0)
    # Uniform-grid trapezoid integral / window duration, implemented without SciPy.
    integral=step*(.5*per_time[0]+np.sum(per_time[1:-1])+.5*per_time[-1])/12
    return {'D':float(np.mean(per_time)),
            'per_operator_mean_absolute':np.mean(np.abs(curves[:7]),axis=1).tolist(),
            'normalized_trapezoidal_integral':float(integral)}


def generate():
    psi,terms,H0,P=matrices()
    H0_centered,target_rms=center_and_rms(H0)
    times=np.arange(120,241)/10
    refined_times=np.arange(240,481)/20
    identity=np.eye(DIM)
    validation={
        'majorana_algebra':checked(max(maxabs(a@b+b@a-(identity if i==j else 0))
             for i,a in enumerate(psi) for j,b in enumerate(psi))),
        'baseline_term_commutation':checked(max(maxabs(a@b-b@a) for a in terms for b in terms)),
        'baseline_centering':checked(maxabs(H0_centered-H0)),
    }
    zero,_,_=correlations(np.zeros((DIM,DIM)),psi,times)
    validation['zero_hamiltonian']=checked(maxabs(zero-1))
    cases=[]
    for lam in LAMBDAS:
        raw=H0+lam*P
        centered,raw_rms=center_and_rms(raw)
        factor=target_rms/raw_rms
        H=centered*factor
        _,matched_rms=center_and_rms(H)
        curves,E,V=correlations(H,psi,times)
        refined,_,_=correlations(H,psi,refined_times)
        initial,_,_=correlations(H,psi,[0.])
        direct_error=0.
        unitary_error=0.
        for k,t in enumerate(times):
            U=(V*np.exp(-1j*E*t))@V.conj().T
            values=np.array([2/DIM*np.trace(U.conj().T@p@U@p).real for p in psi])
            direct_error=max(direct_error,maxabs(values-curves[:,k]))
            unitary_error=max(unitary_error,maxabs(U.conj().T@U-identity))
        independent_error=0.
        independent_U_error=0.
        for t in (0.,12.,15.,18.,21.,24.):
            U=independent_exp(-1j*H*t)
            spectral_U=(V*np.exp(-1j*E*t))@V.conj().T
            independent_U_error=max(independent_U_error,maxabs(U-spectral_U))
            values=np.array([2/DIM*np.trace(U.conj().T@p@U@p).real for p in psi])
            predicted,_,_=correlations(H,psi,[t])
            independent_error=max(independent_error,maxabs(values-predicted[:,0]))
        coarse_stats=statistics(curves,.1)
        fine_stats=statistics(refined,.05)
        # Independent centered RMS via variance of eigenvalues.
        eigen_rms=float(np.sqrt(np.mean((E-np.mean(E))**2)))
        checks={
            'hermiticity':checked(maxabs(H-H.conj().T)),
            'centered_trace':checked(abs(np.trace(H)/DIM)),
            'matched_RMS_matrix':checked(abs(matched_rms-target_rms)),
            'matched_RMS_eigenvalues':checked(abs(eigen_rms-target_rms)),
            'initial_all_one':checked(maxabs(initial-1)),
            'spectator_commutes':checked(maxabs(H@psi[7]-psi[7]@H)),
            'spectator_conserved':checked(maxabs(curves[7]-1)),
            'direct_U_vs_spectral_all_times':checked(direct_error),
            'unitarity_all_times':checked(unitary_error),
            'independent_exponential_U':checked(independent_U_error),
            'independent_exponential_correlations':checked(independent_error),
            'nested_grid_agreement':checked(maxabs(curves-refined[:,::2])),
        }
        cases.append({
            'perturbation_lambda':lam,
            'raw_centered_RMS_energy':raw_rms,
            'rescaling_factor':factor,
            'matched_centered_RMS_energy':matched_rms,
            'matched_spectral_width':float(E[-1]-E[0]),
            'traces':np.round(curves,14).tolist(),
            'mean_absolute_active7':np.round(np.mean(np.abs(curves[:7]),axis=0),14).tolist(),
            'primary':coarse_stats,
            'refinement_step_0_05':fine_stats,
            'refinement_changes':{
                'D_refined_minus_primary':fine_stats['D']-coarse_stats['D'],
                'trapezoid_refined_minus_coarse':fine_stats['normalized_trapezoidal_integral']-coarse_stats['normalized_trapezoidal_integral'],
            },
            'validation':checks,
        })
    baseline=cases[0]
    for case in cases:
        contrast=case['primary']['D']-baseline['primary']['D']
        refined_contrast=case['refinement_step_0_05']['D']-baseline['refinement_step_0_05']['D']
        case['contrast_D_minus_baseline']=contrast
        case['refined_contrast_D_minus_baseline']=refined_contrast
        case['contrast_refinement_change']=refined_contrast-contrast
        case['hypothesis_comparison']='baseline' if case['perturbation_lambda']==0 else ('supports_direction' if contrast<0 else 'fails_direction')
        case['refined_direction_matches_primary']=(contrast<0)==(refined_contrast<0)
    failures=[c['perturbation_lambda'] for c in cases[1:] if c['contrast_D_minus_baseline']>=0]
    supports=[c['perturbation_lambda'] for c in cases[1:] if c['contrast_D_minus_baseline']<0]
    conclusion='supported_on_tested_grid' if not failures else 'fails_for_at_least_one_tested_lambda'
    return {
        'schema_version':1,
        'metadata':{
            'title':'Fixed-window mean absolute autocorrelation with matched centered RMS energy',
            'state':'beta=0; rho=I16/16',
            'beta':0,
            'normalization':'psi_i=gamma_i/sqrt(2); {psi_i,psi_j}=delta_ij I; hbar=1',
            'majorana_spin_mapping':'gamma_(2j-1)=Z^(tensor j-1) tensor X tensor I^(tensor 4-j); gamma_(2j) replaces X with Y; j=1..4, first qubit leftmost',
            'hamiltonian_terms':[{'coefficient':c,'indices':list(idx)} for c,idx in TERMS],
            'perturbation_indices':list(PERTURBATION),
            'matching':'H_matched(lambda)=[H_raw(lambda)-Tr(H_raw(lambda))/16 I] * RMS(H0)/RMS(H_raw(lambda)); RMS(H)=sqrt(Tr[(H-Tr(H)/16 I)^2]/16)',
            'target_centered_RMS_energy':target_rms,
            'time_units':'Time remains in inverse baseline coefficient-energy units, hbar=1; H0 is unchanged. No physical seconds assigned.',
            'observable':'G_i(t)=(2/16) Re Tr[exp(+iHt) psi_i exp(-iHt) psi_i]',
            'spectator':'psi8 commutes with both models and G8=1; it is plotted in supplied traces but excluded from D.',
            'array_layout':'traces[operator_index_zero_based][time_index]; mean_absolute_active7[time_index]',
            'numpy_version':np.__version__,
            'sources':[
                {'url':'https://arxiv.org/abs/2302.07897','use':'Published learned H0 coefficients and normalization; critical analysis'},
                {'url':'https://arxiv.org/abs/2303.15423','use':'Authors response, including the quartic perturbation with coefficient 0.3'},
            ],
        },
        'fixed_protocol':{
            'declared_before_this_calculation':True,
            'not_preregistered':True,
            'prior_exposure':'Earlier unscaled lambda=0 and .3 calculations were inspected. This is a newly stated diagnostic hypothesis, not blind validation or a novelty claim.',
            'active_operator_indices':[1,2,3,4,5,6,7],
            'lambdas':list(LAMBDAS),
            'window':[12.,24.],
            'step':.1,
            'endpoints':'Both included, 121 samples',
            'D_definition':'Arithmetic mean across all 121 selected times and seven active operators of abs(G_i(t)); no postselection',
            'contrast':'D(lambda)-D(0)',
            'directional_hypothesis':'D(lambda)<D(0) for every tested nonzero lambda',
            'fixed_choices':'No tuning of window, operators, lambdas or matching rule after obtaining these results',
            'refinement':'Same endpoints and operators, step .05 (241 samples). Refinement changes are numerical discretization sensitivity, not statistical confidence intervals.',
            'quadrature_note':'Primary D is a finite-grid arithmetic mean. Separately supplied normalized trapezoidal integrals approximate a continuous-time average and are not substituted for primary D.',
        },
        'times':times.tolist(),
        'cases':cases,
        'validation':validation,
        'result':{
            'hypothesis_status':conclusion,
            'direction_supported_lambdas':supports,
            'direction_failed_lambdas':failures,
            'all_refined_directions_match':all(c['refined_direction_matches_primary'] for c in cases),
            'interpretation': 'The fixed-grid directional hypothesis '+('holds for all tested nonzero lambdas.' if not failures else 'does not hold for all tested nonzero lambdas.')+' These deterministic finite-system results quantify one chosen window and observable. Refinement is a numerical sensitivity check, not statistical uncertainty.',
            'limitations':[
                'Matching centered RMS removes this particular overall energy-scale difference, not every confound; spectral width and spectrum shape may still differ.',
                'Adding a quartic term also changes detailed couplings, conserved quantities and operator dynamics. This is not causal proof isolating noncommutation.',
                'No conclusion about gravitation, emergent geometry, thermalization, teleportation or experimental replication follows.',
                'The tested window and operators were fixed for this calculation, but prior unscaled data were inspected. No preregistration, blind test or novel-physics claim is made.',
                'No ensemble sampling or statistical confidence intervals are used. Numerical refinement changes are not rigorous bounds on continuum error.',
            ],
        },
    }


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output-dir',type=Path,default=Path(__file__).resolve().parent)
    args=parser.parse_args()
    result=generate()
    result['metadata']['generator_sha256']=hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
    args.output_dir.mkdir(parents=True,exist_ok=True)
    path=args.output_dir/'matched-test.json'
    path.write_text(json.dumps(result,indent=2,allow_nan=False)+'\n',encoding='utf-8')
    print(f'Wrote {path}')
    print('lambda | RMS factor | D(.1) | contrast | D(.05) | refined contrast')
    for c in result['cases']:
        print(f"{c['perturbation_lambda']:.1f} | {c['rescaling_factor']:.12f} | {c['primary']['D']:.12f} | {c['contrast_D_minus_baseline']:+.12f} | {c['refinement_step_0_05']['D']:.12f} | {c['refined_contrast_D_minus_baseline']:+.12f}")
    print(json.dumps(result['result'],indent=2))


if __name__=='__main__':
    main()