"""MAXSER-PR with a full-rank zero-investment common-prediction-shift basis.

Annual raw-feature LightGBMs use only matured past labels. Five delete-block
fits affect the portfolio basis, while the main prediction remains unchanged.
POET K=10; no-intercept elastic net; alpha=0,.01,...,.10; both CV selections
maximize mean-fold squared Sharpe. Raw weights are returned without gross scaling.
"""
from __future__ import annotations
from functools import lru_cache
from typing import Callable
import math
import warnings
import numpy as np
import pandas as pd
import lightgbm as lgb
from numpy.typing import ArrayLike, NDArray
from scipy.special import roots_jacobi
from numba import njit
from joblib import Parallel, delayed, parallel_config

TARGET_DAILY_SIGMA = 0.04 / np.sqrt(252.0)
_DOUBLE_EPS = np.finfo(np.float64).eps
_SQRT_DOUBLE_EPS = np.sqrt(_DOUBLE_EPS)
FloatArray = NDArray[np.float64]
IntArray = NDArray[np.int64]
ALPHA_GRID = np.arange(11, dtype=np.float64) / 100.0
N_LAMBDA = 100
LAMBDA_MIN_RATIO = 0.01
ALPHA_ZERO_REFERENCE = 0.01
CV_FOLDS = 10
CV_SEED = 123
_BLAS_LIMITER = None


def _as_finite_float64(value: object, *, ndim: int, name: str) -> np.ndarray:
    array = np.asarray(value, dtype=np.float64)
    if array.ndim != ndim or not np.isfinite(array).all():
        raise ValueError(f"Invalid {name}")
    return array

def spatial_factors(Y: object, rmax: int) -> np.ndarray:
    matrix = _as_finite_float64(Y, ndim=2, name="spatial-factor input")
    nr, nc = matrix.shape
    rmax = int(rmax)
    if nr < 3 or nc < 2 or rmax < 1 or (rmax + 1 > min(nr - 1, nc)):
        raise ValueError("Invalid rmax")
    # Reproduce the spatial-factor step used by the baseline POET path.
    shifted = matrix - matrix[0]
    squared_norm = np.einsum("ij,ij->i", shifted, shifted)
    distance2 = squared_norm[:, None] + squared_norm[None, :] - 2.0 * (shifted @ shifted.T)
    np.fill_diagonal(distance2, np.inf)
    bad_i, bad_j = np.where(np.triu(distance2 <= 0.0, k=1))
    for i, j in zip(bad_i.tolist(), bad_j.tolist()):
        exact = float(np.sum((matrix[i] - matrix[j]) ** 2))
        if exact <= 0.0:
            raise ValueError("Duplicate spatial-factor rows")
        distance2[i, j] = distance2[j, i] = exact
    with np.errstate(divide="ignore", invalid="ignore"):
        adjacency = 1.0 / distance2
    np.fill_diagonal(adjacency, 0.0)
    laplacian = np.diag(np.sum(adjacency, axis=1)) - adjacency
    eigenvalues, eigenvectors = np.linalg.eigh((laplacian + laplacian.T) / 2.0)
    order = np.arange(eigenvalues.size - 1, -1, -1)
    eigenvalues, eigenvectors = (eigenvalues[order], eigenvectors[:, order])
    keep = eigenvalues > eigenvalues[0] * 1e-13
    weighted = (
        np.sqrt(np.maximum(eigenvalues[keep], 0.0))[:, None] * eigenvectors[:, keep].T @ shifted
    )
    weighted *= np.sqrt(2.0 / (nr * (nr - 1.0)))
    _, singular_values, right_t = np.linalg.svd(weighted, full_matrices=False)
    values = singular_values[: rmax + 1] ** 2
    rank = int(np.argmax(values[:rmax] / values[1 : rmax + 1])) + 1
    return right_t[:rank].T.copy()

def gamma_path(gamma_design: object, mu: object) -> np.ndarray:
    design = _as_finite_float64(gamma_design, ndim=2, name="gamma design")
    mean = _as_finite_float64(mu, ndim=1, name="synthetic mean")
    nr, nc = design.shape
    if mean.size != nc:
        raise ValueError("Invalid synthetic mean")
    # Remove the estimated common-factor space before measuring idiosyncratic scale.
    vectors = spatial_factors(design, min(8, nc - 1, nr - 2))
    centered = design - mean
    projected = centered - centered @ vectors @ vectors.T
    projected_sq = np.einsum("ij,ij->i", projected, projected)
    if np.any(projected_sq <= 0.0):
        raise ValueError("Degenerate projected residual norm")
    centered_sq = np.einsum("ij,ij->i", centered, centered)
    eta = 1.0 / float(np.mean(centered_sq / projected_sq))
    return np.sqrt(projected_sq) / np.sqrt(eta) / np.sqrt(float(nc))

def poet_covariance(standardized_residual: object, K: int = 3) -> np.ndarray:
    residual = _as_finite_float64(standardized_residual, ndim=2, name="POET input")
    Y = residual.T.copy()
    p, n = Y.shape
    K = int(K)
    if K < 0 or K >= min(p, n):
        raise ValueError("Invalid POET input")
    Y -= np.mean(Y, axis=1, keepdims=True)
    # The low-rank component uses K principal factors; the diagonal keeps residual risk.
    if K > 0:
        gram = Y.T @ Y
        gram = (gram + gram.T) / 2.0
        eigenvalues, eigenvectors = np.linalg.eigh(gram)
        order = np.arange(eigenvalues.size - 1, eigenvalues.size - K - 1, -1)
        F = np.sqrt(float(n)) * eigenvectors[:, order]
        loadings = Y @ F / float(n)
        uhat = Y - loadings @ F.T
        lowrank = loadings @ loadings.T
    else:
        uhat = Y
        lowrank = np.zeros((p, p), dtype=np.float64)
    variance = np.einsum("ij,ij->i", uhat, uhat) / float(n)
    if not np.isfinite(variance).all() or np.any(variance <= 0.0):
        raise ValueError("Invalid POET residual variance")
    covariance = lowrank.copy()
    covariance[np.diag_indices(p)] += variance
    covariance = (covariance + covariance.T) / 2.0
    if not np.isfinite(covariance).all():
        raise ValueError("Nonfinite POET covariance")
    return covariance

def zero_invest_theta(means: object, covariance: object, gamma: object) -> np.ndarray:
    mean_matrix = _as_finite_float64(means, ndim=2, name="theta means")
    sigma = _as_finite_float64(covariance, ndim=2, name="theta covariance")
    scale = _as_finite_float64(gamma, ndim=1, name="theta gamma")
    if (
        sigma.shape != (mean_matrix.shape[1], mean_matrix.shape[1])
        or scale.size != mean_matrix.shape[0]
    ):
        raise ValueError("Invalid theta dimensions")
    try:
        lower = np.linalg.cholesky(sigma)
    except np.linalg.LinAlgError as exc:
        raise ValueError(f"POET covariance is not positive definite: {exc}") from exc
    # Project whitened means away from the budget vector, imposing 1'w = 0.
    whitened = np.linalg.solve(lower, mean_matrix.T)
    budget = np.linalg.solve(lower, np.ones(mean_matrix.shape[1], dtype=np.float64))
    denominator = float(budget @ budget)
    coefficient = budget @ whitened / denominator
    projected = whitened - budget[:, None] * coefficient[None, :]
    theta = np.einsum("ij,ij->j", projected, projected) / scale**2
    if not np.isfinite(theta).all() or np.any(theta < 0.0):
        raise ValueError("Invalid theta path")
    return theta

