// crates/ml/src/cuda_pipeline/pnl_aggregation_kernel.cu
//
// SP5 Layer D Task D1 (rewrite, 2026-05-02): PnL aggregation producer kernel.
//
// **Rewrite supersedes the original D1 (commit `5ee795f14`).** The original
// kernel was authored against a brief that assumed per-trade event arrays
// (`trade_pnls`, `trade_durations`, `trade_directions`) exist in production.
// They do not. Production exposes per-bar `step_returns` + `done_flags`
// arrays plus already-aggregated summary scalars (`sum_returns`,
// `sum_sq_returns`, `n_trades`) on `TradeStats` (see
// `gpu_experience_collector.rs:60-86`). The host-side oracle —
// `compute_epoch_financials` in `crates/ml/src/trainers/dqn/financials.rs`
// — drives every output of this kernel from those exact inputs. The previous
// D4 agent found the per-trade vs per-bar mismatch as a structural blocker
// before any code changes; this rewrite reconciles the kernel with the host
// oracle bit-for-bit so D4 atomic wiring becomes straightforward.
//
// ──────────────────────────────────────────────────────────────────────────
// Output formulas — verbatim mirror of `compute_epoch_financials`
// ──────────────────────────────────────────────────────────────────────────
//
//   pnl_total      Log-space compounded total return over the rollout.
//                  `exp(Σ ln(max(1.0 + step_returns[i], 1e-10))) − 1`,
//                  matching financials.rs:80-97 (line-by-line). Log-space sum
//                  avoids f64 overflow on long rollouts; the 1e-10 floor
//                  prevents `ln(0)` when a single bar realises a −100% step.
//                  Empty input ⇒ 0.0.
//
//   pnl_mean       Per-trade mean = `sum_returns / n_trades`. Matches
//                  financials.rs:122. **Per-trade, not per-bar** — the
//                  experience collector accumulates `sum_returns` /
//                  `sum_sq_returns` from per-trade exit events, not per-bar
//                  step returns, so dividing by `num_bars` here would change
//                  the metric semantics. `n_trades == 0` ⇒ 0.0 (financials
//                  guards via `total_trades > 1` for the variance/Sharpe
//                  branch; we mirror the same guard).
//
//   pnl_var        Per-trade population variance =
//                  `max(0, (sum_sq_returns / n_trades) − mean²)`. Matches
//                  financials.rs:123-124. Variance is non-negative by
//                  definition; the floor handles catastrophic cancellation
//                  in the one-pass form. `n_trades <= 1` ⇒ 0.0 (financials
//                  branch on `total_trades > 1`).
//
//   pnl_max_dd     Per-bar equity-curve max drawdown over the **last
//                  MAX_DD_WINDOW=10_000 bars** with episode-boundary resets.
//                  Matches financials.rs:163-194 line-by-line:
//                    1. `equity = peak = initial_capital`
//                    2. For each bar i in [num_bars − 10_000, num_bars):
//                         `equity *= 1.0 + step_returns[i]`
//                         `if equity > peak: peak = equity`
//                         `if done_flags[i] > 0.5:        // post-update`
//                             `equity = peak = initial_capital`
//                         `dd = (peak > 1e-10) ? (peak − equity) / peak : 0`
//                         `max_dd = max(max_dd, dd)`
//                    3. Final cap: `max_dd = min(max_dd, 1.0)`.
//                  The episode-boundary reset fires AFTER processing the
//                  done bar's return — the loss on the done bar belongs to
//                  the just-ending episode; the next bar starts fresh.
//                  Empty input ⇒ 0.0.
//
// ──────────────────────────────────────────────────────────────────────────
// Block layout
// ──────────────────────────────────────────────────────────────────────────
//
//   - Single block, 256 threads, grid_dim=(1,1,1) — matches every other SP5
//     producer kernel (`feedback_no_atomicadd.md`; in-block tree-reduce only).
//   - Pass 1 (parallel): each thread accumulates a strided slice of the
//     log-space sum `Σ ln(max(1+r, 1e-10))` across the windowed returns.
//     In-block tree-reduction in shared memory folds the per-thread partial
//     sums into one scalar.
//   - Pass 2 (serial, tid == 0): walk the last MAX_DD_WINDOW bars to compute
//     `max_dd`. The equity / peak pair is inherently sequential — peak is
//     a running max of a cumulative-product, with sequence-breaking resets
//     at episode boundaries — so a single-thread serial scan is correct by
//     construction. Cost is bounded at MAX_DD_WINDOW = 10⁴ iterations,
//     well under the per-step producer-cadence budget on the host smoke.
//   - Final writes by tid == 0 with one `__threadfence_system()` so the
//     mapped-pinned host_ptr / dev_ptr alias sees the values before the
//     `apply_pearls_ad_kernel` consumer reads them on the same stream.
//
// Producer-only — D1 is additive. The host-side aggregation in
// `training_loop.rs:4862-4916` continues to run unchanged. D4 (the atomic
// Layer D commit) wires D1+D2+D3 to call sites in a single migration per
// `feedback_no_partial_refactor.md`.
//
// References:
//   - `feedback_no_cpu_compute_strict.md` (host EMA / reduction → GPU)
//   - `feedback_no_atomicadd.md`
//   - `feedback_trust_code_not_docs.md` (oracle is the code, not the brief)
//   - `pearl_first_observation_bootstrap.md` (Pearl A sentinel chain)

