"""Read the original failed realization; cross-check its exact raw random inputs."""
from pathlib import Path
import argparse
import hashlib
import json
import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.stats import kstest,chi2_contingency


def audit(directory):
    data=np.load(directory/'original-projection-inputs.npz')
    radial=data['radial_uniforms'];cosine=data['direction_cosines'];azimuth=data['azimuths'];positions=data['positions']
    uniforms=np.column_stack([radial,(cosine+1)/2,azimuth/(2*np.pi)])
    state=json.loads((directory/'original-projection-rng-state.json').read_text())
    generator=np.random.Generator(np.random.PCG64());generator.bit_generator.state=state
    regenerated=[generator.uniform(1e-10,1-1e-10,len(radial)),generator.uniform(-1,1,len(radial)),generator.uniform(0,2*np.pi,len(radial))]
    replay=[bool(np.array_equal(a,b)) for a,b in zip(regenerated,[radial,cosine,azimuth])]
    assert all(replay)
    # Independently solve enclosed number probability by adaptive radial integration.
    chosen=np.linspace(0,len(radial)-1,128,dtype=int)
    independent_radii=[]
    for uniform in radial[chosen]:
        cumulative=lambda radius:quad(lambda x:3*x*x/(1+x*x)**2.5,0,radius,epsabs=1e-12)[0]
        high=1.
        while cumulative(high)<uniform:high*=2
        independent_radii.append(brentq(lambda r:cumulative(r)-uniform,0,high,xtol=1e-12))
    independent_radii=np.array(independent_radii)
    # Polar-angle construction uses the same raw inputs with no shared sampler call.
    theta=np.arccos(cosine[chosen])
    independent_positions=independent_radii[:,None]*np.column_stack([np.sin(theta)*np.cos(azimuth[chosen]),np.sin(theta)*np.sin(azimuth[chosen]),np.cos(theta)])
    position_error=float(np.max(np.abs(independent_positions-positions[chosen])))
    assert position_error<1e-10
    radius=np.linalg.norm(positions,axis=1)
    projected=np.hypot(positions[:,0],positions[:,1])
    same_input_algebra=np.sqrt(radius**2*(1-cosine**2))
    projected_error=float(np.max(np.abs(projected-same_input_algebra)));assert projected_error<1e-10
    uniform_tests=[]
    for index,name in enumerate(['radial','direction-cosine','azimuth']):
        result=kstest(uniforms[:,index],'uniform')
        uniform_tests.append({'input':name,'n':len(radial),'ksD':float(result.statistic),'pValueDiagnostic':float(result.pvalue)})
    contingency=[]
    for a,b in [(0,1),(0,2),(1,2)]:
        bins=np.histogram2d(uniforms[:,a],uniforms[:,b],bins=16,range=[[0,1],[0,1]])[0]
        statistic,pvalue,dof,_=chi2_contingency(bins)
        contingency.append({'columns':[a,b],'binsPerAxis':16,'chiSquared':float(statistic),'dof':int(dof),'pValueDiagnostic':float(pvalue)})
    conditional=[]
    for cut in [.1,.5,1,2,5]:
        probability=np.ones(len(radius));outside=radius>cut
        probability[outside]=1-np.sqrt(1-(cut/radius[outside])**2)
        conditional.append({'radius':cut,'actualNumberCDF':float(np.mean(projected<cut)),
                            'expectedForTheseRadiiWithIsotropicAngles':float(np.mean(probability)),
                            'unconditionalNumberCDF':cut*cut/(1+cut*cut)})
    result=kstest(projected,lambda r:r*r/(1+r*r))
    return {'status':'Original failed sample preserved; same-input independent geometry agrees; statistical anomaly remains explicit',
      'source':str(directory/'original-projection-inputs.npz'),'originalSeed':20261104,'rngStateReplayExact':replay,
      'inputArraySha256':{name:hashlib.sha256(data[name].tobytes()).hexdigest() for name in data.files},
      'rawInputUniformity':uniform_tests,'rawInputPearsonCorrelation':np.corrcoef(uniforms,rowvar=False).tolist(),
      'rawInputPairwiseContingency':contingency,'independentAdaptiveRadialRoots':len(chosen),
      'independentSameInputPositionError':position_error,'projectedRadiusIndependentError':projected_error,
      'sameSampleProjectedKSD':float(result.statistic),'sameSampleProjectedKSPValueDiagnostic':float(result.pvalue),
      'conditionalOnActualRadii':conditional,'numberCDFNotMassWeighted':True,
      'weights':'Every saved point has weight1/n. No stellar masses enter this projection validation; imposed mass segregation is a different model.',
      'remaining':'Uniformity/independence tests are diagnostics, not proofs. Do not replace this failed sample. An independently predeclared seed cohort is a separate experiment, not a retroactive pass.'}


if __name__=='__main__':
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--inputs',type=Path,default=Path('results'))
    parser.add_argument('--report',type=Path,default=Path('results/projection-audit.json'))
    args=parser.parse_args();result=audit(args.inputs);args.report.parent.mkdir(parents=True,exist_ok=True);args.report.write_text(json.dumps(result,indent=2)+'\n')
    print(json.dumps({k:result[k] for k in ['status','rngStateReplayExact','rawInputUniformity','rawInputPearsonCorrelation','independentSameInputPositionError','sameSampleProjectedKSD','sameSampleProjectedKSPValueDiagnostic']},indent=2))