def _r_numeric_gradient(
    function: Callable[[np.ndarray], float], parameters: np.ndarray
) -> np.ndarray:
    gradient = np.empty_like(parameters)
    for index in range(parameters.size):
        plus = parameters.copy()
        minus = parameters.copy()
        plus[index] += 0.001
        minus[index] -= 0.001
        gradient[index] = (function(plus) - function(minus)) / 0.002
    return gradient

def _r_vmmin(
    function: Callable[[np.ndarray], float], initial: np.ndarray, maxit: int = 100
) -> np.ndarray:
    parameters = np.asarray(initial, dtype=np.float64).copy()
    n = parameters.size
    inverse_hessian = np.eye(n, dtype=np.float64)
    value = float(function(parameters))
    if not np.isfinite(value):
        raise ValueError("Initial AR(1) likelihood is not finite")
    minimum = value
    gradient = _r_numeric_gradient(function, parameters)
    gradient_count = 1
    iteration = 1
    last_reset = gradient_count
    count = 0
    relative_tolerance = np.sqrt(_DOUBLE_EPS)
    while True:
        if last_reset == gradient_count:
            inverse_hessian = np.eye(n, dtype=np.float64)
        old_parameters = parameters.copy()
        old_gradient = gradient.copy()
        direction = -(inverse_hessian @ gradient)
        gradient_projection = float(direction @ gradient)
        if gradient_projection < 0.0:
            step_length = 1.0
            accepted = False
            while True:
                parameters = old_parameters + step_length * direction
                count = int(np.count_nonzero(10.0 + old_parameters == 10.0 + parameters))
                if count < n:
                    value = float(function(parameters))
                    accepted = bool(
                        np.isfinite(value)
                        and value <= minimum + gradient_projection * step_length * 0.0001
                    )
                    if not accepted:
                        step_length *= 0.2
                if count == n or accepted:
                    break
            enough = abs(value - minimum) > relative_tolerance * (abs(minimum) + relative_tolerance)
            if not enough:
                count = n
                minimum = value
            if count < n:
                minimum = value
                gradient = _r_numeric_gradient(function, parameters)
                gradient_count += 1
                iteration += 1
                displacement = step_length * direction
                gradient_delta = gradient - old_gradient
                curvature = float(displacement @ gradient_delta)
                if curvature > 0.0:
                    transformed_delta = inverse_hessian @ gradient_delta
                    multiplier = 1.0 + float(transformed_delta @ gradient_delta) / curvature
                    inverse_hessian += (
                        multiplier * np.outer(displacement, displacement)
                        - np.outer(transformed_delta, displacement)
                        - np.outer(displacement, transformed_delta)
                    ) / curvature
                else:
                    last_reset = gradient_count
            elif last_reset < gradient_count:
                count = 0
                last_reset = gradient_count
        else:
            count = 0
            if last_reset == gradient_count:
                count = n
            else:
                last_reset = gradient_count
        if iteration >= maxit:
            break
        if gradient_count - last_reset > 2 * n:
            last_reset = gradient_count
        if count == n and last_reset == gradient_count:
            break
    return parameters

def ar1_ml_forecast(series: object) -> float:
    x = _as_finite_float64(series, ndim=1, name="AR(1) series")
    n = x.size
    if n < 2:
        raise ValueError("AR(1) series is too short")
    initial_mean = float(np.mean(x))
    residual = x - initial_mean
    mean_scale = 10.0 * np.sqrt(float(residual @ residual) / (n - 1.0)) / np.sqrt(float(n))
    if not np.isfinite(mean_scale) or mean_scale <= 0.0:
        raise ValueError("Degenerate AR(1) initialization scale")
    parameter_scale = np.array([1.0, mean_scale], dtype=np.float64)

    def objective(scaled_parameters: np.ndarray) -> float:
        raw_phi, mean = scaled_parameters * parameter_scale
        phi = float(np.tanh(raw_phi))
        if abs(phi) >= 1.0:
            return float("inf")
        centered = x - mean
        innovations = centered[1:] - phi * centered[:-1]
        sum_squares = (1.0 - phi * phi) * centered[0] ** 2 + float(innovations @ innovations)
        if sum_squares <= 0.0 or not np.isfinite(sum_squares):
            return float("inf")
        return 0.5 * (np.log(sum_squares / n) - np.log1p(-(phi * phi)) / n)

    # Use the R-compatible vmmin port so the AR(1) forecast matches the baseline.
    raw_phi, mean = (
        _r_vmmin(objective, np.array([0.0, initial_mean / mean_scale], dtype=np.float64))
        * parameter_scale
    )
    phi = float(np.tanh(raw_phi))
    return float(mean + phi * (x[-1] - mean))

def _limit_blas_threads() -> None:
    global _BLAS_LIMITER
    if _BLAS_LIMITER is not None:
        return
    try:
        from threadpoolctl import threadpool_limits

        _BLAS_LIMITER = threadpool_limits(limits=1, user_api="blas")
    except ImportError:
        pass

def _finite_matrix(value: ArrayLike, name: str, min_rows: int = 1) -> FloatArray:
    answer = np.asarray(value, dtype=np.float64)
    if answer.ndim != 2 or answer.shape[0] < min_rows or answer.shape[1] < 1:
        raise ValueError(f"{name} must be a two-dimensional numeric matrix")
    if not np.isfinite(answer).all():
        raise ValueError(f"{name} must be finite")
    return answer

def _finite_vector(value: ArrayLike, name: str) -> FloatArray:
    answer = np.asarray(value, dtype=np.float64)
    if answer.ndim == 2 and answer.shape[1] == 1:
        answer = answer[:, 0]
    if answer.ndim != 1 or answer.size < 1 or (not np.isfinite(answer).all()):
        raise ValueError(f"{name} must be a finite numeric vector")
    return answer

def transform_helmert_design(X: ArrayLike) -> FloatArray:
    matrix = _finite_matrix(X, "X")
    n_assets = matrix.shape[1]
    if n_assets < 2:
        raise ValueError("Helmert design needs at least two assets")
    # Helmert coordinates span the zero-sum asset subspace without adding a constraint.
    p = n_assets - 1
    denominator = np.sqrt(
        np.arange(1, p + 1, dtype=np.float64) * np.arange(2, p + 2, dtype=np.float64)
    )
    design = np.empty((matrix.shape[0], p), dtype=np.float64)
    prefix = matrix[:, 0].copy()
    for j in range(p):
        design[:, j] = ((j + 1) * matrix[:, j + 1] - prefix) / denominator[j]
        if j + 1 < p:
            prefix += matrix[:, j + 1]
    return design