extern "C" __global__ void pnl_aggregation_update(
    const float* __restrict__ step_returns,    /* [num_bars] per-bar log/arith returns */
    const float* __restrict__ done_flags,      /* [num_bars] 1.0 = end of episode */
    int          num_bars,
    int          n_trades,                      /* per-trade count for mean / variance */
    float        sum_returns,                   /* per-trade Σ returns (already reduced) */
    float        sum_sq_returns,                /* per-trade Σ returns² (already reduced) */
    float        initial_capital,               /* equity-curve starting capital */
    float*       __restrict__ scratch_buf,      /* producer_step_scratch_buf */
    int          pnl_total_idx,
    int          pnl_mean_idx,
    int          pnl_var_idx,
    int          pnl_max_dd_idx
) {
    /* Single-block contract — match SP5/SP4 producer kernel grid guard. */
    if (blockIdx.x != 0) return;

    extern __shared__ float smem[];           /* blockDim.x floats (log-growth bank) */
    const int tid   = (int)threadIdx.x;
    const int block = (int)blockDim.x;

    /* Empty-input fast path. Thread 0 zeros all four output slots so the
     * downstream apply_pearls chain reads a defined value every step. */
    if (num_bars <= 0) {
        if (tid == 0) {
            scratch_buf[pnl_total_idx]  = 0.0f;
            scratch_buf[pnl_mean_idx]   = 0.0f;
            scratch_buf[pnl_var_idx]    = 0.0f;
            scratch_buf[pnl_max_dd_idx] = 0.0f;
            __threadfence_system();
        }
        return;
    }

    /* ── Pass 1: log-space cumulative growth over all bars ─────────────────
     * `Σ ln(max(1 + step_returns[i], 1e-10))` — financials.rs:84-89. The
     * 1e-10 floor matches the host clamp ("a -100% bar gives a finite -23
     * log contribution rather than -inf"). The result is exponentiated into
     * `pnl_total = exp(Σ) - 1` after the reduction. */
    float local_log_growth = 0.0f;
    for (int i = tid; i < num_bars; i += block) {
        float g = 1.0f + step_returns[i];
        if (g < 1e-10f) g = 1e-10f;
        local_log_growth += logf(g);
    }
    smem[tid] = local_log_growth;
    __syncthreads();

    /* In-block tree reduction (no atomicAdd). */
    for (int s = block / 2; s > 0; s >>= 1) {
        if (tid < s) {
            smem[tid] += smem[tid + s];
        }
        __syncthreads();
    }

    /* Thread 0 holds the reduced log-growth. Compute total / mean / var
     * from the (already-reduced) per-trade sums. financials.rs:90-94 +
     * 122-124. Variance non-negativity floor at 0 matches the
     * `var.max(0.0)` host clamp. */
    float total  = 0.0f;
    float mean   = 0.0f;
    float var    = 0.0f;
    if (tid == 0) {
        const float log_growth = smem[0];
        if (isfinite(log_growth)) {
            total = expf(log_growth) - 1.0f;
        }
        if (n_trades > 0) {
            const float n = (float)n_trades;
            mean = sum_returns / n;
        }
        if (n_trades > 1) {
            const float n        = (float)n_trades;
            const float e_xsq    = sum_sq_returns / n;
            const float mean_sq  = mean * mean;
            float v              = e_xsq - mean_sq;
            if (v < 0.0f) v = 0.0f;
            var = v;
        }
    }

    /* ── Pass 2: max drawdown over last MAX_DD_WINDOW bars ────────────────
     * Verbatim port of financials.rs:163-194. equity / peak start at
     * initial_capital; both reset at done_flags[i] > 0.5 AFTER processing
     * the done bar's return. dd uses peak > 1e-10 guard. Final max_dd
     * capped at 1.0 (full liquidation upper bound). */
    constexpr int MAX_DD_WINDOW = 10000;
    float max_dd = 0.0f;
    if (tid == 0) {
        const int dd_start = num_bars > MAX_DD_WINDOW
            ? num_bars - MAX_DD_WINDOW
            : 0;
        float equity = initial_capital;
        float peak   = initial_capital;
        for (int i = dd_start; i < num_bars; i++) {
            equity *= 1.0f + step_returns[i];
            if (equity > peak) peak = equity;
            /* Reset AFTER processing the done bar's return — the loss on
             * the done bar belongs to the ending episode; the next bar
             * starts fresh. financials.rs:175-183. */
            if (done_flags[i] > 0.5f) {
                equity = initial_capital;
                peak   = initial_capital;
            }
            const float dd = (peak > 1e-10f)
                ? (peak - equity) / peak
                : 0.0f;
            if (dd > max_dd) max_dd = dd;
        }
        if (max_dd > 1.0f) max_dd = 1.0f;     /* financials.rs:194 */

        /* Single PCIe-visible barrier across all four writes. */
        scratch_buf[pnl_total_idx]  = total;
        scratch_buf[pnl_mean_idx]   = mean;
        scratch_buf[pnl_var_idx]    = var;
        scratch_buf[pnl_max_dd_idx] = max_dd;
        __threadfence_system();
    }
}
