// crates/ml/src/cuda_pipeline/alpha_kill_criteria.cu
//
// Phase E.0 kill-criteria diagnostic producer kernel (Task 9, 2026-05-15).
//
// Writes 4 raw observations into a contiguous scratch_out[0..4] block.
// Downstream `apply_pearls_ad_kernel` (n_slots=4) chained on the same stream
// applies Pearl A first-observation bootstrap + Pearl D Wiener-optimal α
// smoothing into ISV slots 539..542 — see alpha_isv_slots.rs for the
// slot reservation and project_phase_e0_foundation_close memory for the
// baseline values that drive slot 541's denominator.
//
// Observations (written to scratch_out, in this exact order so the
// chained apply_pearls launcher can read them via a single
// scratch_idx_base + n_slots=4 stride):
//
//   scratch_out[0]  →  ISV[539]  Q_SPREAD_EMA_INDEX
//                       mean_over_batch( std(Q, axis=action) / |mean(Q)| )
//                       — collapses to ~0 when the network is in degenerate
//                       constant-Q mode (kill criterion: must be ≥ 0.05)
//
//   scratch_out[1]  →  ISV[540]  ACTION_ENTROPY_EMA_INDEX
//                       H(empirical action distribution over recent rollout)
//                       in nats; kill criterion ≥ 0.5 · ln(9) ≈ 1.10
//
//   scratch_out[2]  →  ISV[541]  RETURN_VS_RANDOM_EMA_INDEX
//                       (rollout_R_mean − ISV[547]) / max(ISV[548], 1e-6)
//                       — signed σ-distance from the Task 7c random-uniform
//                       baseline; kill criterion ≥ 0 (i.e. policy ≥ baseline)
//
//   scratch_out[3]  →  ISV[542]  EARLY_Q_MOVEMENT_EMA_INDEX
//                       |q_early_norm − q_init_norm| / max(|q_init_norm|, 1e-6)
//                       — relative L2-norm drift of early-step Q from init;
//                       kill criterion ≥ 0.01 (some learning happened)
//
// All four are scalars; the kernel is single-block, single-thread, cheap
// arithmetic per `pearl_fused_per_group_statistics_oracle` (one fused
// kernel, no atomicAdd, no host branches).
//
// Inputs:
//   q_values         device [batch, n_actions] — Q from the most-recent
//                    forward pass; reduced in-kernel.
//   action_counts    device [n_actions] — empirical action histogram from
//                    the rollout (host writes via mapped-pinned counter
//                    array or a producer kernel — either way, device-
//                    visible by the time this launches).
//   scalar_inputs    mapped-pinned [3] — host writes between rollouts:
//                      [0] rollout_reward_mean
//                      [1] q_init_norm  (captured once, post-init)
//                      [2] q_early_step_norm (re-captured each rollout)
//   isv              device ISV bus (read-only here); we read slots
//                    `random_baseline_mean_slot` (547) and
//                    `random_baseline_std_slot` (548).
//
// Output:
//   scratch_out      mapped-pinned (or device) [4] — consumed by the
//                    chained apply_pearls_ad_kernel launch.
//
// __threadfence_system() guarantees scratch_out writes are visible
// before the Pearls applicator reads them.

#include <cuda_runtime.h>

extern "C" __global__ void alpha_kill_criteria_compute_kernel(
    const float* __restrict__ q_values,
    const int*   __restrict__ action_counts,
    const float* __restrict__ scalar_inputs,
    const float* __restrict__ isv,
    int          random_baseline_mean_slot,
    int          random_baseline_std_slot,
    int          batch,
    int          n_actions,
    float*       __restrict__ scratch_out
) {
    // Single-block, single-thread: 4 scalar reductions of small batches.
    // No cross-thread parallelism opportunity is meaningful here.
    if (blockIdx.x != 0 || threadIdx.x != 0) return;

    // ---- (1) Q-spread: mean_b( std_a(Q[b, :]) / |mean_a(Q[b, :])| ) ----
    //
    // Per-row two-pass reduction. n_actions = 9 in Phase E so the inner
    // loops are tiny; batch ≤ ~256 typically. Single thread is fast enough.
    float q_spread = 0.0f;
    if (batch > 0 && n_actions > 0) {
        float q_spread_sum = 0.0f;
        const float n_a_inv = 1.0f / (float)n_actions;
        for (int b = 0; b < batch; ++b) {
            const int row = b * n_actions;
            float mean = 0.0f;
            for (int a = 0; a < n_actions; ++a) {
                mean += q_values[row + a];
            }
            mean *= n_a_inv;
            float var = 0.0f;
            for (int a = 0; a < n_actions; ++a) {
                const float d = q_values[row + a] - mean;
                var += d * d;
            }
            var *= n_a_inv;
            const float std_v = sqrtf(var);
            // |mean| + ε guards against 0/0 when Q is centered exactly at 0.
            // ε = 1e-9 matches EPS_DIV from sp4_wiener_ema.rs (Adam-numerical
            // category, not a tuning knob).
            q_spread_sum += std_v / (fabsf(mean) + 1.0e-9f);
        }
        q_spread = q_spread_sum / (float)batch;
    }

    // ---- (2) Action entropy: H(empirical action dist), in nats ----
    int total = 0;
    for (int a = 0; a < n_actions; ++a) {
        total += action_counts[a];
    }
    float entropy = 0.0f;
    if (total > 0) {
        const float inv_total = 1.0f / (float)total;
        for (int a = 0; a < n_actions; ++a) {
            const int c = action_counts[a];
            if (c > 0) {
                const float p = (float)c * inv_total;
                entropy -= p * logf(p);
            }
        }
    }

    // ---- (3) Return-vs-random: σ-distance from baseline ----
    //
    // Baseline mean/std live in ISV slots 547/548 (TrainingPersist anchors
    // populated once at Task 7c run; never reset). Their typical values are
    // negative-ish (mean ≈ -5185, std ≈ 4953 for the H=600 baseline) so the
    // signed division is meaningful — positive = policy outperforms random.
    const float rollout_R       = scalar_inputs[0];
    const float random_R_mean   = isv[random_baseline_mean_slot];
    const float random_R_std    = isv[random_baseline_std_slot];
    const float return_vs_random =
        (rollout_R - random_R_mean) / fmaxf(random_R_std, 1.0e-6f);

    // ---- (4) Early-Q-movement: |Δ Q-norm| / |Q-init-norm| ----
    //
    // q_init_norm is captured once at training start (post-weight-init,
    // pre-first-update) and held constant for the run. q_early_step_norm
    // is recomputed each rollout. Both are L2 norms of the Q-output for
    // a fixed canary state batch held by the host. The denominator's
    // ε guards against q_init_norm ≈ 0 in adversarial init scenarios.
    const float q_init       = scalar_inputs[1];
    const float q_early      = scalar_inputs[2];
    const float early_movement =
        fabsf(q_early - q_init) / fmaxf(fabsf(q_init), 1.0e-6f);

    scratch_out[0] = q_spread;
    scratch_out[1] = entropy;
    scratch_out[2] = return_vs_random;
    scratch_out[3] = early_movement;

    // Per the val_sharpe_delta_compute_kernel pattern: fence so the
    // chained apply_pearls_ad_kernel reads what we wrote here, not stale
    // values. Same fence the consumer kernels expect upstream.
    __threadfence_system();
}
