"""Finite conditional least-squares profiles; no Bayesian/particle-model verdict."""
import math
import numpy as np
from scipy.optimize import least_squares
from rotation_model import forward,observation_coordinate,log_likelihood

BOUNDS=([.05,3,math.log10(.05)],[1.5,10,math.log10(200)])
STARTS=[[.2,6,1.3],[.5,7,.7],[.8,7.5,.2],[1.2,6.5,1],[.1,8.5,-.5],[1.4,5,2]]

def fit(record,kind,distance_ratio=1.,inclination_degrees=None,fixed_mass_to_light=None):
    observed,error=observation_coordinate(record)
    def unpack(parameters):
        full=parameters if fixed_mass_to_light is None else [fixed_mass_to_light,*parameters]
        return float(full[0]),10**float(full[1]),10**float(full[2])
    def residual(parameters):
        prediction=forward(record,*unpack(parameters),kind,distance_ratio,inclination_degrees)
        return (prediction['referenceProjectedModelKms']-observed)/error
    attempts=[];solutions=[]
    for start in STARTS:
        candidate=start if fixed_mass_to_light is None else start[1:]
        bounds=BOUNDS if fixed_mass_to_light is None else (BOUNDS[0][1:],BOUNDS[1][1:])
        solution=least_squares(residual,candidate,bounds=bounds,xtol=1e-10,ftol=1e-10,gtol=1e-10,max_nfev=1500)
        attempts.append({'success':bool(solution.success),'status':int(solution.status),'cost':float(solution.cost),'nfev':int(solution.nfev)})
        if solution.success:solutions.append(solution)
    if not solutions:raise RuntimeError('No completed conditional optimization')
    best=min(solutions,key=lambda value:value.cost);mass_to_light,density,scale=unpack(best.x)
    prediction=forward(record,mass_to_light,density,scale,kind,distance_ratio,inclination_degrees)
    lower=np.asarray(BOUNDS[0] if fixed_mass_to_light is None else BOUNDS[0][1:]);upper=np.asarray(BOUNDS[1] if fixed_mass_to_light is None else BOUNDS[1][1:])
    near=bool(np.any(np.minimum(best.x-lower,upper-best.x)<1e-4))
    return {'haloKind':kind,'massToLight':mass_to_light,'haloDensitySolarPerKpc3':density,'haloScaleKpc':scale,
            'distanceRatio':distance_ratio,'inclinationDegrees':inclination_degrees or record['metadata']['inclinationDegrees'],
            'chiSquared':float(2*best.cost),'workingReducedChiSquared':float(2*best.cost/(len(observed)-len(best.x))),
            'logWorkingGaussianLikelihood':log_likelihood(record,prediction),'boundaryBest':near,
            'jacobianSingularValues':np.linalg.svd(best.jac,compute_uv=False).tolist(),'attempts':attempts,
            'predictedRecord':{key:values.tolist() for key,values in prediction.items()},
            'interpretation':'Finite conditional working-likelihood optimum; no calibrated coverage, evidence or theory/particle verdict'}

def analyze_record(record,profiles=True,nuisance=True):
    result={'fixedDataLabel':record['label'],'fits':{},'profileMassToLight':{},'sharedNuisanceSensitivity':{}}
    for kind in ['nfw','core']:
        result['fits'][kind]=fit(record,kind)
        if profiles:
            result['profileMassToLight'][kind]=[fit(record,kind,fixed_mass_to_light=float(value)) for value in np.arange(.1,1.01,.1)]
        if nuisance:
            metadata=record['metadata'];distance=metadata['distanceMpc'];distance_sigma=metadata['distanceErrorMpc'];inclination=metadata['inclinationDegrees'];inclination_sigma=metadata['inclinationErrorDegrees']
            result['sharedNuisanceSensitivity'][kind]=[fit(record,kind,distance_ratio=(distance+d_step*distance_sigma)/distance,
                inclination_degrees=inclination+i_step*inclination_sigma) for d_step in [-1,0,1] for i_step in [-1,0,1]]
    return result