def recover_helmert_weights(beta: ArrayLike) -> FloatArray:
    coefficients = np.asarray(beta, dtype=np.float64)
    vector_input = coefficients.ndim == 1
    if vector_input:
        coefficients = coefficients[:, None]
    if (
        coefficients.ndim != 2
        or coefficients.shape[0] < 1
        or coefficients.shape[1] < 1
        or (not np.isfinite(coefficients).all())
    ):
        raise ValueError("beta must be a finite vector or matrix with at least one row")
    # Map coordinate coefficients back to asset weights whose sum is zero.
    p = coefficients.shape[0]
    denominator = np.sqrt(
        np.arange(1, p + 1, dtype=np.float64) * np.arange(2, p + 2, dtype=np.float64)
    )
    scaled = coefficients / denominator[:, None]
    weights = np.zeros((p + 1, coefficients.shape[1]), dtype=np.float64)
    weights[p, :] = p * scaled[p - 1, :]
    suffix = np.zeros(coefficients.shape[1], dtype=np.float64)
    for j in range(p - 1, -1, -1):
        suffix += scaled[j, :]
        weights[j, :] = -suffix
        if j > 0:
            weights[j, :] += j * scaled[j - 1, :]
    return weights[:, 0] if vector_input else weights

def build_lambda_grid(
    X: ArrayLike,
    y: ArrayLike,
    alpha_grid: ArrayLike = ALPHA_GRID,
    n_lambda: int = N_LAMBDA,
    lambda_min_ratio: float = LAMBDA_MIN_RATIO,
    alpha_zero_reference: float = ALPHA_ZERO_REFERENCE,
) -> tuple[FloatArray, ...]:
    matrix = _finite_matrix(X, "X")
    response = _finite_vector(y, "y")
    alphas = _finite_vector(alpha_grid, "alpha_grid")
    if response.size != matrix.shape[0]:
        raise ValueError("X and y row counts differ")
    if not isinstance(n_lambda, (int, np.integer)) or n_lambda < 2:
        raise ValueError("n_lambda must be an integer >= 2")
    if not 0.0 < lambda_min_ratio < 1.0:
        raise ValueError("lambda_min_ratio must lie in (0, 1)")
    if not 0.0 < alpha_zero_reference <= 1.0:
        raise ValueError("alpha_zero_reference must lie in (0, 1]")
    if (alphas < 0).any() or (alphas > 1).any():
        raise ValueError("alpha_grid must lie in [0, 1]")
    gradient = matrix.T @ response / matrix.shape[0]
    spread = float(np.max(gradient) - np.min(gradient))
    centered_norm = float(np.linalg.norm(gradient - np.mean(gradient)))
    paths: list[FloatArray] = []
    for alpha in alphas:
        # alpha=0 reuses the alpha=.01 scale only to define a finite ridge lambda grid.
        path_alpha = alpha_zero_reference if alpha == 0.0 else float(alpha)
        if spread == 0.0 or centered_norm == 0.0:
            paths.append(np.zeros(n_lambda, dtype=np.float64))
            continue
        upper = centered_norm / path_alpha
        lower = spread / (2.0 * path_alpha) * lambda_min_ratio
        path = np.exp(np.linspace(math.log(upper), math.log(lower), n_lambda))
        path[0] = upper
        paths.append(path)
    return tuple(paths)

class _RMersenneTwister:
    _N = 624
    _M = 397
    _MASK32 = 4294967295
    _MATRIX_A = 2567483615
    _UPPER_MASK = 2147483648
    _LOWER_MASK = 2147483647

    def __init__(self, seed: int):
        if not isinstance(seed, (int, np.integer)) or not 0 <= int(seed) <= 2147483647:
            raise ValueError("seed must be an R-compatible nonnegative integer")
        state_seed = int(seed) & self._MASK32
        for _ in range(50):
            state_seed = 69069 * state_seed + 1 & self._MASK32
        raw = np.empty(self._N + 1, dtype=np.uint32)
        for j in range(self._N + 1):
            state_seed = 69069 * state_seed + 1 & self._MASK32
            raw[j] = state_seed
        self._state = raw[1:].copy()
        self._index = self._N

    def _uint32(self) -> int:
        if self._index >= self._N:
            state = self._state
            for kk in range(self._N - self._M):
                y = int(state[kk]) & self._UPPER_MASK | int(state[kk + 1]) & self._LOWER_MASK
                state[kk] = int(state[kk + self._M]) ^ y >> 1 ^ (self._MATRIX_A if y & 1 else 0)
            for kk in range(self._N - self._M, self._N - 1):
                y = int(state[kk]) & self._UPPER_MASK | int(state[kk + 1]) & self._LOWER_MASK
                state[kk] = (
                    int(state[kk + self._M - self._N]) ^ y >> 1 ^ (self._MATRIX_A if y & 1 else 0)
                )
            y = int(state[self._N - 1]) & self._UPPER_MASK | int(state[0]) & self._LOWER_MASK
            state[self._N - 1] = int(state[self._M - 1]) ^ y >> 1 ^ (self._MATRIX_A if y & 1 else 0)
            self._index = 0
        y = int(self._state[self._index])
        self._index += 1
        y ^= y >> 11
        y ^= y << 7 & 2636928640
        y ^= y << 15 & 4022730752
        y ^= y >> 18
        return y & self._MASK32

    def uniform(self) -> float:
        value = self._uint32() * (1.0 / 4294967296.0)
        if value <= 0.0:
            return 0.5 / 4294967295.0
        if 1.0 - value <= 0.0:
            return 1.0 - 0.5 / 4294967295.0
        return value

    def uniform_index(self, n: int) -> int:
        if not isinstance(n, (int, np.integer)) or n < 1:
            raise ValueError("n must be a positive integer")
        bits = int(math.ceil(math.log2(n)))
        mask = (1 << bits) - 1
        while True:
            value = 0
            for _ in range(0, bits + 1, 16):
                value = 65536 * value + int(math.floor(self.uniform() * 65536.0))
            value &= mask
            if value < n:
                return value

def r_sample(n: int, seed: int = CV_SEED) -> IntArray:
    if not isinstance(n, (int, np.integer)) or n < 1:
        raise ValueError("n must be a positive integer")
    rng = _RMersenneTwister(seed)
    pool = np.arange(n, dtype=np.int64)
    result = np.empty(n, dtype=np.int64)
    remaining = n
    for i in range(n):
        j = rng.uniform_index(remaining)
        result[i] = pool[j]
        remaining -= 1
        pool[j] = pool[remaining]
    return result

def build_cv_folds(n: int, n_folds: int = CV_FOLDS, seed: int = CV_SEED) -> tuple[IntArray, ...]:
    if not isinstance(n_folds, (int, np.integer)) or n_folds < 2 or n_folds > n // 2:
        raise ValueError("n_folds must lie between 2 and floor(n / 2)")
    permutation = r_sample(n, seed)
    return tuple((permutation[k::n_folds].copy() for k in range(n_folds)))

