"""Independent quadrature/rejection/flux/projection references and retained failures."""
import math
import json
import numpy as np
from scipy.integrate import quad
from scipy.stats import kstest
from .population import sample_imf, inverse_power, combine_magnitudes, sample_plummer, pair_components
from .inputs import model_at_age


def validate_limits(output=None):
    rng=np.random.default_rng(20261104)
    slopes=[1.,1.3,2.3,3.2]
    checks=[]
    for slope in slopes:
        normalization=quad(lambda mass:mass**(-slope),.5,8,epsabs=1e-12)[0]
        uniforms=np.linspace(.0001,.9999,501)
        masses=inverse_power(uniforms,.5,8,slope)
        cdf=np.array([quad(lambda mass:mass**(-slope),.5,mass,epsabs=1e-12)[0]/normalization for mass in masses])
        maximum=float(np.max(np.abs(cdf-uniforms)));assert maximum<2e-12
        checks.append({'slope':slope,'maximumIndependentCDFError':maximum})
    def density(mass):
        return np.where(np.asarray(mass)<.5,np.asarray(mass)**(-1.3),.5**(2.3-1.3)*np.asarray(mass)**(-2.3))
    normalization=quad(lambda mass:float(density(mass)),.1,.5)[0]+quad(lambda mass:float(density(mass)),.5,8)[0]
    inverse=sample_imf(rng.random(100000))
    expected_low=quad(lambda mass:float(density(mass)),.1,.5)[0]/normalization
    counts=int((inverse<.5).sum());expected=100000*expected_low
    assert abs(counts-expected)<5*math.sqrt(100000*expected_low*(1-expected_low))
    # Independent rejection draw in log mass, not the inverse-CDF implementation.
    accepted=[]
    maximum=max(.1*float(density(.1)),.5*float(density(.5)))
    while sum(len(chunk) for chunk in accepted)<20000:
        candidate=np.exp(rng.uniform(np.log(.1),np.log(8),10000))
        accepted.append(candidate[rng.random(len(candidate))*maximum<candidate*density(candidate)])
    rejection=np.concatenate(accepted)[:20000]
    from scipy.stats import ks_2samp
    comparison=ks_2samp(inverse,rejection);assert comparison.statistic<.02
    first=rng.uniform(-8,19,(300,3));second=rng.uniform(-8,19,(300,3))
    combined=combine_magnitudes(first,second)
    independent=np.array([[-2.5*math.log10(math.fsum([10**(-.4*a),10**(-.4*b)])) for a,b in zip(x,y)] for x,y in zip(first,second)])
    flux_error=float(np.max(np.abs(combined-independent)));assert flux_error<2e-13
    equal=combine_magnitudes(first,first)-first
    equal_error=float(np.max(np.abs(equal+2.5*np.log10(2))));assert equal_error<1e-13
    # Pairing keeps the component list and total initial mass exactly; f is per system.
    pairing=[]
    for fraction in [0,.2,.5,1]:
        systems=pair_components(inverse[:1024],fraction,rng)
        identities=np.concatenate(systems)
        assert sorted(identities.tolist())==list(range(1024))
        pairing.append({'requestedSystemBinaryFraction':fraction,'realizedFraction':sum(len(s)==2 for s in systems)/len(systems),'systems':len(systems),'components':len(identities)})
    before_projection_state=rng.bit_generator.state
    radial_uniforms=rng.uniform(1e-10,1-1e-10,200000)
    direction_cosines=rng.uniform(-1,1,200000)
    azimuths=rng.uniform(0,2*np.pi,200000)
    positions=sample_plummer(radial_uniforms,direction_cosines,azimuths)
    if output is not None:
        np.savez_compressed(output/'original-projection-inputs.npz',radial_uniforms=radial_uniforms,direction_cosines=direction_cosines,azimuths=azimuths,positions=positions)
        (output/'original-projection-rng-state.json').write_text(json.dumps(before_projection_state,indent=2)+'\n')
    radii=np.linalg.norm(positions,axis=1)
    inverse_radius_error=float(np.max(np.abs(radii**3/(1+radii**2)**1.5-radial_uniforms)))
    assert inverse_radius_error<2e-14
    projected=np.hypot(positions[:,0],positions[:,1])
    projection=[]
    for radius in [.1,.5,1,2,5]:
        # Independently integrate the three-dimensional density along the line of sight.
        surface=quad(lambda z:2*3/(4*np.pi)*(1+radius**2+z*z)**(-2.5),0,np.inf,epsabs=1e-12)[0]
        analytic_surface=1/np.pi/(1+radius**2)**2
        assert abs(surface-analytic_surface)<2e-12
        analytic_cdf=radius**2/(1+radius**2)
        sampled_cdf=float((projected<=radius).mean())
        projection.append({'radiusOverScale':radius,'independentSurfaceDensity':surface,'analyticSurfaceDensity':analytic_surface,'sampledProjectedCDF':sampled_cdf,'analyticProjectedCDF':analytic_cdf})
    ks_projection=kstest(projected,lambda radius:radius**2/(1+radius**2))
    dkw_error=math.sqrt(math.log(2/.001)/(2*len(projected)))
    # Retain the initial arbitrary 0.003 point tolerance failure; do not reroll.
    initial_failures=[row for row in projection if abs(row['sampledProjectedCDF']-row['analyticProjectedCDF'])>=.003]
    plummer_potential=quad(lambda radius:-3*radius**4/(1+radius**2)**4,0,np.inf,epsabs=1e-12)[0]
    assert abs(plummer_potential+3*np.pi/32)<1e-12
    spatial_normalization=quad(lambda radius:3*radius**2/(1+radius**2)**2.5,0,np.inf,epsabs=1e-12)[0]
    projected_normalization=quad(lambda radius:2*radius/(1+radius**2)**2,0,np.inf,epsabs=1e-12)[0]
    assert abs(spatial_normalization-1)<1e-12 and abs(projected_normalization-1)<1e-12
    assert abs(5*np.log10(1000/10)-10)<1e-14
    return {'status':'PASS deterministic mathematical references; original seeded sampling gate separately failed and retained','powerLawCDF':checks,'imfLowMassCount':counts,
            'independentExpectedCount':expected,'inverseVsRejectionKS':float(comparison.statistic),
            'binaryMaximumMagnitudeError':flux_error,'equalMassBinaryMaximumError':equal_error,
            'equalMassBrighteningMag':float(-2.5*np.log10(2)),'pairing':pairing,'projection':projection,
            'inversePlummerRadialCDFError':inverse_radius_error,
            'retainedInitialPointToleranceFailures':initial_failures,
            'projectionSimultaneousSamplingCheck':{'n':len(projected),'confidence':.999,'dkwSupremumBand':dkw_error,'measuredKS':float(ks_projection.statistic),'passesStatedBand':bool(ks_projection.statistic<dkw_error),'notDeterministicCertificate':True},
            'independentPlummerPotentialGMaUnits':plummer_potential,
            'independentSpatialNumberNormalization':spatial_normalization,'independentProjectedNumberNormalization':projected_normalization,
            'boundary':'Fixed seeded finite checks, not an empirical universal IMF or cluster-dynamics validation.'}


