// crates/ml/src/cuda_pipeline/alpha_grad_compute_kernel.cu
//
// SP14 Layer B Task B.4 (2026-05-05): Earned Gradient Flow consumer kernel —
// the heart of the EGF pearl. Single-thread state-machine kernel that reads
// driver signals from the global ISV bus, runs Schmitt-trigger Gate 1,
// computes adaptive k_aux/k_q/β from Welford variance EMAs, evaluates two
// sigmoids, multiplies with a host-supplied warmup gate, applies a β-rate-
// limiter on the result, and writes 7 outputs back to ISV.
//
// ── Single-thread launch ────────────────────────────────────────────────
// The body runs only on (threadIdx.x == 0, blockIdx.x == 0) — there is no
// parallelism to exploit (this is a state machine consuming O(1) inputs
// and producing O(1) outputs). Caller MUST still launch with
// blockDim ≥ (1,1,1) and gridDim ≥ (1,1,1); a 32-thread block is fine
// because the early-return masks all but lane 0. No shared memory, no
// reductions, no atomicAdd (per `feedback_no_atomicadd.md`).
//
// ── Driver inputs (read-only on this kernel) ────────────────────────────
//   ISV[TARGET_DIR_ACC = 372]            target directional accuracy (SP13)
//   ISV[AUX_DIR_ACC_SHORT_EMA = 373]     fast EMA of aux dir-acc (SP13)
//   ISV[AUX_DIR_ACC_LONG_EMA = 374]      slow EMA of aux dir-acc (SP13) — used
//                                        as the running mean for the var_aux
//                                        Welford EMA (SP14 B.11 var_aux producer)
//   ISV[Q_DISAGREEMENT_SHORT_EMA = 383]  fast EMA of Q-aux disagreement (SP14 B.3)
//   ISV[VAR_AUX = 388]                   Welford variance EMA of aux dir-acc
//                                        (SP14 B.11: this kernel now produces
//                                        the var_aux update — was a ghost
//                                        feature in B.4, see "RESOLVED" block
//                                        below)
//   ISV[VAR_Q   = 389]                   Welford variance EMA of Q disagreement (SP14 B.3)
//   ISV[GATE1_OPEN = 391]                persistent Schmitt state (read+written)
//   ISV[VAR_ALPHA  = 390]                Welford variance EMA of α_grad_raw (read+written)
//   ISV[ALPHA_SMOOTHED = 393]            α_grad_smoothed previous value
//
// ── Outputs (written by this kernel) ────────────────────────────────────
//   ISV[K_AUX_ADAPTIVE = 385]    adaptive Gate-1 sigmoid steepness
//   ISV[K_Q_ADAPTIVE   = 386]    adaptive Gate-2 sigmoid steepness
//   ISV[BETA_ADAPTIVE  = 387]    adaptive β rate-limiter coefficient
//   ISV[VAR_AUX        = 388]    updated aux_dir_acc variance EMA (SP14 B.11)
//   ISV[VAR_ALPHA      = 390]    updated α_grad_raw variance EMA
//   ISV[GATE1_OPEN     = 391]    updated Schmitt state (0 = closed, 1 = open)
//   ISV[ALPHA_RAW      = 392]    α_grad_raw = gate1 × gate2 (warmup_gate removed
//   ISV[ALPHA_SMOOTHED = 393]    α_grad_smoothed = β × prev + (1-β) × raw
//
// ── Slot indices: must match `sp14_isv_slots.rs` ────────────────────────
// The kernel hardcodes ISV slot indices via `const int` locals below.
// They MUST match the constants declared in
// `crates/ml/src/cuda_pipeline/sp14_isv_slots.rs` (and 372/373 from
// `sp13_isv_slots.rs`). Any edit to the .rs constants requires a matching
// edit here. The plan originally documented 381/383/384/385/386/387/388/
// 389/390/391 for SP14 slots — the actual values are +2 because SP13
// closeout added HOLD_RATE_TARGET=381 + HOLD_RATE_OBSERVED_EMA=382 after
// the plan was written. See `sp14_isv_slots.rs` lines 17-20.
//
// ── Schmitt-trigger hysteresis intentional discontinuity ────────────────
// At the moment Gate 1 transitions (closed → open or open → closed), the
// sigmoid argument switches between `(aux_short - threshold_open)` and
// `(aux_short - threshold_close)`, producing a small DISCONTINUITY in
// α_grad_raw. This is intentional, not a bug — Schmitt hysteresis
// requires asymmetric thresholds between the rising and falling edges to
// suppress chatter in noisy signals. The discontinuity is smoothed
// downstream by the β rate-limiter, and consumers downstream see the
// β-smoothed signal (α_grad_smoothed = ISV[393]), so the discontinuity
// is invisible to the gradient flow it modulates.
//
// ── Bounded modifier per pearl_bounded_modifier_outputs_require_… ───────
// α_grad_raw = sigmoid(·) × sigmoid(·) is structurally
// bounded to [0, 1] — the two sigmoids guarantee [0, 1] each, the warmup
// gate is supplied in [0, 1] by the host. No runtime clamp needed. Per
// `pearl_bounded_modifier_outputs_require_structural_activation.md`.
//
// ── Numerical stability: sigmoid arg clipping ───────────────────────────
// Sigmoid arguments are clipped to [-30, 30] before `__expf`. At |arg| >
// 30, fp32 sigmoid saturates to 0 or 1 with bit-equal output, so the
// clip is a precision-neutral overflow guard (`__expf(30) ≈ 1.07e13`,
// well within fp32 range; `__expf(38) ≈ 3.18e16` and `__expf(89)`
// overflows). The clip prevents NaN propagation from pathological
// k_aux × (aux_short - threshold) products if k_aux is corrupted.
//
// ── RESOLVED: var_aux producer wired in B.11 (2026-05-05) ───────────────
// B.4 originally READ ISV[388] (AUX_DIR_ACC_VARIANCE_EMA) but did NOT
// write it, so var_aux stayed at sentinel 0.0 and k_aux degenerated to
// the constant K_BASE_AUX. B.11 closes the gap inline: the kernel now
// also computes a Welford-style variance EMA on `aux_dir_acc_short` with
// `aux_dir_acc_long` as the running mean (mirrors how var_q in
// `q_disagreement_update_kernel` uses the fast EMA `new_short` as its
// mean reference). The compute lives next to the existing var_alpha
// update (same Welford pattern, same `alpha_var` blend coefficient) so
// the kernel stays a single state-machine launch with no companion
// kernel needed. Per `feedback_no_partial_refactor.md`, the kernel
// signature is unchanged; we add reads of ISV[374] and a write to
// ISV[388] inside the existing dispatch.
//
// var_q (slot 389) IS written — by `q_disagreement_update_kernel` (B.3),
// which lands the Welford variance EMA in the same launch as the mean
// EMAs. With var_aux now produced too, both adaptive k_aux and k_q are
// fully functional from B.11 onward.