def solve_ridge_path(Z: ArrayLike, y: ArrayLike, lambdas: ArrayLike) -> FloatArray:
    design = _finite_matrix(Z, "Z", min_rows=2)
    response = _finite_vector(y, "y")
    penalty = _finite_vector(lambdas, "lambdas")
    if response.size != design.shape[0] or (penalty <= 0).any():
        raise ValueError("Invalid ridge path inputs")
    _limit_blas_threads()
    # Solve the no-intercept ridge path
    kernel = design @ design.T / design.shape[0]
    values, vectors = np.linalg.eigh((kernel + kernel.T) / 2.0)
    values = np.maximum(values, 0.0)
    dual = vectors @ ((vectors.T @ response)[:, None] / (values[:, None] + penalty[None, :]))
    beta = design.T @ dual / design.shape[0]
    if not np.isfinite(beta).all():
        raise RuntimeError("nonfinite ridge path")
    return beta

def solve_positive_enet_path(
    Z: ArrayLike, y: ArrayLike, lambdas: ArrayLike, alpha: float, *, max_iter: int = 1000000
) -> FloatArray:
    design = _finite_matrix(Z, "Z", min_rows=2)
    response = _finite_vector(y, "y")
    penalty = _finite_vector(lambdas, "lambdas")
    if response.size != design.shape[0] or not 0.0 < alpha < 1.0:
        raise ValueError("Invalid positive elastic-net inputs")
    if (penalty <= 0).any() or np.any(np.diff(penalty) > 0.0):
        raise ValueError("elastic-net lambdas must be positive and descending")
    from scipy.linalg import cho_factor, cho_solve

    _limit_blas_threads()
    design = np.asfortranarray(design)
    t, p = design.shape
    q = 1.0 - float(alpha)
    active = np.zeros(p, dtype=bool)
    signs = np.zeros(p, dtype=np.float64)
    kernel = np.zeros((t, t), dtype=np.float64)
    signed_sum = np.zeros(t, dtype=np.float64)
    beta = np.empty((p, penalty.size), dtype=np.float64)
    gradient_zero = design.T @ response / t
    eps = np.finfo(np.float64).eps
    gamma_t = t * eps / (1.0 - t * eps)
    gradient_zero_error = gamma_t * (np.abs(design).T @ np.abs(response)) / t + 4.0 * eps * np.abs(
        gradient_zero
    )
    # Active-set updates solve the same no-intercept elastic-net objective as glmnet.
    for path_index, lam in enumerate(penalty):
        if path_index and path_index % 25 == 0:
            active_design = design[:, active]
            kernel = active_design @ active_design.T
        threshold = lam * alpha
        ridge = lam * q
        zero_limit = threshold + gradient_zero_error + 4.0 * abs(np.spacing(threshold))
        if not np.any(active) and np.all(np.abs(gradient_zero) <= zero_limit):
            beta[:, path_index] = 0.0
            continue
        seen: set[tuple[bytes, bytes]] = set()
        for _ in range(max_iter):
            key = (active.tobytes(), signs.tobytes())
            if key in seen:
                raise RuntimeError(f"elastic-net active set cycled at lambda {path_index}")
            seen.add(key)
            if not np.any(active):
                dual_residual = response
            else:
                system = kernel / (t * ridge)
                system.flat[:: t + 1] += 1.0
                right = response + alpha / q * signed_sum
                factor = cho_factor(system, lower=True, overwrite_a=True, check_finite=False)
                dual_residual = cho_solve(factor, right, check_finite=False)
            v = design.T @ dual_residual / t
            new_active = np.abs(v) > threshold
            new_signs = np.where(new_active, np.sign(v), 0.0)
            if np.array_equal(new_active, active) and np.array_equal(new_signs, signs):
                break
            added, dropped = (new_active & ~active, active & ~new_active)
            if np.any(added):
                changed = design[:, added]
                kernel += changed @ changed.T
            if np.any(dropped):
                changed = design[:, dropped]
                kernel -= changed @ changed.T
            active, signs = (new_active, new_signs)
            signed_sum = design[:, active] @ signs[active]
        else:
            raise RuntimeError(f"elastic-net active set exceeded max_iter at lambda {path_index}")
        soft = np.sign(v) * np.maximum(np.abs(v) - threshold, 0.0)
        beta[:, path_index] = soft / ridge
    if not np.isfinite(beta).all():
        raise RuntimeError("nonfinite elastic-net path")
    return beta

def solve_elastic_net_path(
    Z: ArrayLike, y: ArrayLike, lambdas: ArrayLike, alpha: float, *, max_iter: int = 1000000
) -> FloatArray:
    if not np.isfinite(alpha) or not 0.0 <= alpha < 1.0:
        raise ValueError("alpha must lie in [0, 1)")
    if alpha == 0.0:
        return solve_ridge_path(Z, y, lambdas)
    return solve_positive_enet_path(Z, y, lambdas, float(alpha), max_iter=max_iter)

def prepare_poet_problem(R, S, mu, K):
    """Construct the synthetic regression, or return None for a zero signal."""
    returns = _as_finite_float64(R, ndim=2, name="MAXSER-PR inputs")
    signals = _as_finite_float64(S, ndim=2, name="MAXSER-PR inputs")
    mean = _as_finite_float64(mu, ndim=1, name="MAXSER-PR inputs")
    t, n = returns.shape
    if signals.shape != returns.shape or mean.size != n or t < 6 or n < 2:
        raise ValueError("Invalid MAXSER-PR inputs")
    if not isinstance(K, int) or K < 0 or K >= min(t, n):
        raise ValueError("POET factor number is inadmissible; no silent clipping")
    if np.all(signals - signals.mean(axis=1, keepdims=True) == 0.0):
        return None
    raw_e = returns - signals
    gamma_e = raw_e - raw_e.mean(axis=0, keepdims=True)
    gamma = np.maximum(gamma_path(gamma_e + mean, mean), _SQRT_DOUBLE_EPS)
    gamma_next = np.sqrt(max(np.exp(ar1_ml_forecast(np.log(gamma**2))), _SQRT_DOUBLE_EPS))
    standardized = raw_e / gamma[:, None]
    x = (standardized - standardized.mean(axis=0, keepdims=True)) * gamma_next + mean
    theta = zero_invest_theta(np.vstack((signals, mean)), poet_covariance(standardized, K=K),
                              np.concatenate((gamma, np.array([gamma_next]))))
    xi = float(np.mean(theta[:t] / (1 + theta[:t])))
    if not np.isfinite(xi) or xi >= 1:
        raise ValueError(f"Invalid MAXSER-PR xi at K={K}")
    if xi <= 1e-14:
        return None
    response = TARGET_DAILY_SIGMA / np.sqrt(xi * (1 - xi))
    return x, np.full(t, response, dtype=float)

