"""Descriptive CMD distances and native-grid discrepancies, never an age posterior."""
import numpy as np
from .inputs import model_at_age


def cmd_distances(observed_colour, observed_g, model, distance_modulus, twins=False):
    colour=model['G_BPmag']-model['G_RPmag']
    g=model['Gmag']+distance_modulus
    # Coordinates in explicitly declared display scales, not measured sigmas.
    coordinates=np.column_stack([colour/.05, g/.2])
    same_phase=model['label'][:-1]==model['label'][1:]
    starts=np.concatenate([coordinates[:-1][same_phase],coordinates])
    ends=np.concatenate([coordinates[1:][same_phase],coordinates])
    if twins:
        twin_coordinates=np.column_stack([colour/.05, (g-2.5*np.log10(2))/.2])
        starts=np.concatenate([starts,twin_coordinates[:-1][same_phase],twin_coordinates])
        ends=np.concatenate([ends,twin_coordinates[1:][same_phase],twin_coordinates])
    vectors=ends-starts;squared_lengths=np.sum(vectors*vectors,axis=1)
    displacement=np.column_stack([observed_colour/.05, observed_g/.2])[:,None,:]-starts[None,:,:]
    fraction=np.clip(np.sum(displacement*vectors[None,:,:],axis=2)/np.maximum(squared_lengths,1e-30),0,1)
    # Exact nearest segment on this native polyline, avoiding point-density bias.
    residual=displacement-fraction[:,:,None]*vectors[None,:,:]
    return np.sqrt(np.min(np.sum(residual*residual,axis=2),axis=1))


def compare_cmd(gaia, masks, grids):
    result=[]
    for selection in ['fixed','covarianceSensitivity']:
        selected=masks[selection]
        for window, mask in [('full-selected',selected),('declared-turnoff-window',selected & (gaia['phot_g_mean_mag']>12) & (gaia['phot_g_mean_mag']<16.5))]:
            colour=gaia['phot_bp_mean_mag'][mask]-gaia['phot_rp_mean_mag'][mask]
            g=gaia['phot_g_mean_mag'][mask]
            for av,grid in grids.items():
                for age in sorted(set(grid['logAge'])):
                    model=model_at_age(grid,age)
                    # Omit only model phases outside the first-slice contract.
                    supported=(model['label']>=0)&(model['label']<=5)
                    model={k:v[supported] for k,v in model.items()}
                    for modulus in [9.5,9.6,9.7]:
                        distance=cmd_distances(colour,g,model,modulus)
                        result.append({'selection':selection,'window':window,'stars':int(mask.sum()),'av':av,
                                       'logAgeYr':float(age),'distanceModulus':modulus,
                                       'medianDisplayDistance':float(np.median(distance)),
                                       'p90DisplayDistance':float(np.quantile(distance,.9)),
                                       'singleAndTwinUnionMedian':float(np.median(cmd_distances(colour,g,model,modulus,True)))})
    summaries=[]
    for selection in ['fixed','covarianceSensitivity']:
        for window in ['full-selected','declared-turnoff-window']:
            rows=[row for row in result if row['selection']==selection and row['window']==window]
            best=dict(min(rows,key=lambda row:row['medianDisplayDistance']))
            best['touchesAgeGridBoundary']=best['logAgeYr'] in [min(row['logAgeYr'] for row in rows),max(row['logAgeYr'] for row in rows)]
            best['touchesDistanceGridBoundary']=best['distanceModulus'] in [9.5,9.7]
            best['interiorAgeOptimumEstablished']=False
            summaries.append(best)
    return {'coordinateScalesMag':{'BPminusRP':.05,'G':.2},'notChiSquaredOrLikelihood':True,
            'notAnAgePosterior':True,'doesNotInferBinaryFraction':True,
            'selectionConditionedBestDisplayMatches':summaries,'comparisons':result,
            'curveMetric':'Nearest same-phase native-polyline segment or original vertex in declared display coordinates. No line across disconnected evolutionary phases.',
            'boundary':'Aperture, faint cuts and fixed kinematic cuts remain; no Gaia completeness or field model. The single/twin union has extra flexibility and is not model evidence.'}


def compare_native_grids(parsec,mist,native):
    grid=model_at_age(parsec,9.6)
    # Compare main-sequence branches at common initial mass, not arbitrary RGB matches.
    pms=grid['label']==1
    mms=mist['phase']==0
    masses=np.linspace(.5,1.1,61)
    rows=[]
    for mass in masses:
        row={'initialMassMsun':float(mass)}
        for name,pkey,mkey in [('logTeff','logTe','log_Teff'),('logL','logL','log_L'),('G','Gmag','Gaia_G_EDR3')]:
            pvalue=np.interp(mass,grid['Mini'][pms],grid[pkey][pms])
            mvalue=np.interp(mass,mist['initial_mass'][mms],mist[mkey][mms])
            row['parsecMinusMist'+name]=float(pvalue-mvalue)
        rows.append(row)
    # The published M67 model has DR2 magnitudes: compare physical HR columns only.
    keep=(native['center_h1']>1e-3)&(native['initial_mass']<=1.2)
    native_delta=[]
    for mass in np.linspace(.86,1.15,30):
        native_delta.append({'initialMassMsun':float(mass),
          'parsecMinusIsoALogL':float(np.interp(mass,grid['Mini'][pms],grid['logL'][pms])-np.interp(mass,native['initial_mass'][keep],native['log_L'][keep])),
          'parsecMinusIsoALogTeff':float(np.interp(mass,grid['Mini'][pms],grid['logTe'][pms])-np.interp(mass,native['initial_mass'][keep],native['log_Teff'][keep]))})
    return {'sameNominalLogAgeYr':9.6,'parsecZ':float(grid['Zini'][0]),'mistZ':.0178798,
            'sameNominalSolarMetallicityIsNotSameComposition':True,'parsecVsMistMainSequence':rows,
            'publishedIsoAAgeGyr':3.95,'publishedIsoAPhotometry':'Gaia DR2, not DR3; not overlaid on DR3 data',
            'physicalHROnlyVsIsoA':native_delta,
            'boundary':'These independent native outputs expose model/composition/bolometric-correction differences. Agreement would not validate an absolute clock; discrepancies are retained.'}