"""Independent limits and actual coordinate/source/mixture refinement. No preflightconstants as producers."""
import signal
signal.signal(signal.SIGALRM,lambda *_:(_ for _ in()).throw(TimeoutError('120s producer cap')));signal.alarm(120)
import argparse,math,hashlib
import numpy as np
from scipy.integrate import quad,solve_ivp,dblquad
from evolution import histories,formed,incoming,delay_cdf,response,simulate
from compare import finite_rows,log_probabilities
from io_helpers import PACKAGE,read_json,write_json,family,fixed_inputs
from pathlib import Path

def main():
 parser=argparse.ArgumentParser();parser.add_argument('--records',type=Path,default=Path('results/records'));parser.add_argument('--reference',type=Path,default=Path('results/reference'));parser.add_argument('--real-data',type=Path);parser.add_argument('--output',type=Path,default=Path('results/validation.json'));args=parser.parse_args();fixed=fixed_inputs();maximum={'bulk':0.,'elementScaled':0.,'abundanceDex':0.,'gasAnalytic':0.,'gasInventory':0.,'elementInventory':0.,'mixtureMeanLogDensity':0.,'convolutionRelative':0.,'sourceIntegral':0.};details=[]
 datasets=[('synthetic',finite_rows(read_json(PACKAGE/'inputs/own-observations.json')))];
 if args.real_data:datasets.append(('Hillconditional',finite_rows(read_json(args.real_data))))
 initialMg=fixed['initialFeFraction']*fixed['netCC']['mg']/fixed['netCC']['fe']
 for p in family():
  path=args.records/(p['name']+'.npz');coarse=read_json(path);reference=read_json(args.reference/(p['name']+'.npz'));a=np.array(coarse['states']);time=np.array(coarse['timeGyr']);bracket=np.array(coarse['abundances'])
  assert reference['coarseRecordSha256']==hashlib.sha256(path.read_bytes()).hexdigest()
  for key in ['bulk','elementScaled','abundanceDex']:maximum[key]=max(maximum[key],reference['refinement'][key])
  for row in reference['observationRefinement']:
   data_path=PACKAGE/'inputs/own-observations.json'if row['data']=='synthetic'else args.real_data
   if data_path is None:continue
   assert row['dataSha256']==hashlib.sha256(data_path.read_bytes()).hexdigest()
   maximum['mixtureMeanLogDensity']=max(maximum['mixtureMeanLogDensity'],row['adjacentBirthBinMeanAbsoluteLogDensity'],row['productionToPositiveGL8MeanAbsoluteLogDensity'],row['GL4to8MeanAbsoluteLogDensity'])
  assert reference['earlyPositiveSupport']['firstFeFraction']>0 and reference['earlyPositiveSupport']['firstMgFraction']>0
  exact=np.array([1+incoming(t,p)-(.6+p['eta'])*formed(t,p)for t in time]);maximum['gasAnalytic']=max(maximum['gasAnalytic'],float(np.max(np.abs(a[:,0]-exact))));maximum['gasInventory']=max(maximum['gasInventory'],float(np.max(np.abs(a[:,0]+a[:,1]+a[:,2]-a[:,3]-1))))
  for element,initial in [(0,initialMg),(1,1e-6)]:maximum['elementInventory']=max(maximum['elementInventory'],float(np.max(np.abs(a[:,4+element]+a[:,6+element]+a[:,8+element]-a[:,10+element]-a[:,12+element]-initial))))
  nodes,weights=np.polynomial.legendre.leggauss(384)
  for t in [.151,.8,3.71,7.93]:
   lo=.15;hi=t;delay=(hi-lo)/2*nodes+(hi+lo)/2;A=.1/(.15**-.1-15**-.1);independent=(hi-lo)/2*np.dot(weights,[histories(t-u,p)[0]*A*u**-1.1 for u in delay]);q=response(t,p);maximum['convolutionRelative']=max(maximum['convolutionRelative'],abs(independent-q)/max(abs(q),1e-14))
  for t in [.2,1,8]:
   source=quad(lambda u:response(u,p),.15,t,points=[x for x in [3.7,4]if .15<x<t],epsabs=1e-10)[0];fubini=quad(lambda birth:histories(birth,p)[0]*delay_cdf(t-birth),0,t-.15,points=[x for x in [3.3,3.7,4.1]if 0<x<t-.15],epsabs=1e-10)[0];maximum['sourceIntegral']=max(maximum['sourceIntegral'],abs(source-fubini))
  drops=np.diff(bracket[:,0]);details.append({'name':p['name'],'finalBrackets':bracket[-1].tolist(),'negativeFeHIntervals':int(np.sum(drops<0)),'minimumGas':coarse['minimumGas'],'totalLockedMass':float(a[-1,1])})
 # Prompt-only analytic controls change no admittedAppcoefficients; they are separate mathematical limits.
 wind=[]
 for eta in [0,1]:
  c=.6+eta;psi=.3;y=fixed['netCC']['fe'];z0=1e-6
  answer=solve_ivp(lambda t,u:[-c*psi,y*psi-c*psi*u[1]/u[0]],(0,1),[1,z0],rtol=1e-12,atol=1e-15);G,Mi=answer.y[:,-1];analytic=z0+y/c*math.log(1/G)
  assert abs(Mi/G-analytic)<1e-10
  wind.append({'eta':eta,'sameTimeGas':float(G),'sameTimeIronFraction':float(Mi/G),'sameGasFractionHalfIron':z0+y/c*math.log(2)})
 k=1.;times=np.linspace(0,1,11);mix=solve_ivp(lambda t,z:[k*(z[1]-z[0]),k*(z[0]-z[1])],(0,1),[.002,0],t_eval=times,rtol=1e-12,atol=1e-15).y;mixError=float(np.max(np.abs((mix[0]-mix[1])-.002*np.exp(-2*k*times))));assert mixError<1e-10;assert np.max(np.abs(mix.sum(axis=0)-.002))<1e-12
 # Linear change from independently noisy absoluteFe/Mg toFeH andMg/Fe fixes the negative cross-term.
 transform=np.array([[1,0],[-1,1]]);covariance=transform@np.diag([.1**2,.08**2])@transform.T;rho=covariance[0,1]/math.sqrt(covariance[0,0]*covariance[1,1]);inverse=np.linalg.inv(covariance);det=np.linalg.det(covariance)
 integral=dblquad(lambda y,x:math.exp(-.5*np.array([x,y])@inverse@np.array([x,y]))/(2*math.pi*math.sqrt(det)),-1,1,lambda _: -1,lambda _:1,epsabs=1e-9,epsrel=1e-9)[0];assert abs(integral-1)<1e-8
 rejections={}
 for name,parameters in [('domain',{**family()[0],'eta':2}),('gas',{**family()[0],'formedMass':.55,'eta':1.5,'inflowPeak':0})]:
  try:simulate(parameters,fixed);rejections[name]=False
  except ValueError:rejections[name]=True
 for invalid in [-1,1,float('nan')]:
  try:log_probabilities(coarse,datasets[0][1],invalid);rejections['rho'+str(invalid)]=False
  except ValueError:rejections['rho'+str(invalid)]=True
 assert all(rejections.values())
 checks={'negativeControls':all(rejections.values()),'bulk':maximum['bulk']<1e-6,'elements':maximum['elementScaled']<.001,'abundances':maximum['abundanceDex']<.001,'gasAnalytic':maximum['gasAnalytic']<1e-10,'gasInventory':maximum['gasInventory']<1e-10,'elementInventory':maximum['elementInventory']<1e-11,'mixtureRefinement':maximum['mixtureMeanLogDensity']<.01,'convolution':maximum['convolutionRelative']<.001,'sourceIntegral':maximum['sourceIntegral']<1e-8,'DTD':abs(quad(lambda t:.1*t**-1.1/(.15**-.1-15**-.1),.15,15)[0]-1)<1e-10}
 checks={key:bool(value)for key,value in checks.items()};maximum={key:float(value)for key,value in maximum.items()}
 result={'status':'PASS'if all(checks.values())else'FAIL','checks':checks,'refusedNegativeControls':rejections,'maximum':maximum,'modelSummaries':details,'windControls':wind,'delayFractionAt8Gyr':delay_cdf(8),'covarianceExample':{'matrix':covariance.tolist(),'ratioSigma':math.sqrt(covariance[1,1]),'rho':rho,'independentIntegral':integral},'mixingLimit':{'ratePerGyr':k,'contrastAtHalfGyr':.002/math.e,'maximumAbsoluteDifference':mixError},'realSampleIncludedInRefinement':bool(args.real_data),'referenceArchivePolicy':'Fine full state actually solved; archived fine cohorts and per-field maxima, not full fine bulk arrays','scope':'Finite model source/state refinement, independent mathematical limits and a discrete mixture definition; no observedpopulationprecisioncalibration'};write_json(args.output,result);print(result['status'],maximum)
 if result['status']!='PASS':raise SystemExit(1)
if __name__=='__main__':main()