"""Conditional nonlinear negative-injection flux fit with full re-estimated KL bases.

The already observed source photons retain their noise when a deterministic
model is subtracted. The adopted covariance stays fixed through this local fit;
its log determinant is included, not optimized away by inflating noise.
"""
import numpy as np
from imaging_reduction import reduce_cube, measurements


def fit_negative_injection(data, templates, angles, wavelengths, weights, covariance, starting_contrasts, algorithm="klip", maximum_iterations=8):
    inverse = np.linalg.inv(covariance)
    logdet = float(np.linalg.slogdet(covariance)[1])
    contrasts = np.asarray(starting_contrasts).copy()
    step = 2e-5

    def residual(candidate):
        corrected = data - 3e7 * candidate[:, None, None, None] * templates
        return measurements(reduce_cube(corrected, angles, wavelengths, algorithm), weights)

    def objective(candidate):
        value = residual(candidate)
        return float(value @ inverse @ value + logdet + len(value) * np.log(2 * np.pi))

    history = []
    status = "iteration-limit"
    for iteration in range(maximum_iterations):
        value = residual(contrasts)
        columns = []
        for channel in range(len(contrasts)):
            direction = np.zeros(len(contrasts)); direction[channel] = step
            columns.append((residual(contrasts + direction) - residual(contrasts - direction)) / (2 * step))
        jacobian = np.asarray(columns).T
        normal = jacobian.T @ inverse @ jacobian
        eigenvalues = np.linalg.eigvalsh(normal)
        if eigenvalues[0] <= eigenvalues[-1] * 1e-8:
            status = "rank-deficient"; break
        update = np.linalg.solve(normal, -jacobian.T @ inverse @ value)
        previous = objective(contrasts)
        scale = 1.0
        while scale >= 1 / 64 and objective(contrasts + scale * update) > previous:
            scale /= 2
        history.append({"iteration": iteration, "negative2loglikelihood": previous, "step_norm": float(np.linalg.norm(update)), "accepted_scale": scale})
        if scale < 1 / 64:
            status = "line-search-limit"; break
        contrasts += scale * update
        if np.max(np.abs(scale * update)) < 1e-9:
            status = "converged"; break
    final = residual(contrasts)
    return {"status": status, "contrast": contrasts.tolist(), "negative2loglikelihood": objective(contrasts), "chi_square": float(final @ inverse @ final), "history": history,
            "interpretation": "Conditional nonlinear negative-injection flux fit at the adopted position. Every evaluation rebuilds the subtraction basis and the target-minus-comparison measurement. Covariance is fixed; this is not a posterior or a calibrated confidence interval."}