"""Resolve all coarse likelihood peaks, then integrate a positive adaptive measure.

The numerical work log and refinement comparison are not a formal interval
certificate. Inadequate work remains an explicit budget-limited result.
"""
import numpy as np
from numpy.polynomial.legendre import leggauss
from scipy.optimize import minimize_scalar


class WorkLimit(Exception):
    pass


def frequency_measure(likelihood, minimum_period=.5, maximum_period=60., samples_per_cycle=16,
                      tolerance=2e-6, maximum_evaluations=80000, quadrature_order=8):
    minimum_frequency = 1 / maximum_period
    maximum_frequency = 2 / minimum_period
    duration = np.ptp(likelihood.times)
    count = int(np.ceil((maximum_frequency - minimum_frequency) * samples_per_cycle * duration)) + 1
    def evaluate(values):
        if likelihood.evaluations + np.size(values) > maximum_evaluations:
            raise WorkLimit('Actual harmonic likelihood work budget reached')
        return likelihood.evaluate(values)
    coarse_frequency = np.linspace(minimum_frequency, maximum_frequency, count)
    try:
        coarse_log = evaluate(coarse_frequency)
        local = np.flatnonzero((coarse_log[1:-1] >= coarse_log[:-2]) & (coarse_log[1:-1] >= coarse_log[2:])) + 1
        peaks = []
        step = coarse_frequency[1] - coarse_frequency[0]
        for index in local:
            fit = minimize_scalar(lambda value: -evaluate([value])[0], bounds=(coarse_frequency[index - 1], coarse_frequency[index + 1]), method='bounded', options={'xatol': 2e-12, 'maxiter': 80})
            if not fit.success:
                raise WorkLimit('One-dimensional peak refinement did not meet its stopping rule')
            center = float(fit.x); height = float(-fit.fun)
            probe = min(step / 8, center - minimum_frequency, maximum_frequency - center)
            shoulders = evaluate([center - probe, center + probe])
            curvature = (shoulders[0] - 2 * height + shoulders[1]) / probe ** 2
            width = np.sqrt(1 / max(-curvature, 1 / (4 * step ** 2)))
            peaks.append({'frequency': center, 'log_evidence': height, 'width': float(width), 'optimizer_success': bool(fit.success)})
        highest = max([float(max(coarse_log))] + [peak['log_evidence'] for peak in peaks])
        retained = [peak for peak in peaks if peak['log_evidence'] >= highest - 40]
        boundaries = [minimum_frequency, maximum_frequency, 2 / maximum_period, 1 / minimum_period]
        for peak in retained:
            boundaries.extend([max(minimum_frequency, peak['frequency'] - 9 * peak['width']), peak['frequency'], min(maximum_frequency, peak['frequency'] + 9 * peak['width'])])
        # Smooth age kernels are integrated accurately even between resolved peaks.
        boundaries.extend(np.exp(np.linspace(np.log(minimum_frequency), np.log(maximum_frequency), int(np.ceil(np.log(maximum_frequency / minimum_frequency) / .16)) + 1)))
        boundaries = sorted(set(float(value) for value in boundaries if minimum_frequency <= value <= maximum_frequency))
        coarse_mass = np.trapezoid(np.exp(coarse_log - highest) / coarse_frequency, coarse_frequency)
        laplace_mass = sum(np.exp(peak['log_evidence'] - highest) * np.sqrt(2 * np.pi) * peak['width'] / peak['frequency'] for peak in retained)
        estimated_mass = max(coarse_mass, laplace_mass, 1e-300)
        local_tolerance = tolerance * estimated_mass / (len(boundaries) - 1)
        rules = [leggauss(quadrature_order), leggauss(2 * quadrature_order)]
        nodes = []; masses = []; largest_error = 0.; leaves = 0
        def rule(left, right, index):
            abscissas, weights = rules[index]
            frequency = (left + right) / 2 + (right - left) / 2 * abscissas
            log_evidence = evaluate(frequency)
            mass = weights * (right - left) / 2 / frequency * np.exp(log_evidence - highest)
            scaled_log = (np.log(frequency) - np.log(minimum_frequency)) / np.log(maximum_frequency / minimum_frequency)
            moments = np.array([np.sum(mass), np.sum(mass / frequency / maximum_period), np.sum(mass * scaled_log), np.sum(mass * scaled_log ** 2)])
            return frequency, mass, moments
        def integrate(left, right, depth=0):
            nonlocal largest_error, leaves
            low = rule(left, right, 0); high = rule(left, right, 1)
            error = float(max(abs(high[2] - low[2])))
            if error > local_tolerance and depth < 16:
                middle = (left + right) / 2
                integrate(left, middle, depth + 1); integrate(middle, right, depth + 1)
            else:
                if error > local_tolerance:
                    raise WorkLimit('Adaptive measure depth limit before target')
                largest_error = max(largest_error, error)
                leaves += 1; nodes.extend(high[0]); masses.extend(high[1])
        for left, right in zip(boundaries[:-1], boundaries[1:]):
            integrate(left, right)
        frequency = np.array(nodes); mass = np.array(masses)
        order = np.argsort(frequency)
        best_peak = max(peaks, key=lambda value: value['log_evidence']) if peaks else {'frequency': float(coarse_frequency[np.argmax(coarse_log)])}
        if likelihood.evaluations + 1 > maximum_evaluations:
            raise WorkLimit('Work limit before final diagnostic')
        diagnostic = likelihood.evaluate([best_peak['frequency']], diagnostics=True)[0]
        return {'status': 'integrated', 'frequency': frequency[order], 'scaled_likelihood_log_frequency_mass': mass[order], 'log_scale': highest,
                'peaks': sorted(retained, key=lambda value: -value['log_evidence']), 'coarse_samples': count, 'all_local_peaks_refined': len(peaks),
                'actual_evidence_evaluations': likelihood.evaluations, 'adaptive_leaves': leaves, 'largest_accepted_moment_error': largest_error,
                'target_relative_moment_tolerance': tolerance, 'best_fit_diagnostic': diagnostic,
                'scope': 'Positive adaptive floating-point quadrature with independent refinement required; no formal global certificate'}
    except WorkLimit as error:
        return {'status': 'budget-limited', 'reason': str(error), 'actual_evidence_evaluations': likelihood.evaluations}


