"""Saved-observation failures and aggregate whole-cluster conditional checks."""
from pathlib import Path
import json
import numpy as np
from scipy.special import ndtr
from period_likelihood import HarmonicLikelihood
from frequency_measure import frequency_measure, period_posterior
from empirical_calibration import EmpiricalCalibration
from conditional_age import age_from_frequency_measure, age_from_reported_log_period
from rotation_reference import flux_series
from inspect_cluster_samples import read_benchmarks

rng = np.random.default_rng(271104)
times = np.sort(rng.uniform(0, 60, 180))
cases = [('coherent', {}), ('short-lifetime', {'lifetime_days': 3., 'peak_day': 30.}),
         ('pole-on', {'inclination_degrees': 0.}), ('differential-latitude', {'latitude_degrees': 60., 'shear': .2}),
         ('two-longitudes', {'longitudes': (0., np.pi)})]
results = []
for name, parameters in cases:
    clean, truth = flux_series(times, **parameters)
    observed = clean + rng.normal(0, .0005, len(times))
    likelihood = HarmonicLikelihood(times, observed, .0005)
    measure = frequency_measure(likelihood)
    assert measure['status'] == 'integrated', measure
    posterior = period_posterior(measure, likelihood)
    flags = posterior['flags']
    result = {'case': name, 'flux_range': float(np.ptp(clean)), 'leading_pattern_period_days': float(1 / measure['peaks'][0]['frequency']),
              'injection_local_period_for_validation_only_days': truth, 'flags': flags,
              'reduced_chi_square': posterior['best_reduced_chi_square'], 'log_bayes_factor': posterior['log_bayes_factor_vs_constant_trend'],
              'age_policy': 'withheld because saved-observation likelihood failed' if flags else 'fixed-scale mixture conditional density only',
              'evidence_evaluations': measure['actual_evidence_evaluations']}
    if name in ['short-lifetime', 'pole-on']:
        assert flags, result
    else:
        assert not flags, result
    results.append(result)
    print(result, flush=True)
# An exact daily window preserves genuine equal-height aliases; no period-control cut hides them.
daily_times = np.arange(61, dtype=float)
daily_clean, _ = flux_series(daily_times)
daily = HarmonicLikelihood(daily_times, daily_clean + rng.normal(0, .0005, len(daily_times)), .0005)
exact = daily.evaluate([.1, 1.1, 2.1, 3.1])
assert np.ptp(exact) < 1e-6, exact
measure = frequency_measure(daily)
assert measure['status'] == 'integrated', measure
posterior = period_posterior(measure, daily)
results.append({'case': 'daily-alias', 'exact_alias_log_evidence_range': float(np.ptp(exact)),
                'retained_likelihood_peak_count': len(measure['peaks']), 'leading_pattern_periods_days': [float(1 / row['frequency']) for row in measure['peaks'][:12]],
                'age_policy': 'Retain all supported aliases and one/two-pattern hypotheses; never select the injected period'})
# Proper generative model prior-predictive check: coverage belongs only to this declared fixed model.
calibration = EmpiricalCalibration(Path('inputs/curtis2020-source/Benchmark_Clusters.txt'))
pit = []; covered = []
for _ in range(400):
    age = rng.uniform(*calibration.age_bounds); temperature = rng.uniform(4400, 5700)
    mean, width = calibration.sequence([temperature], [age])
    log_period = rng.normal(mean[0, 0], width[0, 0])
    posterior = age_from_reported_log_period(calibration, log_period, 0., temperature, sigma_temperature=0., age_points=801)
    quantiles = posterior['quantiles025_16_50_84_975_gyr']
    covered.append([quantiles[0] <= age <= quantiles[4], quantiles[1] <= age <= quantiles[3]])
    pit.append(float(np.interp(age, posterior['age_gyr'], posterior['cdf'])))
coverage = np.mean(covered, axis=0)
ordered = np.sort(pit); ks = max(np.max(np.arange(1, 401) / 400 - ordered), np.max(ordered - np.arange(400) / 400))
assert .91 < coverage[0] < .98 and .60 < coverage[1] < .76 and ks < .09, (coverage, ks)
# Only aggregate diagnostics leave private source cache. Reported periods have no supplied errors.
stars = read_benchmarks(Path('inputs/curtis2020-source/Benchmark_Clusters.txt'))
holdouts = []
for cluster, model in [('NGC 6819', calibration), ('NGC 752', calibration), ('NGC 6811', EmpiricalCalibration(Path('inputs/curtis2020-source/Benchmark_Clusters.txt'), excluded_cluster='NGC 6811'))]:
    selected = [row for row in stars if row['cluster'] == cluster and model.temperature_bounds[0] <= row['teff_k'] <= model.temperature_bounds[1]]
    summaries = []
    for row in selected:
        posterior = age_from_reported_log_period(model, np.log(row['period_days']), 0., row['teff_k'], sigma_temperature=0.)
        quantiles = posterior['quantiles025_16_50_84_975_gyr']
        summaries.append([quantiles[2], quantiles[4] - quantiles[0], quantiles[0] <= row['age_gyr'] <= quantiles[4], bool(posterior['flags'])])
    assert selected
    summaries = np.array(summaries)
    holdouts.append({'whole_cluster': cluster, 'stars_in_actual_support': len(selected), 'fixed_source_age_gyr': selected[0]['age_gyr'],
                     'median_conditional_age_gyr': float(np.median(summaries[:, 0])), 'median_95_percent_width_gyr': float(np.median(summaries[:, 1])),
                     'fraction_including_fixed_source_age': float(np.mean(summaries[:, 2])), 'fraction_boundary_notices': float(np.mean(summaries[:, 3])),
                     'conditioning': 'Exact reported catalog period and temperature; no missing measurement errors invented. Source age is held fixed, not independently certified truth.'})
report = {'status': 'PASS retained observation failures, proper generative normalization and aggregate whole-cluster audit',
          'saved_flux_cases': results, 'prior_predictive_fixed_model': {'draws': 400, 'coverage95_68': coverage.tolist(), 'cdf_uniformity_ks': float(ks),
          'qualification': 'Model-generated prior-predictive check, not field-star coverage or absolute age accuracy'}, 'whole_cluster_holdouts': holdouts,
          'limits': 'Cluster holdouts retain common source extinction/membership/period-selection assumptions. No author rows or star scatter are published; boundary accumulations remain warnings.'}
Path('results/failure-domain-results.json').write_text(json.dumps(report, indent=2) + '\n')
print(json.dumps(report, indent=2))