"""Independent finite limits, proper likelihood, actual fit/source and failure controls."""
import argparse,copy,hashlib,json,math,signal,time
from pathlib import Path
import numpy as np
from scipy.linalg import cho_factor,cho_solve
from scipy.optimize import brentq
from analyze import load_saved
from conditional_fit import fit
from verify_analytic_limits import verify as limits
from verify_nuisance_mapping import verify as geometry
from rotation_model import forward,log_likelihood
from physical_units import GRAVITY_KPC_KMS2_PER_SOLAR_MASS as GRAVITY

def clean(value):
 if isinstance(value,dict):return {key:clean(item) for key,item in value.items()}
 if isinstance(value,(list,tuple)):return [clean(item) for item in value]
 if isinstance(value,np.ndarray):return clean(value.tolist())
 if isinstance(value,np.generic):return value.item()
 return value

def main():
 parser=argparse.ArgumentParser();parser.add_argument('--input',type=Path,required=True);parser.add_argument('--analysis',type=Path,required=True);parser.add_argument('--output',type=Path,required=True);args=parser.parse_args();started=time.monotonic()
 def cap(signum,frame):raise TimeoutError('120s numerical cap')
 signal.signal(signal.SIGALRM,cap);signal.setitimer(signal.ITIMER_REAL,120)
 record=load_saved(args.input);analysis=json.loads((args.analysis/'analysis.json').read_text());assert analysis['inputSha256']==hashlib.sha256(args.input.read_bytes()).hexdigest()
 independent_limits=limits();mapping=geometry(record)
 saved_checks=[]
 for kind,solution in analysis['fits'].items():
  prediction=forward(record,solution['massToLight'],solution['haloDensitySolarPerKpc3'],solution['haloScaleKpc'],kind,solution['distanceRatio'],solution['inclinationDegrees'])
  max_difference=max(float(np.max(np.abs(prediction[key]-np.asarray(value)))) for key,value in solution['predictedRecord'].items())
  likelihood_difference=abs(log_likelihood(record,prediction)-solution['logWorkingGaussianLikelihood'])
  saved_checks.append({'kind':kind,'wholeSavedPredictionMaximumAbsoluteDifference':max_difference,'sameSavedLikelihoodAbsoluteDifference':likelihood_difference,'sixActualAttempts':len(solution['attempts']),'boundaryBest':solution['boundaryBest'],'passed':max_difference<=1e-8 and likelihood_difference<=1e-8 and len(solution['attempts'])==6})
 poisoned=copy.deepcopy(record);poisoned['injectedTruth']={'massToLight':float('nan'),'density':-1,'scale':-99,'fakeOptimum':'must never be read'}
 poisoned_fit=fit(poisoned,'nfw');truth_difference=abs(poisoned_fit['logWorkingGaussianLikelihood']-analysis['fits']['nfw']['logWorkingGaussianLikelihood'])
 count=9;sigma=2.5;common=3.;covariance=sigma*sigma*np.eye(count)+common*common*np.ones((count,count));residual=np.linspace(-3,5,count)
 numeric=float(residual@cho_solve(cho_factor(covariance,lower=True),residual));analytic=float(residual@residual/sigma**2-common**2*residual.sum()**2/(sigma**2*(sigma**2+count*common**2)))
 covariance_error=abs(numeric-analytic);mean_variance=float(covariance.sum()/count**2)
 feedback_checks=[]
 for removed_fraction in [.05,.1,.2]:
  energy=removed_fraction-.5
  turning_radius=brentq(lambda radius:.5/radius**2-(1-removed_fraction)/radius-energy,1.0001,5.,xtol=1e-13)
  analytic_radius=1/(1-2*removed_fraction);error=abs(turning_radius-analytic_radius)
  feedback_checks.append({'removedMassFraction':removed_fraction,'rootOfTurningEnergyEquation':turning_radius,'analyticApocentre':analytic_radius,'absoluteDifference':error,'gate':1e-12,'passed':error<=1e-12})
 toy={'flatSpeed200At10KpcHaloSphereMassSolar':200**2*10/GRAVITY,'flatSpeed200At20KpcHaloSphereMassSolar':200**2*20/GRAVITY,
      'signedGasToySpeedKms':math.sqrt(-20**2+.5*80**2+100**2),'incorrectSquaredGasToySpeedKms':math.sqrt(20**2+.5*80**2+100**2),
      'distance3p16To3p32Ratio':3.32/3.16,'inclination63To60VelocityRatio':math.sin(math.radians(63))/math.sin(math.radians(60)),
      'wrongScaledPerfectFitLogPreference72Over54For48':48*math.log(math.sin(math.radians(72))/math.sin(math.radians(54))),
      'feedbackInstantPointMassDropFraction':.1,'feedbackPointMassApocentreOverInitialRadius':1/(1-2*.1),
      'deepMONDIllustrativeMassSolar':1e10,'deepMONDSpeedKms':(1.3271244e20*1e10*1.2e-10)**.25/1000}
 result={'status':'PASS' if independent_limits['allRequiredPassed'] and mapping['allRequiredPassed'] and all(item['passed'] for item in saved_checks) and truth_difference<=1e-8 and covariance_error<=1e-10 and all(item['passed'] for item in feedback_checks) else 'FAIL',
 'finiteAnalyticLimits':independent_limits,'distanceInclinationAndLikelihood':mapping,'savedRecordChecks':saved_checks,'poisonedUnusedTruthLogLikelihoodDifference':truth_difference,
 'independentCorrelatedErrorExample':{'points':count,'independentSigmaKms':sigma,'commonSigmaKms':common,'meanVarianceKmsSquared':mean_variance,'incorrectlyIndependentMeanVarianceKmsSquared':(sigma*sigma+common*common)/count,'choleskyQuadratic':numeric,'rankOneAnalyticQuadratic':analytic,'absoluteDifference':covariance_error,'gate':1e-10},
 'independentInstantMassDropTurningRadii':feedback_checks,'workedNumbers':toy,'wallSeconds':time.monotonic()-started,'scope':'Finite mathematical/data-map/refit integrity checks, not calibrated statistical coverage, global optimum or theory selection'}
 args.output.parent.mkdir(parents=True,exist_ok=True)
 with args.output.open('x') as stream:stream.write(json.dumps(clean(result),indent=2,allow_nan=False)+'\n')
 signal.setitimer(signal.ITIMER_REAL,0);print(json.dumps({'status':result['status'],'output':str(args.output),'wallSeconds':result['wallSeconds'],'workedNumbers':toy},indent=2))
 if result['status']!='PASS':raise SystemExit(2)
if __name__=='__main__':main()