class ProtectedBasis:
    def __init__(self, dictionary, prediction, max_rank=64, strength=1.0):
        self.n_assets = len(prediction)
        self.strength = strength
        centered = prediction - np.mean(prediction)
        norm = np.linalg.norm(centered)
        has_signal = norm > 32 * np.finfo(float).eps * max(np.linalg.norm(prediction), np.finfo(float).tiny)
        q = centered / norm if has_signal else None
        residual = np.asarray(dictionary, dtype=float).copy()
        if not np.isfinite(residual).all():
            raise ValueError("nonfinite basis dictionary")
        residual -= residual.mean(axis=0)
        if q is not None:
            residual -= q[:, None] * (q @ residual)[None, :]
        left, values, _ = np.linalg.svd(residual, full_matrices=False)
        tolerance = max(residual.shape) * np.finfo(float).eps * (values[0] if len(values) else 0)
        take = min(max_rank - int(has_signal), self.n_assets - 1 - int(has_signal), int((values > tolerance).sum()))
        # Reorthogonalize to protect the exact signal even with nearly singular dictionaries.
        columns = [] if q is None else [q]
        for v in left[:, :take].T:
            v = v - v.mean()
            for u in columns:
                v -= u * (u @ v)
            vnorm = np.linalg.norm(v)
            if vnorm > 1e-12:
                columns.append(v / vnorm)
        self.directions = np.column_stack(columns) if columns else np.empty((self.n_assets, 0))
        self.rank = len(columns)
        self.scale = 1 / np.sqrt(1 + self.rank * ((1 + strength)**2 - 1) / (self.n_assets - 1))
        self.coordinates = transform_helmert_design(self.directions.T)

    def transform(self, x):
        return self.scale * (transform_helmert_design(x) +
                             self.strength * (x @ self.directions) @ self.coordinates)

    def recover(self, beta):
        w = recover_helmert_weights(beta)
        return self.scale * (w + self.strength * self.directions @ (self.directions.T @ w))

class MomentBasis:
    """B=s B_SHAP (I+3 B_SHAP' F F' B_SHAP/lambda_max)^(-1/2)."""
    def __init__(self, original, factor, maximum_rank=5):
        self.original, self.n_assets = original, original.n_assets
        factor = np.asarray(factor, dtype=float)
        if factor.ndim != 2 or factor.shape[0] != self.n_assets or not np.isfinite(factor).all():
            raise ValueError("Invalid moment factor")
        coordinates = original.transform(factor.T).T
        left, singular, _ = np.linalg.svd(coordinates, full_matrices=False)
        top = float(singular[0]) if singular.size else 0.
        keep = singular > max(coordinates.shape)*np.finfo(float).eps*top
        self.v = left[:, keep]
        if self.v.shape[1] > maximum_rank:
            raise RuntimeError("Uncertainty rank exceeds estimator contract")
        self.d = 1/np.sqrt(1 + 3*(singular[keep]/top)**2) if top else np.empty(0)
        recovered = original.recover(self.v) if len(self.d) else np.empty((self.n_assets, 0))
        p = self.n_assets-1
        trace = p + float(np.sum((self.d**2-1)*np.sum(recovered**2, axis=0)))
        if not np.isfinite(trace) or trace <= 0:
            raise RuntimeError("Nonpositive moment-basis trace")
        self.scale = float(np.sqrt(p/trace))

    def transform(self, x):
        z = self.original.transform(x)
        return self.scale*(z + ((z @ self.v)*(self.d-1)) @ self.v.T)

    def recover(self, beta):
        beta = np.asarray(beta)
        product = self.v.T @ beta
        correction = (self.d-1)*product if beta.ndim == 1 else (self.d-1)[:, None]*product
        return self.scale*self.original.recover(beta + self.v @ correction)

def kkt_error(design, y, beta, lambdas, alpha):
    gradient = design.T @ (design @ beta - y[:, None]) / len(y)
    gradient += beta * (lambdas * (1 - alpha))[None, :]
    threshold = (lambdas * alpha)[None, :]
    error = np.where(beta != 0, np.abs(gradient + threshold * np.sign(beta)),
                     np.maximum(np.abs(gradient) - threshold, 0))
    scale = max(float(np.max(np.abs(design.T @ y / len(y)))), np.finfo(float).tiny)
    return float(error.max() / scale)

def solve_path(design, y, lambdas, alpha):
    if alpha == 0:
        beta = solve_ridge_path(design, y, lambdas)
        error = kkt_error(design, y, beta, lambdas, alpha)
        if not np.isfinite(error) or error > 1e-7:
            raise RuntimeError(f"Ridge KKT failed: {error}")
        return beta
    if alpha < 1:
        try:
            beta = solve_elastic_net_path(design, y, lambdas, alpha)
            error = kkt_error(design, y, beta, lambdas, alpha)
            if error <= 1e-7:
                return beta
        except RuntimeError:
            pass
    from sklearn.linear_model import enet_path
    with warnings.catch_warnings(record=True):
        _, beta, _ = enet_path(
            np.asfortranarray(design), np.asarray(y), l1_ratio=alpha,
            alphas=np.asarray(lambdas), tol=1e-12, max_iter=1000000,
            precompute=False, selection="cyclic")
    error = kkt_error(design, y, beta, lambdas, alpha)
    if not np.isfinite(beta).all() or not np.isfinite(error) or error > 1e-7:
        raise RuntimeError(f"coordinate-descent KKT failed: alpha={alpha}, error={error}")
    return beta

def fit_cv(X, y, basis, alpha_grid=ALPHA_GRID, *, n_lambda=N_LAMBDA,
           n_folds=CV_FOLDS, seed=CV_SEED):
    alphas = np.asarray(alpha_grid, dtype=float)
    if (alphas < 0).any() or (alphas > 1).any() or (np.diff(alphas) <= 0).any():
        raise ValueError("alphas must increase within [0,1]")
    count = n_lambda
    folds = build_cv_folds(len(X), n_folds, seed)
    grids = build_lambda_grid(X, y, alphas, count)
    scores = np.full((len(alphas), count, len(folds)), -np.inf)
    coverage = np.zeros(len(X), dtype=int)
    for k, omit in enumerate(folds):
        coverage[omit] += 1
        keep = np.ones(len(X), dtype=bool); keep[omit] = False
        design = basis.transform(X[keep])
        for a, alpha in enumerate(alphas):
            beta = solve_path(design, y[keep], grids[a], float(alpha))
            held = X[omit] @ basis.recover(beta)
            with np.errstate(divide="ignore", invalid="ignore"):
                score = (held.mean(axis=0) / held.std(axis=0, ddof=1)) ** 2
            score[~np.isfinite(score)] = -np.inf
            scores[a, :, k] = score
    if not np.all(coverage == 1):
        raise RuntimeError("CV rows must occur exactly once")
    mean = np.full((len(alphas), count), -np.inf)
    valid = np.all(np.isfinite(scores), axis=2)
    for a in range(len(alphas)):
        mean[a, valid[a]] = scores[a, valid[a]].mean(axis=1)
        if not np.isfinite(mean[a]).any():
            raise RuntimeError(f"no finite CV candidate at alpha={alphas[a]}")
    indices = mean.argmax(axis=1)
    a = int(mean[np.arange(len(alphas)), indices].argmax()); j = int(indices[a])
    beta = solve_path(basis.transform(X), y, grids[a][j:j+1], float(alphas[a]))[:, 0]
    return basis.recover(beta)

