from pathlib import Path
import pandas as pd
import numpy as np

ROOT = Path(__file__).resolve().parents[1]
D = ROOT / 'derived'

def med(x): return float(np.median(np.asarray(x, dtype=float)))

def section(name):
    print('\n' + '='*72)
    print(name)
    print('='*72)

section('2011 primary raw target')
t = pd.read_csv(D/'URUANA_2011_RAW_GRADIENT_TARGET_LINES.csv')
print('n_lines =', len(t))
print('raw negative =', int((t.raw_slope_nT_km < 0).sum()), '/', len(t))
print('IGRF positive =', int((t.igrf_slope_nT_km > 0).sum()), '/', len(t))
print('median raw slope = %.6f nT/km' % med(t.raw_slope_nT_km))
print('median IGRF slope = %.6f nT/km' % med(t.igrf_slope_nT_km))
assert len(t)==14 and (t.raw_slope_nT_km<0).all() and (t.igrf_slope_nT_km>0).all()

section('2011 same-line self-control / heading / dates')
s = pd.read_csv(D/'URUANA_2011_HEADING_DATE_SELF_CONTROL.csv')
print('target negative =', int((s.raw_slope_nT_km<0).sum()), '/', len(s))
print('north control positive =', int((s.north_MAGBRU_slope>0).sum()), '/', len(s))
print('median target raw = %.6f' % med(s.raw_slope_nT_km))
print('median north raw = %.6f' % med(s.north_MAGBRU_slope))
print('median paired change = %.6f' % med(s.raw_slope_change))
print('minimum paired change = %.6f' % float(s.raw_slope_change.min()))
print('headings =', s.heading_acquisition.value_counts().to_dict())
print('dates =', sorted(s.date.astype(str).unique().tolist()))
assert len(s)==14 and (s.raw_slope_nT_km<0).all() and (s.north_MAGBRU_slope>0).all() and (s.raw_slope_change>0).all()


section('2011 secondary channel and bin-phase/height sensitivity')
print('MAGMIC target negative =', int((s.mic_slope_nT_km<0).sum()), '/', len(s))
print('MAGMIC north positive =', int((s.north_MAGMIC_slope>0).sum()), '/', len(s))
print('median MAGMIC target = %.6f' % med(s.mic_slope_nT_km))
print('median MAGMIC north = %.6f' % med(s.north_MAGMIC_slope))
mic_change = s.north_MAGMIC_slope - s.mic_slope_nT_km
print('median/min MAGMIC paired change = %.6f / %.6f' % (med(mic_change), float(mic_change.min())))
assert (s.mic_slope_nT_km<0).all() and (s.north_MAGMIC_slope>0).all() and (mic_change>0).all()

bp = pd.read_csv(D/'URUANA_2011_BIN_PHASE_CHANNEL_HEIGHT_SENSITIVITY_SUMMARY.csv')
for channel in ['MAGCOM','MAGNIV','MAGMIC']:
    z=bp[(bp.channel==channel)&(bp.estimator=='ols')]
    print(channel, 'target negative counts by phase =', z[z.region=='target'].sort_values('phase_m').n_negative.tolist(),
          '; north positive counts =', z[z.region=='north'].sort_values('phase_m').n_positive.tolist())
    assert (z[z.region=='target'].n_negative==14).all() and (z[z.region=='north'].n_positive==14).all()
for channel in ['MAGBRU','MAGCOR']:
    z=bp[(bp.channel==channel)&(bp.estimator=='ols')]
    print(channel, 'target negative counts by phase =', z[z.region=='target'].sort_values('phase_m').n_negative.tolist(),
          '; north positive counts =', z[z.region=='north'].sort_values('phase_m').n_positive.tolist())
    assert z[z.region=='target'].n_negative.min()==13 and z[z.region=='target'].n_negative.max()==14
    assert (z[z.region=='north'].n_positive==14).all()
expected_height={'GPSALT':13,'BARO':13,'ALTURA':4,'MDT':10}
for channel, npos in expected_height.items():
    z=bp[(bp.region=='target')&(bp.channel==channel)&(bp.estimator=='ols')]
    print(channel, 'target positive counts by phase =', z.sort_values('phase_m').n_positive.tolist())
    assert (z.n_positive==npos).all()

