"""Conditional proper-Gaussian harmonic inference; no source calibration posterior assumed."""
import numpy as np
from scipy.optimize import minimize_scalar


def harmonic_design(phase):
    angle = 2*np.pi*(phase - .5)
    return np.array([np.ones(len(phase)), np.cos(angle), np.sin(angle), np.cos(2*angle), np.sin(2*angle)]).T


def conditional_posterior(times, values, errors, period, epoch, prior_sd=.02, quality_mask=None):
    phase = ((times - epoch)/period) % 1
    selected = (np.isfinite(values) & np.isfinite(errors) & (errors > 0)
                & (np.minimum(phase, 1-phase) > .045) & (np.abs(phase-.5) > .045))
    if quality_mask is not None:
        selected &= ~np.asarray(quality_mask, dtype=bool)
    design = harmonic_design(phase[selected])
    whitened = design/errors[selected, None]
    information = np.einsum('ni,nj->ij', whitened, whitened) + np.eye(5)/prior_sd**2
    covariance = np.linalg.inv(information)
    mean = np.linalg.solve(information, np.einsum('ni,n->i', whitened, values[selected]/errors[selected]))
    residual = values[selected] - np.einsum('ni,i->n', design, mean)
    return dict(phase=phase, selected=selected, design=design, mean=mean, covariance=covariance, residual=residual)


def phase_summary(coefficients):
    def flux(angle):
        return (coefficients[0] + coefficients[1]*np.cos(angle) + coefficients[2]*np.sin(angle)
                + coefficients[3]*np.cos(2*angle) + coefficients[4]*np.sin(2*angle))
    angles = np.linspace(-np.pi, np.pi, 257)
    grid = flux(angles)
    maximum = np.argmax(grid)
    minimum = np.argmin(grid)
    step = 2*np.pi/256
    peak = minimize_scalar(lambda value: -flux(value), bounds=(angles[maximum]-step, angles[maximum]+step), method='bounded', options={'xatol': 1e-12})
    trough = minimize_scalar(flux, bounds=(angles[minimum]-step, angles[minimum]+step), method='bounded', options={'xatol': 1e-12})
    return [-float(peak.x)*180/np.pi, float(-peak.fun*1e6), float(trough.fun*1e6)]


def interval(values):
    return [float(value) for value in np.percentile(values, [16, 50, 84])]


def wild_block_sensitivity(posterior, times, errors, seed=2121):
    indices = np.flatnonzero(posterior['selected'])
    # Never carry a residual block over a rejected integration or a temporal gap.
    blocks = np.zeros(len(indices), dtype=int)
    block, run = 0, 0
    for index in range(1, len(indices)):
        if indices[index] != indices[index-1]+1 or (times[indices[index]]-times[indices[index-1]])*86400 > 16 or run >= 63:
            block += 1
            run = 0
        else:
            run += 1
        blocks[index] = block
    design = posterior['design']
    selected_errors = errors[posterior['selected']]
    influence = np.einsum('ij,nj->in', posterior['covariance'], design/selected_errors[:, None]**2)
    rng = np.random.default_rng(seed)
    peaks = []
    for _ in range(512):
        signs = rng.choice([-1., 1.], size=block+1)
        change = np.einsum('in,n->i', influence, posterior['residual']*signs[blocks])
        peaks.append(phase_summary(posterior['mean']+change)[0])
    return dict(kind='64-integration wild-block residual sensitivity, not posterior draws',
                seed=seed, draws=512, blocks=block+1,
                phase_peak_before_secondary_degrees_16_50_84=interval(peaks))