@lru_cache(maxsize=8)
def quadrature(depth, extra_nodes=0):
    offsets = np.zeros((depth + 1, depth + 1, 2), dtype=np.int32)
    lengths = np.zeros_like(offsets)
    nodes, weights = [], []
    for c in range(depth + 1):
        for h in range(depth + 1 - c):
            for kind in range(2):
                if (kind == 0 and h == 0) or (kind == 1 and c == 0):
                    continue
                exponent = c if kind == 0 else c - 1
                degree = h - 1 if kind == 0 else h
                q = (degree + 2) // 2 + extra_nodes
                roots, raw_weights = roots_jacobi(q, exponent, 0)
                offsets[c, h, kind], lengths[c, h, kind] = len(nodes), q
                nodes.extend((roots + 1) / 2)
                weights.extend(raw_weights / (2.0 ** (exponent + 1)))
    return (np.asarray(nodes, dtype=float), np.asarray(weights, dtype=float),
            offsets.ravel(), lengths.ravel())

class StableTreeSHAP:
    def __init__(self, dump):
        if dump.get("objective", "regression").split()[0] != "regression":
            raise ValueError("only frozen regression boosters are supported")
        if dump.get("average_output", False) or dump.get("num_tree_per_iteration", 1) != 1:
            raise ValueError("unsupported booster aggregation")
        self.nfeature = len(dump["feature_names"])
        columns = {key: [] for key in ["feature", "left", "right", "slot", "missing",
                                       "default_left", "threshold", "count", "value"]}
        roots, max_depth, max_tree_depth = [], 0, 0

        def flatten(node, path, root_count, tree_depth=0):
            nonlocal max_depth, max_tree_depth
            index = len(columns["feature"])
            for column in columns.values():
                column.append(0)
            max_depth = max(max_depth, len(path))
            max_tree_depth = max(max_tree_depth, tree_depth)
            if "leaf_index" in node or "split_feature" not in node:
                if node.get("leaf_coeff"):
                    raise ValueError("linear leaves are unsupported")
                columns["feature"][index] = -1
                columns["count"][index] = node.get("leaf_count", root_count)
                columns["value"][index] = float(node["leaf_value"])
                return index
            if node["decision_type"] != "<=":
                raise ValueError("categorical split is unsupported")
            feature = int(node["split_feature"])
            if not 0 <= feature < self.nfeature:
                raise ValueError("feature index out of range")
            columns["feature"][index] = feature
            columns["slot"][index] = list(path).index(feature) if feature in path else len(path)
            columns["missing"][index] = {"None": 0, "Zero": 1, "NaN": 2}[node["missing_type"]]
            columns["default_left"][index] = int(node["default_left"])
            columns["threshold"][index] = float(node["threshold"])
            columns["count"][index] = node["internal_count"]
            branch = dict(path)
            branch[feature] = True
            left = flatten(node["left_child"], branch, root_count, tree_depth + 1)
            right = flatten(node["right_child"], branch, root_count, tree_depth + 1)
            columns["left"][index], columns["right"][index] = left, right
            if columns["count"][index] <= 0 or columns["count"][left] <= 0 or columns["count"][right] <= 0:
                raise ValueError("nonpositive training count")
            if columns["count"][index] != columns["count"][left] + columns["count"][right]:
                raise ValueError("inconsistent node training counts")
            return index

        for tree in dump["tree_info"]:
            root = tree["tree_structure"]
            root_count = root.get("internal_count", root.get("leaf_count", 1)) or 1
            roots.append(flatten(root, {}, root_count))
        self.roots = np.asarray(roots, dtype=np.int32)
        self.depth = max_depth
        self.tree_depth = max_tree_depth
        self.columns = {key: np.asarray(value, dtype=np.float64 if key in ("threshold", "count", "value") else np.int32)
                        for key, value in columns.items()}

    def predict_contrib(self, features, extra_nodes=0):
        features = np.ascontiguousarray(features, dtype=np.float64)
        if features.ndim != 2 or features.shape[1] != self.nfeature or np.isinf(features).any():
            raise ValueError("invalid feature matrix")
        nodes, weights, offsets, lengths = quadrature(self.depth, extra_nodes)
        c = self.columns
        return _shap_kernel(features, self.roots, c["feature"], c["left"], c["right"],
            c["slot"], c["missing"], c["default_left"], c["threshold"], c["count"], c["value"],
            nodes, weights, offsets, lengths, self.depth, self.tree_depth)


@njit(fastmath=False, cache=False, boundscheck=True)
def _shap_kernel(x, roots, feature, left, right, slot, missing, default_left,
                 threshold, count, value, nodes, weights, offsets, lengths, depth, tree_depth):
    nrow, nfeature = x.shape
    phi = np.zeros((nrow, nfeature), dtype=np.float64)
    stride = depth + 1
    path_feature = np.zeros(stride, dtype=np.int32)
    one = np.zeros(stride, dtype=np.int32)
    zero = np.zeros(stride, dtype=np.float64)
    hot = np.zeros(stride, dtype=np.int32)
    cold = np.zeros(stride, dtype=np.int32)
    prefix = np.zeros(stride + 1, dtype=np.float64)
    factor = np.zeros(stride, dtype=np.float64)
    # Repeated feature splits increase traversal depth without increasing path rank.
    stack_size = tree_depth + 1
    stack_node = np.zeros(stack_size, dtype=np.int32)
    stack_depth = np.zeros(stack_size, dtype=np.int32)
    stage = np.zeros(stack_size, dtype=np.int32)
    prior_zero = np.zeros(stack_size, dtype=np.float64)
    prior_one = np.zeros(stack_size, dtype=np.int32)
    go_left = np.zeros(stack_size, dtype=np.int32)
    zero_missing_threshold = float(np.float32(1e-35))
    for row in range(nrow):
        for root in roots:
            level = 0
            stack_node[0], stack_depth[0], stage[0] = root, 0, 0
            while level >= 0:
                node, d = stack_node[level], stack_depth[level]
                if feature[node] < 0:
                    h, c, cold_product = 0, 0, 1.0
                    for j in range(d):
                        if one[j]:
                            hot[h] = j
                            h += 1
                        else:
                            cold[c] = j
                            c += 1
                            cold_product *= zero[j]
                    v = value[node] * cold_product
                    if h > 0:
                        key = 2 * (c * stride + h)
                        for q in range(offsets[key], offsets[key] + lengths[key]):
                            t = nodes[q]
                            prefix[0] = 1.0
                            for j in range(h):
                                z = zero[hot[j]]
                                factor[j] = z + (1.0 - z) * t
                                prefix[j + 1] = prefix[j] * factor[j]
                            suffix = 1.0
                            for j in range(h - 1, -1, -1):
                                p = hot[j]
                                phi[row, path_feature[p]] += v * weights[q] * prefix[j] * suffix * (1.0 - zero[p])
                                suffix *= factor[j]
                    if c > 0:
                        integral = 0.0
                        key = 2 * (c * stride + h) + 1
                        for q in range(offsets[key], offsets[key] + lengths[key]):
                            product, t = 1.0, nodes[q]
                            for j in range(h):
                                z = zero[hot[j]]
                                product *= z + (1.0 - z) * t
                            integral += weights[q] * product
                        contribution = -v * integral
                        for j in range(c):
                            phi[row, path_feature[cold[j]]] += contribution
                    level -= 1
                    continue
                p = slot[node]
                if stage[level] == 0:
                    fresh = p == d
                    if fresh:
                        path_feature[p], zero[p], one[p] = feature[node], 1.0, 1
                    prior_zero[level], prior_one[level] = zero[p], one[p]
                    fval = x[row, feature[node]]
                    if np.isnan(fval) and missing[node] != 2:
                        fval = 0.0
                    is_missing = ((missing[node] == 1 and abs(fval) <= zero_missing_threshold)
                                  or (missing[node] == 2 and np.isnan(fval)))
                    goes_left = bool(default_left[node]) if is_missing else fval <= threshold[node]
                    go_left[level] = int(goes_left)
                    zero[p] = prior_zero[level] * count[left[node]] / count[node]
                    one[p] = int(bool(prior_one[level]) and goes_left)
                    stage[level] = 1
                    stack_node[level + 1] = left[node]
                    stack_depth[level + 1] = d + int(fresh)
                    stage[level + 1] = 0
                    level += 1
                elif stage[level] == 1:
                    zero[p] = prior_zero[level] * count[right[node]] / count[node]
                    one[p] = int(bool(prior_one[level]) and not bool(go_left[level]))
                    stage[level] = 2
                    stack_node[level + 1] = right[node]
                    stack_depth[level + 1] = d + int(p == d)
                    stage[level + 1] = 0
                    level += 1
                else:
                    zero[p], one[p] = prior_zero[level], prior_one[level]
                    level -= 1
    return phi


