"""Replay retained numerical observations. Requires Python and NumPy.
Download evidence.json beside this file, then run: python replay.py
This checks numerical consistency, not hardware provenance or integral generation.
"""
import json
from pathlib import Path
import numpy as np

def assemble(values):
    return sum(v for k,v in values.items() if k.startswith('D')) - 2*sum(v for k,v in values.items() if k.startswith('M'))

def main():
    data=json.loads(Path(__file__).with_name('evidence.json').read_text())
    references={k:float(np.linalg.eigvalsh(h)[0]) for k,h in data['hamiltonians'].items()}
    assert abs(assemble(references)-data['reference_energy_hartree']) < 1e-12
    for observation, expected in zip(data['observations'],data['results']):
        energies={}; shots=0
        for phase in observation['phases'].values():
            samples=phase['samples']; shots+=sum(s['shots'] for s in samples)
            size=len(data['hamiltonians'][next(iter(phase['fragments']))]); n=int(np.log2(size))
            counts=np.array([sum(s['counts'].get(format(i,f'0{n}b'),0) for s in samples) for i in range(size)],dtype=float)
            assert int(counts.sum())==sum(s['shots'] for s in samples)
            p=counts/counts.sum()
            features=np.array([[1]+[2*int(b)-1 for b in format(i,f'0{n}b')] for i in range(size)])
            for fragment,fit in phase['fragments'].items():
                c=np.tanh(features@np.array(fit['weights']))*np.sqrt(p)
                h=np.array(data['hamiltonians'][fragment]); energy=float(c@h@c/(c@c))
                assert abs(energy-fit['energy'])<1e-12
                assert abs(energy-references[fragment])<=data['fragment_tolerance_hartree']
                energies[fragment]=energy
        if 'late_tail' in observation:
            shots+=sum(s['shots'] for s in observation['late_tail']['samples'])
        total=assemble(energies); error=abs(total-data['reference_energy_hartree'])
        assert shots==expected['shots']
        assert abs(total-expected['energy_hartree'])<1e-12
        assert abs(error-expected['absolute_error_hartree'])<1e-12
        print(f"{observation['label']}: {shots:,} shots; energy {total:.12f} Ha; error {error:.12f} Ha; PASS")

if __name__=='__main__':
    main()