extern "C" __global__
void alpha_grad_compute_kernel(
    /* Global ISV bus. Slots 385, 386, 387, 388, 390, 391, 392, 393 are
     * written; slots 372, 373, 374, 383, 388, 389 are read (slot 388 is
     * read-modify-write per the Welford pattern; slot 374 added in B.11
     * as the var_aux mean reference). All other slots untouched.
     */
    float* __restrict__ isv,
    /* EMA decay coefficient for the α_grad_raw variance update (typical
     * 0.05; matches the `alpha_var` argument to
     * `q_disagreement_update_kernel`).
     *
     * SP14 L1+L2 follow-up (2026-05-05): the original `warmup_gate`
     * parameter was deleted per `feedback_isv_for_adaptive_bounds`. The
     * variance-driven k_aux/k_q already provide warmup behavior (high
     * variance during cold-start collapses k → K_MIN → flat sigmoid →
     * gate ≈ 0.5; once EMAs settle, k → K_BASE → sharp sigmoid → gates
     * respond correctly). Adding a separate hardcoded step-counter
     * warmup multiplier was double-counting and tuning-driven.
     */
    float alpha_var)
{
    if (threadIdx.x != 0 || blockIdx.x != 0) return;

    // ── ISV slot indices (must match sp14_isv_slots.rs and sp13_isv_slots.rs) ──
    const int TARGET_DIR_ACC           = 372;
    const int AUX_DIR_ACC_SHORT_EMA    = 373;
    const int AUX_DIR_ACC_LONG_EMA     = 374;  // SP14 B.11: var_aux mean reference
    const int Q_DISAGREEMENT_SHORT_EMA = 383;
    const int K_AUX_ADAPTIVE           = 385;
    const int K_Q_ADAPTIVE             = 386;
    const int BETA_ADAPTIVE            = 387;
    const int VAR_AUX                  = 388;
    const int VAR_Q                    = 389;
    const int VAR_ALPHA                = 390;
    const int GATE1_OPEN               = 391;
    const int ALPHA_RAW                = 392;
    const int ALPHA_SMOOTHED           = 393;

    // ── Structural constants (must match sp14_isv_slots.rs) ────────────
    const float K_BASE_AUX             = 20.0f;
    const float K_BASE_Q               = 15.0f;
    const float K_MIN                  = 1.0f;
    const float VARIANCE_REF_AUX       = 0.01f;
    const float VARIANCE_REF_Q         = 0.05f;
    const float VARIANCE_REF_ALPHA     = 0.01f;
    const float BETA_BASE              = 0.5f;
    const float BETA_MAX               = 0.95f;
    const float SCHMITT_BAND           = 0.03f;
    const float Q_DISAGREEMENT_BASELINE = 0.5f;

    // ── Read drivers ───────────────────────────────────────────────────
    const float target              = isv[TARGET_DIR_ACC];
    const float aux_short           = isv[AUX_DIR_ACC_SHORT_EMA];
    const float aux_long            = isv[AUX_DIR_ACC_LONG_EMA];   // SP14 B.11
    const float q_dis               = isv[Q_DISAGREEMENT_SHORT_EMA];
    const float var_aux_prev        = isv[VAR_AUX];                 // SP14 B.11
    const float var_q               = isv[VAR_Q];
    const float var_alpha_prev      = isv[VAR_ALPHA];
    const float gate1_state_prev    = isv[GATE1_OPEN];
    const float alpha_smoothed_prev = isv[ALPHA_SMOOTHED];

    // ── Welford-style variance EMA on aux_dir_acc_short (SP14 B.11) ────
    // Reference is `aux_long` (slow EMA — best estimate of the running
    // mean). Under stable aux dir-acc the diff stays small and var_aux
    // decays toward 0; under chatter the diff grows and var_aux rises,
    // softening the Gate 1 sigmoid (k_aux below) so noisy aux signals
    // don't whipsaw the EGF gate. Uses the same `alpha_var` blend
    // coefficient as the var_alpha update further down — single shared
    // smoothing rate per Pearl-D Wiener-α convention. Closes the
    // ghost-feature gap left by B.4 (slot 388 had no producer; pre-B.11
    // the adaptive k_aux mechanism was degenerate at K_BASE_AUX = 20).
    const float diff_aux  = aux_short - aux_long;
    const float var_aux   = alpha_var * (diff_aux * diff_aux)
                          + (1.0f - alpha_var) * var_aux_prev;

    // ── Adaptive sigmoid steepness (B.2.5) ─────────────────────────────
    // Higher variance → smaller k → flatter sigmoid (more diffusion at
    // the edge). Floor at K_MIN prevents the sigmoid from collapsing to
    // a flat 0.5 line under unbounded noise. `var_aux` is the FRESHLY-
    // updated value above (this step's contribution included), so k_aux
    // tracks current-step noise rather than lagging by one launch.
    const float k_aux = fmaxf(K_BASE_AUX / (1.0f + var_aux / VARIANCE_REF_AUX), K_MIN);
    const float k_q   = fmaxf(K_BASE_Q   / (1.0f + var_q   / VARIANCE_REF_Q),   K_MIN);

    // ── Schmitt-trigger Gate 1 state update (B.2.4) ────────────────────
    // Open at target+0.03 (rising edge); close at target-0.03 (falling
    // edge). Persistent across calls via ISV[GATE1_OPEN].
    const float threshold_open  = target + SCHMITT_BAND;
    const float threshold_close = target - SCHMITT_BAND;
    float gate1_state_new = gate1_state_prev;
    if (gate1_state_prev < 0.5f && aux_short > threshold_open) {
        gate1_state_new = 1.0f;
    } else if (gate1_state_prev >= 0.5f && aux_short < threshold_close) {
        gate1_state_new = 0.0f;
    }

    // ── Sigmoid arg for Gate 1 (uses appropriate threshold per state) ──
    // Open state: distance from CLOSE threshold (must drop further to
    //   flip); the kernel evaluates the soft sigmoid above the close
    //   threshold so α stays high while the gate is open.
    // Closed state: distance from OPEN threshold (must rise further to
    //   flip); the sigmoid stays near 0 below the open threshold.
    // The intentional discontinuity at the transition is documented in
    // the header comment block above.
    const float gate1_arg_raw = (gate1_state_new >= 0.5f)
        ? (aux_short - threshold_close)
        : (aux_short - threshold_open);
    const float gate1_arg = fmaxf(-30.0f, fminf(30.0f, k_aux * gate1_arg_raw));
    const float gate1     = 1.0f / (1.0f + __expf(-gate1_arg));

    // ── Sigmoid for Gate 2 (B.2.3) ─────────────────────────────────────
    // q_dis above the analytic random-alignment baseline (0.5) is the
    // "Q-head agrees more with aux than chance" signal; sigmoid maps it
    // into [0, 1]. No Schmitt on Gate 2 — q_dis is already an EMA so
    // its high-frequency noise is filtered by upstream smoothing.
    const float gate2_arg_raw = q_dis - Q_DISAGREEMENT_BASELINE;
    const float gate2_arg     = fmaxf(-30.0f, fminf(30.0f, k_q * gate2_arg_raw));
    const float gate2         = 1.0f / (1.0f + __expf(-gate2_arg));

    // ── α_grad_raw (B.2.1) ─────────────────────────────────────────────
    // gate1 ∈ [0, 1], gate2 ∈ [0, 1], warmup_gate ∈ [0, 1] → product is
    // structurally bounded to [0, 1]. No runtime clamp per
    // `pearl_bounded_modifier_outputs_require_structural_activation.md`.
    // SP14 L1+L2 follow-up: warmup_gate multiplier removed. Variance-driven
    // k_aux/k_q (from var_aux, var_q EMAs) already provide warmup via the
    // sigmoid steepness collapse during cold-start (high variance → flat
    // sigmoid → gate ≈ 0.5 regardless of input). Per
    // `feedback_isv_for_adaptive_bounds`, adaptive bounds live in ISV
    // (here: variance EMAs); only Invariant-1 numerical-stability anchors
    // are constants. The separate hardcoded WARMUP_STEPS_FALLBACK ramp was
    // tuned, not signal-driven.
    const float alpha_raw = gate1 * gate2;

    // ── Welford-style variance EMA on α_grad_raw ───────────────────────
    // Reference is α_smoothed (the rate-limited tracking estimate); under
    // stable α_raw the diff stays small and var_alpha decays toward 0.
    // Under chatter (Schmitt hysteresis ringing, sigmoid edge oscillation)
    // diff grows and var_alpha rises — driving β toward β_max via the
    // adaptive-β formula below.
    const float diff_alpha    = alpha_raw - alpha_smoothed_prev;
    const float var_alpha_new = alpha_var * (diff_alpha * diff_alpha)
                              + (1.0f - alpha_var) * var_alpha_prev;

    // ── Adaptive β (B.2.8) ─────────────────────────────────────────────
    // Floor BETA_BASE prevents β collapsing below light-smoothing under
    // very-stable α; ceiling BETA_MAX prevents lockup under pathological
    // noise. The variance ratio drives the slide between the two.
    const float beta = fmaxf(BETA_BASE,
                       fminf(BETA_MAX,
                             BETA_BASE + var_alpha_new / VARIANCE_REF_ALPHA));

    // ── α_grad_smoothed (rate-limited) ─────────────────────────────────
    // β = 0 → no smoothing (α_smoothed = α_raw); β = 1 → frozen
    // (α_smoothed = α_smoothed_prev). The kernel writes the smoothed
    // value to ISV[393] for downstream consumers (the backward path's
    // gradient scaler in B.7+).
    const float alpha_smoothed = beta * alpha_smoothed_prev + (1.0f - beta) * alpha_raw;

    // ── Write back ─────────────────────────────────────────────────────
    isv[K_AUX_ADAPTIVE]  = k_aux;
    isv[K_Q_ADAPTIVE]    = k_q;
    isv[BETA_ADAPTIVE]   = beta;
    isv[VAR_AUX]         = var_aux;          // SP14 B.11: ghost-feature gap closed
    isv[VAR_ALPHA]       = var_alpha_new;
    isv[GATE1_OPEN]      = gate1_state_new;
    isv[ALPHA_RAW]       = alpha_raw;
    isv[ALPHA_SMOOTHED]  = alpha_smoothed;
}
