"""Conditional sensitivity checks; geometry fixed, all station nuisances profiled.

Hypothetical exposure lengths and a common ephemeris timing offset do not supply
missing observing metadata or convert headerless photometry to BJD/TDB.
"""
from pathlib import Path
import argparse
import json
import platform
import numpy as np
from event_model import EventReference, load_observations, profile_magnitudes, projected_earth, source_trajectory


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--materials', type=Path, required=True)
    parser.add_argument('--baseline-results', type=Path, required=True)
    parser.add_argument('--parallax-results', type=Path, required=True)
    parser.add_argument('--output', type=Path, required=True)
    arguments = parser.parse_args()
    folder = arguments.materials
    observations = load_observations(folder / 'observations', json.loads((folder / 'observation-sources.json').read_text()))
    library = folder / ('event-reference.dylib' if platform.system() == 'Darwin' else 'event-reference.so')
    reference = EventReference(library)
    baseline = json.loads(arguments.baseline_results.read_text())['models']
    parallax = json.loads(arguments.parallax_results.read_text())['models']
    models = [min((row for row in baseline if row['model'] == 'wide'), key=lambda row: row['chi_square']),
              min((row for row in parallax if row['model'] == 'wide-annual-parallax'), key=lambda row: row['chi_square'])]
    rows = []
    for model in models:
        parameters = np.array(model['parameters'])
        annual = len(parameters) == 9
        def predict(gamma, exposure_minutes=0, time_shift=0, stride=1):
            nodes, weights = np.polynomial.legendre.leggauss(9 if exposure_minutes else 1)
            result = np.zeros(len(observations['time']))
            for node, weight in zip(nodes, weights / 2):
                times = observations['time'] + node * exposure_minutes / (2 * 1440)
                earth = projected_earth(folder, times, stride=stride, time_shift_days=time_shift) if annual else None
                source = source_trajectory(times, parameters, earth)
                result += weight * reference.binary(source, parameters[3], np.exp(parameters[4]), np.exp(parameters[6]), gamma)
            return result
        original = predict(observations['gamma'])
        original_errors, _, _ = profile_magnitudes(observations, original)
        choices = [('six-hour ephemeris / instantaneous exposure', observations['gamma'], 0, 0, 1)]
        if annual:
            choices += [('twelve-hour ephemeris', observations['gamma'], 0, 0, 2),
                        ('Earth epochs shifted -0.007 day', observations['gamma'], 0, -.007, 1),
                        ('Earth epochs shifted +0.007 day', observations['gamma'], 0, .007, 1)]
        moa = observations['station'] == next(index for index, name in enumerate(observations['station_names']) if 'moa' in name.lower())
        alternate = observations['gamma'].copy(); alternate[moa] = .626
        choices += [('unknown MOA coefficient changed from I-like to R-like', alternate, 0, 0, 1)]
        choices += [('hypothetical common ' + str(minutes) + '-minute rectangular exposure', observations['gamma'], minutes, 0, 1) for minutes in [5, 15, 30]]
        for name, gamma, exposure, shift, stride in choices:
            magnification = predict(gamma, exposure, shift, stride)
            errors, stations, _ = profile_magnitudes(observations, magnification)
            rows.append({'model': model['model'], 'case': name, 'chi_square': float(errors @ errors),
                         'chi_square_change_from_fixed_geometry_baseline': float(errors @ errors - original_errors @ original_errors),
                         'maximum_absolute_source_magnitude_change': float(np.max(abs(2.5 * np.log10(magnification / original)))),
                         'station_nuisances': stations})
    report = {'status': 'COMPLETED conditional sensitivity; lens geometry fixed, station baselines and blend fractions refitted',
              'qualification': 'Five original headerless tables; HJD-like epochs retained. Common +/-0.007-day Earth-epoch shift is a sensitivity probe, not a timing conversion. Exposure lengths are hypothetical; actual values are not supplied by these files. These checks are not a reoptimization of all lens parameters.',
              'cases': rows}
    arguments.output.write_text(json.dumps(report, indent=2) + '\n')
    print(json.dumps({'status': report['status'], 'cases': [{key: row[key] for key in ('model', 'case', 'chi_square_change_from_fixed_geometry_baseline', 'maximum_absolute_source_magnitude_change')} for row in rows]}, indent=2))


if __name__ == '__main__':
    main()