section('2011 mirrored raw controls')
c = pd.read_csv(D/'URUANA_2011_RAW_GRADIENT_MIRRORED_CONTROLS.csv')
target = c.loc[c.k==0].iloc[0]
non = c.loc[c.k!=0]
print('target inversion fraction =', target.raw_inversion_fraction)
print('max non-target inversion fraction =', non.raw_inversion_fraction.max())
print('target median raw slope =', target.median_raw_slope_nT_km)
print('most-negative non-target median raw slope =', non.median_raw_slope_nT_km.min())

section('Municipal-filtered transition controls / matched support')
m = pd.read_csv(D/'URUANA_MUNICIPAL_FILTERED_CONTROL_SUMMARY.csv')
u = m.loc[m.k==0].iloc[0]
ext = m.loc[m.municipality=='outside']
print('Uruana raw transition count = %d/%d' % (u.raw_2011_negpos_count, u.raw_2011_n))
print('max external count = %d/%d' % (ext.raw_2011_negpos_count.max(), int(ext.raw_2011_n.max())))
bearing = m.loc[m.transition_all_3_windows==True].copy()
print('transition-bearing blocks =', len(bearing))
print('Uruana Dmed = %.6f km; Cmean = %.6f' % (u.matched_median_front_distance_km, u.matched_mean_sign_agreement))
print('rank Dmed (ascending) =', int(bearing.matched_median_front_distance_km.rank(method='min').loc[bearing.k==0].iloc[0]), '/', len(bearing))
print('rank Cmean (descending) =', int(bearing.matched_mean_sign_agreement.rank(method='min', ascending=False).loc[bearing.k==0].iloc[0]), '/', len(bearing))
assert u.raw_2011_negpos_count==14 and ext.raw_2011_negpos_count.max()<=9

section('1971-72 historical scalar non-inversion')
h = pd.read_csv(D/'URUANA_1009_HISTORICAL_GRADIENT_HOLDOUT.csv')
print('residual slope range = [%.6f, %.6f]' % (h.residual_slope_nT_per_km.min(), h.residual_slope_nT_per_km.max()))
print('IGRF slope range = [%.6f, %.6f]' % (h.igrf_slope_nT_per_km.min(), h.igrf_slope_nT_per_km.max()))
print('reconstructed total range = [%.6f, %.6f]' % (h.reconstructed_total_slope_nT_per_km.min(), h.reconstructed_total_slope_nT_per_km.max()))
assert (h.residual_slope_nT_per_km<0).all() and (h.igrf_slope_nT_per_km>0).all() and (h.reconstructed_total_slope_nT_per_km>0).all()

section('1971-72 1DV robustness')
v = pd.read_csv(D/'URUANA_1009_NUMERIC_1DV_ROBUSTNESS_450.csv')
# detect nearest-zero column robustly
cols=[x for x in v.columns if 'zero' in x.lower() and ('dist' in x.lower() or 'km' in x.lower())]
if not cols:
    cols=[x for x in v.columns if 'nearest' in x.lower()]
zc=cols[0]
print('rows =',len(v),'zero distance column =',zc)
print('<=8 km =',int((v[zc]<=8).sum()),'/',len(v))
print('min/median/max =',float(v[zc].min()),med(v[zc]),float(v[zc].max()))
assert len(v)==450 and (v[zc]<=8).all()

section('2011 matched-resolution stress')
r = pd.read_csv(D/'URUANA_2011_MATCHED_RESOLUTION_TOTAL_GRADIENT_Y6_8.csv')
for fwhm,rep in [(2,'q2.5'),(6,'q2.5')]:
    z=r[(r.fwhm_km==fwhm)&(r.representation==rep)]
    if len(z): print('FWHM %s km / %s: inverted %d/%d; median slope %.6f' % (fwhm,rep,z.inverted.astype(bool).sum(),len(z),med(z.slope)))

section('Modern scale addition')
e = pd.read_csv(D/'URUANA_EMM2017_RESOLUTION_LADDER_SUMMARY.csv')
print(e.to_string(index=False))

print('\nALL ASSERTED VERIFICATION CHECKS PASSED')