def validate_grid_interpolation(grids):
    low=model_at_age(grids[0.],9.4);high=model_at_age(grids[0.],9.6);middle=model_at_age(grids[0.],9.5)
    masses=np.linspace(.15,1.4,501);errors=[];phase_changes=0
    for mass in masses:
        if mass>min(low['Mini'][-1],high['Mini'][-1],middle['Mini'][-1]):continue
        def phase(grid):
            lower=np.searchsorted(grid['Mini'],mass,side='right')-1
            upper=min(lower+1,len(grid['Mini'])-1)
            return int(grid['label'][lower]) if grid['label'][lower]==grid['label'][upper] else None
        if None in [phase(grid) for grid in [low,high,middle]] or len({phase(grid) for grid in [low,high,middle]})>1:
            phase_changes+=1;continue
        row={'mass':float(mass),'phase':phase(middle)}
        for band in ['Gmag','G_BPmag','G_RPmag']:
            values=[np.interp(mass,grid['Mini'],grid[band]) for grid in [low,high,middle]]
            interpolated=-2.5*np.log10((10**(-.4*values[0])+10**(-.4*values[1]))/2)
            row[band+'Error']=float(interpolated-values[2])
        errors.append(row)
    # Av .1 is withheld from a source-native 0/.2 interpolation at equal mass.
    extinction_errors=[]
    for age in np.unique(grids[0.]['logAge']):
        a,b,c=[model_at_age(grids[av],age) for av in [0.,.1,.2]]
        for mass in np.linspace(.15,min(1.2,a['Mini'][-1]),100):
            for band in ['Gmag','G_BPmag','G_RPmag']:
                values=[np.interp(mass,grid['Mini'],grid[band]) for grid in [a,b,c]]
                extinction_errors.append(abs((values[0]+values[2])/2-values[1]))
    worst=max(abs(row[band+'Error']) for row in errors for band in ['Gmag','G_BPmag','G_RPmag'])
    return {'withheldLogAge':9.5,'brackets':[9.4,9.6],'comparedSamePhaseMasses':len(errors),'phaseBoundaryRefusals':phase_changes,
            'maximumSamePhaseMagnitudeError':worst,'rows':errors,
            'sourceExtinctionWithheldAv':.1,'maximumAvInterpolationErrorMag':max(extinction_errors),
            'continuousAgeInterpolationAccepted':worst<.02 and phase_changes==0,
            'fallbackContract':'Use exact native age nodes where age interpolation fails or phases/support differ. Do not call the age grid a continuous certified stellar clock.'}