# Independent annual fits, followed by independent monthly portfolios.
TRAINING_JOBS = 4
TRAINING_THREADS = 8
PORTFOLIO_JOBS = 32
SEED = 42
REFIT_EVERY = 12
MIN_TRAIN_MONTHS = 36
SAMPLE_WINDOW_DAYS = 756
MIN_DAILY_OBS = 250
TARGET_COL = "ret_exc_lead1m"
LGBM_PARAMS = dict(objective="regression", n_estimators=1000, learning_rate=0.01,
                   num_leaves=255, min_child_samples=2000, random_state=SEED,
                   n_jobs=TRAINING_THREADS, deterministic=True, force_row_wise=True,
                   device_type="cpu", verbosity=-1)


def is_test_row(flag):
    if pd.api.types.is_bool_dtype(flag):
        return flag.fillna(False)
    if pd.api.types.is_numeric_dtype(flag):
        return flag.fillna(0).ne(0)
    return flag.astype(str).str.strip().str.lower().isin(["1", "1.0", "true", "t"])


def prepare_panel(chars, features):
    names = [str(name) for name in features["features"]]
    if not names or len(names) != len(set(names)):
        raise ValueError("features must contain distinct column names")
    columns = ["id", "eom", TARGET_COL, *names]
    columns += [name for name in ["eom_ret", "ctff_test"] if name in chars.columns]
    panel = chars.loc[:, columns].copy()
    panel["eom"] = pd.to_datetime(panel.eom)
    expected = panel.eom + pd.offsets.MonthEnd(1)
    if "eom_ret" in panel:
        supplied = pd.to_datetime(panel.eom_ret)
        if (supplied.notna() & supplied.ne(expected)).any():
            raise ValueError("Forward-label dates must end in the following month")
        panel["eom_ret"] = supplied.fillna(expected)
    else:
        panel["eom_ret"] = expected
    panel = panel.sort_values(["eom", "id"], kind="stable").reset_index(drop=True)
    if panel.duplicated(["eom", "id"]).any() or panel[["id", "eom"]].isna().any().any():
        raise ValueError("Missing or duplicate stock-month keys")
    if not np.equal(panel.id.to_numpy(), np.floor(panel.id.to_numpy())).all():
        raise ValueError("CTF identifiers must be integers")
    panel["id"] = panel.id.astype(np.int64)
    return panel, names


def prediction_schedule(panel):
    month, month_ix = np.unique(panel.eom.to_numpy(dtype="datetime64[D]"), return_inverse=True)
    starts = np.searchsorted(month_ix, np.arange(len(month)), side="left")
    ends = np.searchsorted(month_ix, np.arange(len(month)), side="right")
    if len(month) <= MIN_TRAIN_MONTHS:
        raise ValueError("At least 37 monthly observations are required")
    flag = is_test_row(panel.ctff_test).to_numpy() if "ctff_test" in panel else np.zeros(len(panel), dtype=bool)
    target = np.flatnonzero(np.bincount(month_ix, weights=flag, minlength=len(month)) > 0)
    if not len(target):
        target = np.arange(MIN_TRAIN_MONTHS, len(month))
    if (target < MIN_TRAIN_MONTHS).any():
        raise ValueError("Test months precede the minimum training history")
    tasks = []
    for t in range(MIN_TRAIN_MONTHS, len(month), REFIT_EVERY):
        stop = min(t + REFIT_EVERY, len(month))
        tasks.append(dict(index=len(tasks), t=t, train_end=int(starts[t]),
                          lo=int(starts[t]), hi=int(ends[stop - 1]),
                          auxiliary=bool(np.isin(target, np.arange(t, stop)).any())))
    return month, month_ix, starts, ends, target, tasks


def train_refit(task, X, y, month_ix, names):
    """A fit's training rows and seeds are independent of process scheduling."""
    _limit_blas_threads()
    end, lo, hi = task["train_end"], task["lo"], task["hi"]
    if not np.isfinite(y[:end]).all():
        raise ValueError("Nonfinite training labels")
    model = lgb.LGBMRegressor(**LGBM_PARAMS)
    model.fit(X[:end], y[:end], feature_name=names)
    prediction = model.booster_.predict(X[lo:hi], num_threads=TRAINING_THREADS).astype(np.float32)
    result = dict(lo=lo, hi=hi, prediction=prediction,
                  model=model.booster_.model_to_string() if task["auxiliary"] else None,
                  auxiliary=None)
    if task["auxiliary"]:
        blocks = np.array_split(np.arange(task["t"]), 5)
        if all(task["t"] - len(block) >= 24 for block in blocks):
            auxiliary = np.empty((hi - lo, 5), dtype=float)
            for k, block in enumerate(blocks):
                keep = ~np.isin(month_ix[:end], block)
                deleted = lgb.LGBMRegressor(**LGBM_PARAMS)
                deleted.fit(X[:end][keep], y[:end][keep], feature_name=names)
                auxiliary[:, k] = (deleted.booster_.predict(X[lo:hi], num_threads=TRAINING_THREADS)
                                   .astype(np.float32).astype(float) * (12 / 252))
            result["auxiliary"] = auxiliary
    print(f"Annual fit {task['index'] + 1} completed", flush=True)
    return result


