"""F5-dev-01: Cobb-Douglas consumer. Closed-form Marshallian demands + independent numerical validation.
Run from workspace root: ~/wp29/venv/bin/python out/analysis.py
Reads only data/parameters.csv; writes out/results.json. Fully deterministic (no randomness; random_seed = null).
"""
import json, math, os
import numpy as np
import pandas as pd
from scipy.optimize import minimize_scalar, minimize

DATA = os.path.join("data", "parameters.csv")
OUT = os.path.join("out", "results.json")

raw = pd.read_csv(DATA)
n_raw = len(raw)
issues = []
if raw["parameter"].duplicated().any():
    issues.append({"issue": "duplicate parameter names", "handling": "first occurrence used"})
par = raw.drop_duplicates("parameter").set_index("parameter")["value"].astype(float)
needed = ["alpha", "p1", "p2", "m"]
missing = [k for k in needed if k not in par.index or not np.isfinite(par[k])]
if missing:
    raise SystemExit(f"missing/non-finite parameters: {missing}")
alpha, p1, p2, m = (float(par[k]) for k in needed)
if not (0 < alpha < 1 and p1 > 0 and p2 > 0 and m > 0):
    raise SystemExit("parameters outside the range where the interior closed form is the unique optimum")
n_analysed = int(sum(raw["parameter"].isin(needed)))


def u(x1, x2):
    return alpha * np.log(x1) + (1 - alpha) * np.log(x2)


# ---------- Closed form (Lagrangian FOCs: alpha/x1 = lam p1, (1-alpha)/x2 = lam p2, budget binds => lam = 1/m) ----------
x1 = alpha * m / p1
x2 = (1 - alpha) * m / p2
v = float(u(x1, x2))
resid = p1 * x1 + p2 * x2 - m
mrs = (alpha / x1) / ((1 - alpha) / x2)
price_ratio = p1 / p2
share1 = p1 * x1 / m

# ---------- Independent numerical validation (does NOT use the closed form) ----------
# Budget line (monotone utility => budget binds): x1 in (0, m/p1), x2 = (m - p1 x1)/p2.
def neg_u_line(x1v):
    return -u(x1v, (m - p1 * x1v) / p2)

ub = m / p1
# (a) dense grid search on the budget line
grid = np.linspace(ub * 1e-6, ub * (1 - 1e-6), 2_000_001)
vals = -neg_u_line(grid)
i = int(np.argmax(vals))
x1_grid, v_grid = float(grid[i]), float(vals[i])
# (b) bounded scalar optimiser (Brent) on the budget line
res_b = minimize_scalar(neg_u_line, bounds=(ub * 1e-9, ub * (1 - 1e-9)), method="bounded",
                        options={"xatol": 1e-12})
x1_brent, v_brent = float(res_b.x), float(-res_b.fun)
# (c) 2-D SLSQP with the inequality budget constraint (checks that the budget binds without imposing it)
res_s = minimize(lambda z: -u(z[0], z[1]), x0=np.array([1.0, 1.0]), method="SLSQP",
                 bounds=[(1e-9, None), (1e-9, None)],
                 constraints=[{"type": "ineq", "fun": lambda z: m - p1 * z[0] - p2 * z[1]}],
                 options={"ftol": 1e-14, "maxiter": 1000})
v_slsqp = float(-res_s.fun)

v_num = max(v_grid, v_brent, v_slsqp)  # best utility found numerically
x1_num = x1_brent
x2_num = (m - p1 * x1_num) / p2

results = {
    "task_id": "F5-dev-01",
    "x1_demand": x1,
    "x2_demand": x2,
    "indirect_utility": v,
    "budget_identity_residual": resid,
    "mrs_at_optimum": mrs,
    "price_ratio": price_ratio,
    "expenditure_share_good1": share1,
    "indirect_utility_numerical": v_num,
    "feasibility": "complete",
    "feasibility_reason": None,
    "data_issues": issues,
    "n_input_rows": {"parameters.csv": n_raw},
    "n_analysed": n_analysed,
    "random_seed": None,
    # diagnostics of the numerical validation
    "numerical_x1_brent": x1_num,
    "numerical_x2_brent": x2_num,
    "numerical_x1_grid": x1_grid,
    "numerical_utility_grid": v_grid,
    "numerical_utility_brent": v_brent,
    "numerical_utility_slsqp": v_slsqp,
    "numerical_x_slsqp": [float(res_s.x[0]), float(res_s.x[1])],
    "numerical_slsqp_budget_slack": float(m - p1 * res_s.x[0] - p2 * res_s.x[1]),
    "utility_gap_numerical_minus_closed_form": v_num - v,
    "max_abs_demand_gap_x1": max(abs(x1_num - x1), abs(x1_grid - x1), abs(float(res_s.x[0]) - x1)),
    "max_abs_demand_gap_x2": max(abs(x2_num - x2), abs(float(res_s.x[1]) - x2)),
}
with open(OUT, "w") as f:
    json.dump(results, f, indent=2)
print(f"x1*={x1:.6f} x2*={x2:.6f} v={v:.10f} v_num={v_num:.10f} residual={resid:.3e} MRS={mrs:.6f} p1/p2={price_ratio:.6f}")