"""Exact Gaussian mean-CVaR problem in integer shares as a MISOCP with Gurobi (validation of the cutting planes)."""
import sys, time
import numpy as np, gurobipy as gp
from scipy.stats import norm
from quantcore.config import get_settings
from quantcore.data.market_data import MarketDataService
from quantcore.optimization.service import OptimizationService, OptimizationRequest
from quantcore.optimization.problem import OptimizationConstraints
from quantcore.models.covariance import CovarianceModel

N_ASSETS = int(sys.argv[1]) if len(sys.argv) > 1 else 30
BETA = float(sys.argv[2]) if len(sys.argv) > 2 else 0.95
W, H = 1_000_000.0, 5
LAMBDA = 0.05 * W
settings = get_settings(); md = MarketDataService(settings); svc = OptimizationService(settings, md)
req = OptimizationRequest(max_assets=N_ASSETS, scenario_model="gbm", solver="highs", cash=W)
symbols = svc.candidate_symbols(req)
mu_ann, _, _ = svc.expected_returns(req, symbols)
r = md.returns(symbols, currency="base", freq="daily").dropna()
cov_h = CovarianceModel.from_dict(req.covariance).estimate(r).values * H
mu_h = (mu_ann / 252).reindex(symbols).values * H
prices = md.latest_prices(symbols, currency="base").reindex(symbols).values
sectors = md.sectors(symbols).fillna("Unknown").tolist()
N = len(symbols); k = norm.pdf(norm.ppf(BETA)) / (1 - BETA)
L = np.linalg.cholesky(cov_h)

for integer in (True, False):
    t0 = time.time()
    m = gp.Model(); m.Params.OutputFlag = 0; m.Params.MIPGap = 1e-9; m.Params.MIPGapAbs = 1e-6
    ub = np.floor(0.10 * W / prices)
    n = m.addMVar(N, vtype="I" if integer else "C", lb=0, ub=ub)
    t = m.addVar(lb=0)
    y = m.addMVar(N, lb=-gp.GRB.INFINITY)            # y = L' v  ->  ||y|| = sqrt(v' Sigma v)
    v = prices * n
    m.addConstr(y == L.T @ v)
    m.addConstr(y @ y <= t * t)                        # second-order cone
    m.addConstr(k * t <= LAMBDA + mu_h @ v)            # Gaussian CVaR <= Lambda
    m.addConstr(prices @ n <= W)
    for sec in sorted(set(sectors)):
        mask = np.array([s == sec for s in sectors], dtype=float)
        m.addConstr((prices * mask) @ n <= OptimizationConstraints().max_sector_weight * W)
    m.setObjective(mu_h @ v, gp.GRB.MAXIMIZE)
    m.optimize()
    nn = n.X; vv = prices * nn; sd = np.sqrt(vv @ cov_h @ vv)
    print(f"{'integer' if integer else 'continuous':10s} N={N} beta={BETA}: status {m.Status}, E[P&L]={mu_h @ vv:.2f} CHF, "
          f"CVaR={-mu_h @ vv + k * sd:.2f}, positions={(nn > 0.5).sum()}, {time.time()-t0:.1f}s")