1"""CES aggregators with analytic first derivatives.23ces(weights, x, rho) = (sum_i w_i x_i^rho)^(1/rho) for rho != 0, and the Cobb-Douglas limit4prod_i x_i^{w_i} (weights normalized to one) for rho == 0. Zero-weight inputs are ignored, which lets the5same code run the pure Restrepo benchmark (no structures, land or energy) and the full model.6"""7from __future__ import annotations89import numpy as np1011_FLOOR = 1e-300121314def ces(weights, x, rho: float):15 w = np.asarray(weights, dtype=float)16 x = np.asarray(x, dtype=float)17 m = w > 018 if not m.any():19 return 0.020 w, x = w[m], np.maximum(x[m], _FLOOR)21 if abs(rho) < 1e-12:22 ws = w / w.sum()23 return float(np.exp(np.sum(ws * np.log(x))))24 return float(np.sum(w * x ** rho) ** (1.0 / rho))252627def ces_grad(weights, x, rho: float):28 """Vector of partial derivatives d ces / d x_i (zero for zero-weight inputs)."""29 w = np.asarray(weights, dtype=float)30 x = np.asarray(x, dtype=float)31 g = np.zeros_like(x)32 m = w > 033 if not m.any():34 return g35 val = ces(w, x, rho)36 xm = np.maximum(x[m], _FLOOR)37 if abs(rho) < 1e-12:38 ws = w[m] / w[m].sum()39 g[m] = ws * val / xm40 else:41 g[m] = w[m] * xm ** (rho - 1.0) * val ** (1.0 - rho)42 return g434445def rho_of_sigma(sigma: float) -> float:46 """rho = (sigma-1)/sigma; sigma = 1 -> 0 (Cobb-Douglas); sigma = inf -> 1 (linear)."""47 if np.isinf(sigma):48 return 1.049 return (sigma - 1.0) / sigma50