// horizon_lambda.cu — ISV-driven per-horizon gradient scaler AND closed-form Kendall σ.
//
// Single source of truth for the two per-horizon controllers:
//
//   1) lambda[h] ∈ [LAMBDA_FLOOR, LAMBDA_CEILING] — multiplier on the
//      trunk grad_h contribution in heads_grn_bwd. Boosts horizons
//      that the model is currently failing to learn (per their
//      unweighted BCE EMA).
//
//   2) log_sigma_h[h] — Kendall σ for the BCE forward kernel. Computed
//      in closed form from the same loss_ema signal (Kendall's
//      gradient equilibrium ⟹ σ_h = sqrt(mean_bce_h)). Replaces
//      Adam-learned σ, eliminating the m/√v normalization conflict
//      (pearl_adam_normalizes_loss_weights).
//
// ISV anchors per `pearl_controller_anchors_isv_driven`: every anchor/
// target/cap derives from a tracked signal, never from hardcoded
// constants outside the bootstrap epsilons:
//
//   loss_ema[h]          — EMA(unweighted BCE per horizon)
//   z_max_ema            — EMA(max |z_h|) across horizons → drives
//                          adaptive Z_SCALE so the OBSERVED-max-z
//                          horizon maps exactly to LAMBDA_CEILING.
//                          With the prior hardcoded Z_SCALE=0.5 the
//                          controller rarely engaged on real data
//                          (max observed lambda ~1.04); adaptive
//                          Z_SCALE uses the full envelope.
//
// Sentinel bootstrap (`prev <= 0`) per
// `pearl_first_observation_bootstrap.md`; permanent floor
// (`max(real, floor)`) per `pearl_blend_formulas_must_have_permanent_floor.md`;
// asymmetric clamp (floor-only on σ, floor + ceiling on λ) per
// `pearl_audit_unboundedness_for_implicit_asymmetry.md`. Z-score
// normalization per `pearl_zscore_normalization_for_magnitude_asymmetric_signals.md`.

#define N_HORIZONS_LAMBDA 3
#define ALPHA_FIXED       0.1f
#define LAMBDA_FLOOR      1.0f
#define LAMBDA_CEILING    2.0f
#define LOG_SIGMA_FLOOR  (-0.6931472f)   // log(0.5) — σ never collapses below 0.5
#define Z_MAX_FLOOR       0.1f           // guards Z_SCALE_ISV when z_max_ema is tiny

extern "C" __global__ void horizon_ema_and_lambda(
    const float* __restrict__ loss_per_horizon,  // [3] — current step UNWEIGHTED BCE
    float* __restrict__ loss_ema,                 // [3] — EMA state (read + write)
    float* __restrict__ z_max_ema,                // [1] — EMA(max|z|) state (read + write)
    float* __restrict__ lambda,                   // [3] — output: λ_h
    float* __restrict__ log_sigma_h               // [3] — output: closed-form Kendall σ
) {
    if (threadIdx.x != 0 || blockIdx.x != 0) return;

    // -------- 1) EMA update on loss_per_horizon, with first-obs bootstrap.
    float new_ema[N_HORIZONS_LAMBDA];
    float sum = 0.0f;
    #pragma unroll
    for (int h = 0; h < N_HORIZONS_LAMBDA; ++h) {
        const float cur = loss_per_horizon[h];
        const float prev = loss_ema[h];
        new_ema[h] = (prev <= 0.0f) ? cur : (prev + ALPHA_FIXED * (cur - prev));
        loss_ema[h] = new_ema[h];
        sum += new_ema[h];
    }
    const float mean = sum / (float)N_HORIZONS_LAMBDA;

    // -------- 2) Z-score normalize across horizons.
    float ssq = 0.0f;
    #pragma unroll
    for (int h = 0; h < N_HORIZONS_LAMBDA; ++h) {
        const float diff = new_ema[h] - mean;
        ssq += diff * diff;
    }
    const float std = sqrtf(ssq / (float)N_HORIZONS_LAMBDA + 1e-12f);
    const float inv_std = (std > 1e-6f) ? (1.0f / std) : 0.0f;

    float z[N_HORIZONS_LAMBDA];
    float z_max = 0.0f;
    #pragma unroll
    for (int h = 0; h < N_HORIZONS_LAMBDA; ++h) {
        z[h] = (new_ema[h] - mean) * inv_std;
        const float az = fabsf(z[h]);
        if (az > z_max) z_max = az;
    }

    // -------- 3) z_max_ema: tracks the typical spread, drives adaptive Z_SCALE.
    //         Sentinel = 0 ⇒ first-obs bootstrap. Fixed α matches loss_ema's
    //         smoothing horizon.
    const float prev_zmax = z_max_ema[0];
    const float new_zmax = (prev_zmax <= 0.0f)
        ? z_max
        : (prev_zmax + ALPHA_FIXED * (z_max - prev_zmax));
    z_max_ema[0] = new_zmax;

    // -------- 4) Adaptive Z_SCALE: maps the historical-max-z to LAMBDA_CEILING.
    //         When z_max_ema is large (spread regime), Z_SCALE_ISV shrinks
    //         to avoid saturating early; when small (tight regime), it
    //         grows to amplify the controller into engagement.
    const float z_anchor = fmaxf(new_zmax, Z_MAX_FLOOR);
    const float z_scale_isv = (LAMBDA_CEILING - LAMBDA_FLOOR) / z_anchor;

    // -------- 5) Per-horizon outputs.
    #pragma unroll
    for (int h = 0; h < N_HORIZONS_LAMBDA; ++h) {
        // λ_h: boost-only, asymmetric clamped.
        const float raw_lambda = 1.0f + z_scale_isv * z[h];
        lambda[h] = fmaxf(LAMBDA_FLOOR, fminf(LAMBDA_CEILING, raw_lambda));

        // log σ_h: Kendall equilibrium from the same loss_ema signal.
        //   ∂L/∂log σ = 0  ⟹  σ_h² = mean_bce_h  ⟹  log σ_h = ½·log(loss_ema_h).
        // Floor at log(0.5) so σ can't collapse to 0 (which would explode w_h).
        const float ema_h = fmaxf(new_ema[h], 1e-12f);
        const float raw_log_sigma = 0.5f * logf(ema_h);
        log_sigma_h[h] = fmaxf(LOG_SIGMA_FLOOR, raw_log_sigma);
    }
}