"""Normalize a fixed-scale conditional age density from likelihood, not peak truth."""
import numpy as np
from numpy.polynomial.legendre import leggauss
from scipy.special import ndtr


def temperature_quadrature(calibration, observed_temperature, sigma_temperature, order=40):
    low, high = calibration.temperature_bounds
    if not low <= observed_temperature <= high:
        return {'status': 'outside-support', 'reason': 'Temperature outside the actual common calibrator support'}
    if sigma_temperature == 0:
        return {'status': 'supported', 'temperatures': np.array([observed_temperature]), 'weights': np.array([1.]), 'retained_temperature_likelihood': 1.}
    assert sigma_temperature > 0
    nodes, weights = leggauss(order)
    integration_low = max(low, observed_temperature - 7 * sigma_temperature)
    integration_high = min(high, observed_temperature + 7 * sigma_temperature)
    temperatures = (integration_low + integration_high) / 2 + nodes * (integration_high - integration_low) / 2
    retained = float(ndtr((high - observed_temperature) / sigma_temperature) - ndtr((low - observed_temperature) / sigma_temperature))
    weights = weights * (integration_high - integration_low) / 2 / (sigma_temperature * np.sqrt(2 * np.pi)) * np.exp(-.5 * ((temperatures - observed_temperature) / sigma_temperature) ** 2) / retained
    return {'status': 'supported', 'temperatures': temperatures, 'weights': weights, 'retained_temperature_likelihood': retained}


def normalize_age(calibration, ages, likelihood, prior='uniform', retained_temperature_likelihood=1.):
    assert prior in ['uniform', 'log-uniform']
    low, high = calibration.age_bounds
    prior_density = np.ones(len(ages)) / (high - low) if prior == 'uniform' else 1 / (ages * np.log(high / low))
    unnormalized = np.maximum(0, likelihood) * prior_density
    integral = np.trapezoid(unnormalized, ages)
    assert integral > 0 and np.isfinite(integral)
    density = unnormalized / integral
    cdf = np.r_[0, np.cumsum(np.diff(ages) * (density[1:] + density[:-1]) / 2)]
    quantiles = np.interp([.025, .16, .5, .84, .975], cdf, ages)
    edge = max(2, int(.05 * (len(ages) - 1)))
    lower_mass = float(np.trapezoid(density[:edge + 1], ages[:edge + 1]))
    upper_mass = float(np.trapezoid(density[-edge - 1:], ages[-edge - 1:]))
    flags = []
    if lower_mass > .15 or upper_mass > .15:
        flags.append('Conditional density accumulates near a calibration boundary; no physical age limit is established')
    if retained_temperature_likelihood < .995:
        flags.append('Temperature likelihood is conditioned on the actual finite calibration support')
    return {'status': 'conditional-density', 'age_gyr': ages, 'density_per_gyr': density, 'cdf': cdf,
            'quantiles025_16_50_84_975_gyr': quantiles, 'prior': prior, 'retained_temperature_likelihood': retained_temperature_likelihood,
            'boundary_mass_lower_upper': [lower_mass, upper_mass], 'flags': flags,
            'qualification': 'Model-conditional fixed-scale density; no external age uncertainty, absolute precision, physical boundary limit or universal field-star applicability'}


def age_from_frequency_measure(calibration, measure, observed_temperature, sigma_temperature=40., harmonic_weights=(.5, .5),
                               age_points=401, temperature_order=40, prior='uniform', period_support=(.5, 60.)):
    assert measure['status'] == 'integrated'
    temperature = temperature_quadrature(calibration, observed_temperature, sigma_temperature, temperature_order)
    if temperature['status'] != 'supported':
        return temperature
    ages = np.linspace(*calibration.age_bounds, age_points)
    mean, width = calibration.sequence(temperature['temperatures'], ages)
    log_frequency = np.log(measure['frequency'])
    mass = measure['scaled_likelihood_log_frequency_mass']
    # Keep every positive floating-point likelihood node: a small search tail
    # can become important when the calibration suppresses a dominant peak.
    useful = mass > 0
    log_frequency = log_frequency[useful]; mass = mass[useful]
    likelihood = np.zeros(len(ages))
    # For u=log(frequency), P=m*exp(-u): p(P|t,T)dP is a normal
    # density in log(P) times du. The period search prior is NOT included.
    for multiplicity, mixture_weight in enumerate(harmonic_weights, 1):
        log_period = np.log(multiplicity) - log_frequency
        mask = (log_period >= np.log(period_support[0])) & (log_period <= np.log(period_support[1]))
        for index, temperature_weight in enumerate(temperature['weights']):
            location = mean[index]; scale = width[index]
            truncation = ndtr((np.log(period_support[1]) - location) / scale) - ndtr((np.log(period_support[0]) - location) / scale)
            kernel = np.exp(-.5 * ((log_period[mask, None] - location[None, :]) / scale[None, :]) ** 2) / (scale[None, :] * np.sqrt(2 * np.pi)) / truncation[None, :]
            likelihood += mixture_weight * temperature_weight * (mass[mask] @ kernel)
    return normalize_age(calibration, ages, likelihood, prior, temperature['retained_temperature_likelihood'])


def age_from_reported_log_period(calibration, log_period, sigma_log_period, observed_temperature,
                                 sigma_temperature=40., age_points=401, temperature_order=40, prior='uniform'):
    """Independent analytic Gaussian convolution for verification/conditional catalog probes."""
    temperature = temperature_quadrature(calibration, observed_temperature, sigma_temperature, temperature_order)
    if temperature['status'] != 'supported':
        return temperature
    ages = np.linspace(*calibration.age_bounds, age_points)
    mean, width = calibration.sequence(temperature['temperatures'], ages)
    total_width = np.sqrt(width ** 2 + sigma_log_period ** 2)
    kernel = np.exp(-.5 * ((log_period - mean) / total_width) ** 2) / (total_width * np.sqrt(2 * np.pi))
    likelihood = temperature['weights'] @ kernel
    return normalize_age(calibration, ages, likelihood, prior, temperature['retained_temperature_likelihood'])