def period_posterior(measure, likelihood, minimum_period=.5, maximum_period=60., harmonic_weights=(.5, .5), prior='log-uniform'):
    assert measure['status'] == 'integrated'
    assert prior in ['log-uniform', 'uniform']
    assert all(value >= 0 for value in harmonic_weights) and abs(sum(harmonic_weights) - 1) < 1e-12
    frequency = measure['frequency']; base_mass = measure['scaled_likelihood_log_frequency_mass']
    periods = []; masses = []; labels = []
    for multiplicity, model_weight in enumerate(harmonic_weights, 1):
        period = multiplicity / frequency
        chosen = (period >= minimum_period) & (period <= maximum_period)
        prior_factor = np.full(np.sum(chosen), 1 / np.log(maximum_period / minimum_period)) if prior == 'log-uniform' else period[chosen] / (maximum_period - minimum_period)
        periods.extend(period[chosen]); masses.extend(base_mass[chosen] * model_weight * prior_factor); labels.extend([multiplicity] * np.sum(chosen))
    periods = np.array(periods); masses = np.array(masses); labels = np.array(labels)
    integral = np.sum(masses)
    assert integral > 0
    masses /= integral
    order = np.argsort(periods); cumulative = np.cumsum(masses[order])
    credible = np.interp([.05, .16, .5, .84, .95], cumulative, periods[order])
    log_bayes_factor = measure['log_scale'] + np.log(integral) - likelihood.constant_log_evidence()
    reduced_chi = measure['best_fit_diagnostic']['residual_chi_square'] / max(1, len(likelihood.times) - len(likelihood.prior_variances))
    flags = []
    if log_bayes_factor < np.log(10): flags.append('Stationary rotation evidence is weak relative to the proper constant-plus-trend model')
    if reduced_chi > 3: flags.append('Stationary harmonic residuals fail the declared noise model; pattern evolution or another signal is unresolved')
    if np.sum(masses[periods <= np.ptp(likelihood.times) / 2]) < .5 - 1e-10: flags.append('Most clock probability spans fewer than two observed rotations')
    return {'period_days': periods, 'probability_mass': masses, 'multiplicity': labels, 'quantiles05_16_50_84_95_days': credible,
            'log_bayes_factor_vs_constant_trend': log_bayes_factor, 'best_reduced_chi_square': reduced_chi,
            'harmonic_model_probabilities': [float(np.sum(masses[labels == value])) for value in [1, 2]], 'flags': flags,
            'prior': prior, 'harmonic_model_prior_weights': list(harmonic_weights), 'qualification': 'Surface-pattern clock with declared one/two-longitude ambiguity, not measured equatorial spin; Bayes factor is not a FAP'}