#!/usr/bin/env python3
"""Independent fit and alert replay of saved synthetic exposures, Python/NumPy.

No Lab import, browser or network. Does not regenerate source/noise. The observed
OGLE replay has its own inputs, full-grid null calibration and injection trials.
"""
import argparse
import hashlib
import json
import math
from pathlib import Path
import numpy as np

HERE = Path(__file__).resolve().parent
ROOT = HERE.parents[2] if len(HERE.parents) >= 3 else HERE
REPOSITORY_INPUT = ROOT/'docs/astronomy/ac07/synthetic-results.json'
DEFAULT_INPUT = REPOSITORY_INPUT if REPOSITORY_INPUT.exists() else HERE/'synthetic-results.json'


def analyze(entry):
    points = entry['points']
    time = np.array([point['timeDays'] for point in points])
    values = np.array([point['value'] for point in points])
    sigma = np.array([point['sigma'] for point in points])
    centered = (values - np.average(values, weights=1/sigma**2)) / sigma
    reference = np.dot(centered, centered)
    best_frequency, best_power = .1, -1
    for frequency in .1 + .005 * np.arange(981):
        phase = 2*np.pi*frequency*time
        design = np.column_stack([np.ones(len(time)), np.cos(phase), np.sin(phase)]) / sigma[:, None]
        coefficients = np.linalg.lstsq(design, values / sigma, rcond=None)[0]
        residual = values/sigma - design @ coefficients
        power = 1 - np.dot(residual, residual) / reference
        if power > best_power + 1e-12:
            best_frequency, best_power = float(frequency), float(power)
    window = float(abs(np.exp(2j*np.pi*time).mean())**2)
    detected = [point for point in entry['transientPoints'] if point['value'] >= 5*point['sigma']]
    first_alert = detected[0]['timeDays'] + (2.5 + entry['settings']['processingMinutes'])/1440 if detected else None
    crossing = 8.68 - (.7/24)*math.sqrt(2*math.log(1.2/(5*.035)))
    result = {'preset': entry['preset'], 'exposures': len(points), 'best_period_days': 1/best_frequency, 'peak_power': best_power, 'daily_window_power': window, 'detected_transient_exposures': len(detected), 'first_alert_days': first_alert, 'alert_delay_from_reference_threshold_hours': None if first_alert is None else (first_alert-crossing)*24}
    if abs(result['best_period_days']-entry['observation']['bestPeriodDays']) > 1e-9:
        raise ValueError('Independent least-squares replay disagrees with saved period')
    return result


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--input', type=Path, default=DEFAULT_INPUT)
    parser.add_argument('--out', type=Path, default=ROOT/'docs/astronomy/ac07/synthetic-replay.json' if REPOSITORY_INPUT.exists() else HERE/'synthetic-replay-rerun.json')
    args = parser.parse_args()
    raw = args.input.read_bytes()
    data = json.loads(raw)
    result = {'input_sha256': hashlib.sha256(raw).hexdigest(), 'kind': 'independent reanalysis of saved synthetic exposures; not OGLE measurements', 'numpy_version': np.__version__, 'method': 'weighted floating-mean sine/cosine/constant least squares using NumPy SVD; first tied maximum within 1e-12; 0.1..5/day step0.005; five-sigma alert at exposure end plus processing delay', 'results': [analyze(entry) for entry in data['results']]}
    args.out.parent.mkdir(parents=True, exist_ok=True)
    args.out.write_text(json.dumps(result, indent=2, allow_nan=False)+'\n')
    print(json.dumps(result['results'], indent=2))


if __name__ == '__main__':
    main()