"""RN-05 v1.3: population sensitivity of an invented loss mixture.

Run: python tail_sensitivity_v1_3.py > tail_sensitivity_results_v1_3.json
Requires NumPy and SciPy. Reads no files or network resources.
All units and probabilities are synthetic. This is not a backtest.
"""
import json
import platform
import numpy as np
import scipy
from scipy.integrate import quad
from scipy.optimize import brentq, minimize, approx_fprime
from scipy.stats import norm

ALPHA, CAP = .99, .80
COMMON = .20
RESIDUAL = np.array([.35, .95, .65])
JUMP = np.array([4., 0., 1.5])


def components(w, p):
    sd = np.sqrt(COMMON**2 * w.sum()**2 + np.sum((RESIDUAL*w)**2))
    multipliers = np.array([-p, 1-p])
    return np.array([1-p, p]), multipliers*(JUMP@w), sd, multipliers


def ru(x, p):
    w, eta = x[:3], x[3]
    probs, means, sd, multipliers = components(w, p)
    z = (eta-means)/sd
    sf, pdf = norm.sf(z), norm.pdf(z)
    value = eta + probs@(sd*pdf+(means-eta)*sf)/(1-ALPHA)
    dsd = (COMMON**2*w.sum() + RESIDUAL**2*w)/sd
    dw = ((probs*sf)@multipliers*JUMP + (probs@pdf)*dsd)/(1-ALPHA)
    deta = 1 - probs@sf/(1-ALPHA)
    return float(value), np.r_[dw, deta]


def population_risk(w, p):
    probs, means, sd, _ = components(w, p)
    eta = brentq(lambda t: probs@norm.cdf((t-means)/sd)-ALPHA,
                 min(means)-12*sd, max(means)+12*sd, xtol=1e-13)
    z = (eta-means)/sd
    es = probs@(means*norm.sf(z)+sd*norm.pdf(z))/(1-ALPHA)
    # Independent numerical integration of loss times the mixture density.
    integrated = sum(prob*quad(lambda t: t*norm.pdf((t-mean)/sd)/sd,
                              eta, np.inf, epsabs=1e-11)[0]
                     for prob, mean in zip(probs, means))/(1-ALPHA)
    assert abs(es-integrated) < 1e-8
    variance = sd**2 + p*(1-p)*(JUMP@w)**2
    return float(eta), float(es), float(np.sqrt(variance)), float(abs(es-integrated))


def kkt(w, gradient):
    interior = (w>1e-7)&(w<CAP-1e-7)
    level = gradient[interior].mean()
    errors = np.where(w<=1e-7, np.maximum(level-gradient, 0),
                      np.where(w>=CAP-1e-7, np.maximum(gradient-level, 0),
                               np.abs(gradient-level)))
    error = float(errors.max())
    assert abs(w.sum()-1)<1e-9 and w.min()>-1e-9 and w.max()<CAP+1e-9
    assert error<2e-6, (w, gradient, error)
    return error


def solve(p, objective):
    sigma = COMMON**2*np.ones((3,3)) + np.diag(RESIDUAL**2) + p*(1-p)*np.outer(JUMP,JUMP)
    if objective=='minimum_variance':
        function = lambda x: (float(x@sigma@x), 2*sigma@x)
        bounds = [(0,CAP)]*3
        starts = [np.ones(3)/3, np.array([.1,.6,.3])]
    else:
        function = lambda x: ru(x,p)
        bounds = [(0,CAP)]*3+[(None,None)]
        starts = [np.r_[w,population_risk(w,p)[0]]
                  for w in [np.ones(3)/3, np.array([.1,.6,.3])]]
    check_point = starts[0]
    numerical = approx_fprime(check_point, lambda x:function(x)[0], 1e-7)
    gradient_error = float(np.max(np.abs(numerical-function(check_point)[1])))
    assert gradient_error<2e-5
    solutions = [minimize(function, x, jac=True, method='SLSQP', bounds=bounds,
                constraints={'type':'eq','fun':lambda x:x[:3].sum()-1,
                             'jac':lambda x:np.r_[np.ones(3),np.zeros(len(x)-3)]},
                options={'ftol':1e-13,'maxiter':1000}) for x in starts]
    assert all(r.success for r in solutions), [r.message for r in solutions]
    assert abs(solutions[0].fun-solutions[1].fun)<1e-9
    r = solutions[0]
    w = r.x[:3]
    gradient = function(r.x)[1]
    residual = kkt(w,gradient[:3])
    eta,es,sd,integration_error = population_risk(w,p)
    if objective=='minimum_population_es':
        assert abs(r.fun-es)<1e-9 and abs(r.x[3]-eta)<2e-6
        assert abs(gradient[3])<2e-6
    return {'event_probability':p,'objective':objective,'weights':w.tolist(),
            'population_sd':sd,'population_es_99':es,'population_var_99':eta,
            'validation':{'budget_residual':float(abs(w.sum()-1)),
              'kkt_residual':residual,'gradient_error':gradient_error,
              'independent_start_objective_difference':float(abs(solutions[0].fun-solutions[1].fun)),
              'tail_quadrature_error':integration_error,'solver_success':True,
              'cap_nonbinding':bool(w.max()<CAP-1e-7)}}


def reproduce():
    rows = [solve(p,objective) for p in [.01,.04,.10]
            for objective in ['minimum_variance','minimum_population_es']]
    for a,b in zip(rows[::2],rows[1::2]):
        assert b['population_es_99']<=a['population_es_99']+1e-9
        assert a['population_sd']<=b['population_sd']+1e-9
    return {'note':'PSIM-RN-05','version':'1.3','last_revised':'2026-10-01',
            'method':'Exact population two-normal mixture; convex RU expectation minimized with SLSQP; no scenario sample.',
            'parameters':{'alpha':ALPHA,'cap':CAP,'common_scale':COMMON,
              'residual_scales':RESIDUAL.tolist(),'event_loadings':JUMP.tolist(),
              'centering':'I-p is re-centred at each event probability'},
            'rows':rows,'all_checks_passed':True,
            'runtime':{'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__}}


if __name__=='__main__':
    print(json.dumps(reproduce(),indent=2,allow_nan=False))