"""Fresh source fields -> historical rate; native profiles -> costs -> endpoints.

The scanned1993ephemeris remains a manual reading check, not an OCR producer or
raw RV refit. Every native row stays local; only aggregate results are written.
"""
from pathlib import Path
import argparse
import hashlib
import json
import re
import subprocess
import sys
import numpy as np
from observed_reference import beta_lyrae_conditional
from envelope_reference import envelope_cost, ce_final_separation
from roche_reference import volume_radius
from orbit_reference import integrate_transfer, native_mesa_closed


def verify_inputs(inputs):
    for source in json.loads(Path('source-manifest.json').read_text())['sources']:
        body = (inputs/source['name']).read_bytes()
        if len(body) != source['bytes'] or hashlib.sha256(body).hexdigest() != source['sha256']:
            raise ValueError('Input origin changed: '+source['name'])


def native_software_reference(inputs):
    text = (inputs/'mesa-jdot-inlist').read_text()
    def field(name):
        found = re.search(r'^\s*'+re.escape(name)+r'\s*=\s*([-+0-9.eEdD]+)', text, re.M)
        if not found: raise ValueError('Native inlist field missing: '+name)
        return float(found[1].lower().replace('d', 'e'))
    donor, accretor, separation = field('m1'), field('m2'), field('initial_separation_in_Rsuns')
    wind, reemission, ring, gamma = [field('mass_transfer_'+name) for name in ['alpha', 'beta', 'delta', 'gamma']]
    retention = 1-wind-reemission-ring
    numerical = integrate_transfer(donor, accretor, separation, .65*donor, retention, 'mesa-mixture',
        wind_fraction=wind, reemission_fraction=reemission, ring_fraction=ring, ring_gamma=gamma)
    closed = native_mesa_closed(donor, accretor, separation, numerical['stripped'], wind, reemission, ring, gamma)
    error = float(np.max(abs(numerical['separation']/closed-1))); assert error < 1e-9
    return {'source': 'MESA25.12.1 fixed inlist and mathematical run_binary_extras.f90 reference',
            'fresh_inlist_fields': {'donor': donor, 'accretor': accretor, 'separation': separation, 'wind': wind,
                                   'reemission': reemission, 'ring': ring, 'ring_lever_arm': gamma, 'retained': retention},
            'maximum_relative_ODE_vs_closed_error': error,
            'boundary': 'Published/native mathematical formula checked; no MESA compilation, stellar evolution run, or translation of its Fortran into the Engine.'}


def reproduce(inputs, output):
    verify_inputs(inputs); output.mkdir(parents=True, exist_ok=True)
    beta = beta_lyrae_conditional(inputs/'beta-lyrae-2013-fulltext.html')
    assert abs(beta['conservative_rate_solar_per_year']/2.179551062119472e-5-1) < 1e-9
    (output/'beta-lyrae-summary.json').write_text(json.dumps(beta, indent=2)+'\n')
    (output/'beta-loss-scenarios.json').write_text(json.dumps(beta['conditional_loss_channel_scenarios'], indent=2)+'\n')
    structures = []; aggregate = []; endpoints = []
    members = json.loads(Path('native-members.json').read_text())
    for member in members:
        profile = inputs/member['profile_file']
        if hashlib.sha256(profile.read_bytes()).hexdigest() != member['sha256']: raise ValueError('Native profile hash changed')
        for threshold in [.001, .01, .1]:
            cost = envelope_cost(profile, threshold)
            half = envelope_cost(profile, threshold, 2)
            cost['stride_two_relative_change'] = half['gravity_cost_units']/cost['gravity_cost_units']-1
            structures.append({'source': member['name'], 'original_sha256': member['sha256'], 'license': member['license'], 'attribution': member['attribution'], **cost})
            donor, core = cost['native_mass_solar'], cost['core_mass_solar']; companion = 1.4
            initial = cost['native_photospheric_radius_solar']/volume_radius(donor/companion, 256)
            core_fraction = volume_radius(core/companion, 256); companion_fraction = volume_radius(companion/core, 256)
            for alpha in [.1, .25, .5, 1.]:
                final = ce_final_separation(donor, core, companion, initial, cost['gravity_cost_units'], alpha)
                overlap = final <= cost['cut_radius_solar']+2e-5
                fits = cost['cut_radius_solar'] < final*core_fraction and 2e-5 < final*companion_fraction
                endpoints.append({'source': member['name'], 'hydrogen_threshold': threshold, 'alpha': alpha,
                                  'initial_separation_solar': initial, 'final_separation_solar': final,
                                  'direct_overlap': overlap, 'both_supplied_radii_fit': fits,
                                  'event': 'supplied-surfaces-overlap' if overlap else 'detached-with-supplied-radii' if fits else 'supplied-surface-outside-roche-lobe'})
            if threshold == .01:
                aggregate.append({'id': 'core-he-4m' if donor < 5 else 'rsg-15m', 'massSolar': donor,
                    'radiusSolar': cost['native_photospheric_radius_solar'], 'coreMassSolar': core,
                    'coreCutRadiusSolar': cost['cut_radius_solar'], 'gravitationalCostUnits': cost['gravity_cost_units'],
                    'gravityLambda': cost['lambda_gravity'], 'license': member['license'], 'attribution': member['attribution']})
    aggregate.sort(key=lambda row: row['massSolar'])
    (output/'structure-summary.json').write_text(json.dumps(structures, indent=2)+'\n')
    (output/'structure-aggregate.json').write_text(json.dumps(aggregate, indent=2)+'\n')
    (output/'envelope-endpoints.json').write_text(json.dumps(endpoints, indent=2)+'\n')
    native = native_software_reference(inputs)
    summary = {'status': 'PASS freshly parsed historical Table1/period trend and native stellar columns; costs/endpoints/ODE reference recalculated',
        'historical_conservative_rate_solar_per_year': beta['conservative_rate_solar_per_year'],
        'loss_scenarios': len(beta['conditional_loss_channel_scenarios']),
        'no_positive_rate_cases': sum(not row['positive_rate_solution'] for row in beta['conditional_loss_channel_scenarios']),
        'native_profile_cuts': len(structures), 'conditional_endpoints': len(endpoints), 'native_mathematical_reference': native,
        'boundary': 'Historical adopted-parameter replay, not raw RV/light-curve fit, current rate, full posterior or CE alpha measurement. Native structure aggregation is not a new MESA evolution. No relaxed core radius, internal/recombination energy or hydrodynamic envelope outcome inferred. Scanned1993PDF is manual-reading input only.'}
    (output/'summary.json').write_text(json.dumps(summary, indent=2)+'\n')
    subprocess.run([sys.executable, 'produce_figures.py', '--results', str(output)], check=True)
    print(json.dumps(summary, indent=2))


if __name__ == '__main__':
    parser = argparse.ArgumentParser(); parser.add_argument('--inputs', type=Path, default=Path('inputs')); parser.add_argument('--output', type=Path, default=Path('results')); args = parser.parse_args(); reproduce(args.inputs, args.output)