def train_models(panel, names):
    schedule = prediction_schedule(panel)
    _, month_ix, _, _, _, tasks = schedule
    X = panel[names].to_numpy(dtype=np.float32)
    group = panel.groupby("eom", sort=False)[TARGET_COL]
    y = np.arctan(((panel[TARGET_COL] - group.transform("mean")) / group.transform("std"))
                  .to_numpy(dtype=np.float32))
    # Four 8-thread LightGBM fits use the platform's 32 CPUs without nesting.
    with parallel_config(backend="loky", n_jobs=TRAINING_JOBS, inner_max_num_threads=TRAINING_THREADS):
        fits = Parallel(max_nbytes="1M", mmap_mode="r")(
            delayed(train_refit)(task, X, y, month_ix, names) for task in tasks)
    predictions = np.full(len(panel), np.nan, dtype=np.float32)
    for fit in fits:
        predictions[fit["lo"]:fit["hi"]] = fit["prediction"]
    return X, predictions, fits, schedule


def assemble_daily_returns(daily_ret, panel_ids):
    keep = np.isin(daily_ret.id.to_numpy(), panel_ids)
    ids = daily_ret.id.to_numpy()[keep]
    days = daily_ret.date.to_numpy()[keep].astype("datetime64[D]")
    values = daily_ret.ret_exc.to_numpy()[keep].astype(np.float64)
    dates, di = np.unique(days, return_inverse=True)
    assets, ai = np.unique(ids, return_inverse=True)
    if not len(dates) or not len(assets):
        raise ValueError("No daily observations for the supplied universe")
    if len(np.unique(di.astype(np.int64) * len(assets) + ai)) != len(values):
        raise ValueError("Duplicate stock-day observations")
    matrix = np.full((len(dates), len(assets)), np.nan, dtype=np.float64)
    matrix[di, ai] = values
    return dates, assets, matrix


def signal_history(dates, ids, historical):
    result = np.zeros((len(dates), len(ids)), dtype=np.float64)
    day_month = dates.astype("datetime64[M]")
    for month in np.unique(day_month):
        values = np.fromiter((historical.get((month, asset), 0.0) for asset in ids),
                             dtype=np.float64, count=len(ids))
        values[~np.isfinite(values)] = 0.0
        result[day_month == month, :] = values
    return result


def fit_portfolio(returns, history, mu, contribution, auxiliary):
    """The same K=10 problem, full-rank B, lambda paths and SR² CV as the frozen run."""
    _limit_blas_threads()
    K = min(10, min(returns.shape) - 1)
    problem = prepare_poet_problem(returns, history, mu, K)
    if problem is None:
        return np.zeros(len(mu))
    protected = ProtectedBasis(contribution, mu, max_rank=64, strength=1.0)
    basis = protected
    if auxiliary is not None:
        direct = auxiliary - mu[:, None]
        direct -= direct.mean(axis=0, keepdims=True)
        basis = MomentBasis(protected, direct / np.sqrt(5))
    weights = fit_cv(*problem, basis)
    if not np.isfinite(weights).all() or abs(weights.sum()) > 1e-10 + 1e-10 * np.abs(weights).sum():
        raise RuntimeError("Invalid zero-investment weights")
    return weights


def prepare_month_inputs(month, rows, prediction, daily_context, historical_context):
    dates, daily_ids, matrix = daily_context
    history_month, history_ids, history_prediction = historical_context
    ids = rows.id.to_numpy()
    end = int(np.searchsorted(dates, month, side="right"))
    start = max(0, end - SAMPLE_WINDOW_DAYS)
    columns = np.clip(np.searchsorted(daily_ids, ids), 0, len(daily_ids) - 1)
    present = daily_ids[columns] == ids
    eligible = np.zeros(len(ids), dtype=bool)
    eligible[np.flatnonzero(present)] = np.isfinite(matrix[start:end, columns[present]]).sum(axis=0) >= MIN_DAILY_OBS
    if end - start < 6 or eligible.sum() < 2:
        return eligible, None, None, None
    selected_ids = ids[eligible].astype(str)
    returns = matrix[start:end, columns[eligible]]
    returns = np.where(np.isfinite(returns), returns, 0.0).astype(np.float64, copy=False)
    needed = np.isin(history_month, np.unique(dates[start:end].astype("datetime64[M]")))
    historical = {(m, str(i)): float(p) * (12 / 252)
                  for m, i, p in zip(history_month[needed], history_ids[needed], history_prediction[needed])}
    history = signal_history(dates[start:end], selected_ids, historical)
    mu = prediction[eligible].astype(float) * (12 / 252)
    return eligible, returns, history, mu


def optimize_month(month, rows, features, prediction, auxiliary, model_text,
                   daily_context, historical_context):
    _limit_blas_threads()
    eligible, returns, history, mu = prepare_month_inputs(month, rows, prediction, daily_context, historical_context)
    weights = np.zeros(len(rows), dtype=np.float64)
    if returns is not None:
        booster = lgb.Booster(model_str=model_text)
        contributor = StableTreeSHAP(booster.dump_model())
        contribution = contributor.predict_contrib(features[eligible])
        auxiliary_eligible = None if auxiliary is None else auxiliary[eligible]
        weights[eligible] = fit_portfolio(returns, history, mu, contribution, auxiliary_eligible)
    output = rows[["id", "eom"]].copy()
    output["w"] = weights
    print(f"Portfolio {str(month)} completed", flush=True)
    return output


def build_portfolios(panel, daily_ret, trained):
    X, prediction, fits, schedule = trained
    months, _, starts, ends, target, _ = schedule
    daily_context = assemble_daily_returns(daily_ret, panel.id.unique())
    valid = np.isfinite(prediction)
    historical_context = (panel.eom_ret.to_numpy(dtype="datetime64[M]")[valid],
                          panel.id.to_numpy()[valid], prediction[valid])
    args = []
    for t in target:
        annual = (int(t) - MIN_TRAIN_MONTHS) // REFIT_EVERY
        fit = fits[annual]
        lo, hi = int(starts[t]), int(ends[t])
        aux = None if fit["auxiliary"] is None else fit["auxiliary"][lo - fit["lo"]:hi - fit["lo"]]
        args.append((months[t], panel.iloc[lo:hi][["id", "eom"]], X[lo:hi],
                     prediction[lo:hi], aux, fit["model"], daily_context, historical_context))
    # Complete months are independent; 32 processes each use one BLAS thread.
    with parallel_config(backend="loky", n_jobs=PORTFOLIO_JOBS, inner_max_num_threads=1):
        pieces = Parallel(max_nbytes="1M", mmap_mode="r")(delayed(optimize_month)(*a) for a in args)
    output = pd.concat(pieces, ignore_index=True)
    if "ctff_test" in panel and is_test_row(panel.ctff_test).any():
        keys = panel.loc[is_test_row(panel.ctff_test), ["id", "eom"]]
        output = output.merge(keys, on=["id", "eom"], how="inner", validate="one_to_one")
    output = output.sort_values(["eom", "id"], kind="stable").reset_index(drop=True)
    output["id"] = output.id.astype(np.int64)
    if output.empty or output.isna().any().any() or not np.isfinite(output.w).all():
        raise RuntimeError("Invalid CTF output")
    return output[["id", "eom", "w"]]


def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame:
    np.random.seed(SEED)
    panel, names = prepare_panel(chars, features)
    trained = train_models(panel, names)
    return build_portfolios(panel, daily_ret, trained)