// Per-window metrics reduction kernel.
// One block per window. Threads cooperate to reduce step_returns.
//
// Output per window [14 floats]:
//   [0] sharpe_ratio (annualized via configurable annualization_factor)
//   [1] total_pnl (multiplicatively compounded fractional return, e.g. 0.05 = +5%)
//   [2] max_drawdown (exact sequential scan, worst peak-to-trough fraction)
//   [3] sortino_ratio
//   [4] win_rate (winning_trades / total_trades, exact via boundary stitching)
//   [5] total_trades (exact position changes, boundary-aware)
//   [6] var_95 (5th-percentile return — Value at Risk at 95% confidence)
//   [7] cvar_95 (mean of returns below VaR — Expected Shortfall)
//   [8] calmar_ratio (annualized mean return / max drawdown)
//   [9] omega_ratio (sum of gains / sum of losses)
//  [10] buy_count  (exposure_idx 5-8 — Long25..Long100)
//  [11] sell_count (exposure_idx 0-3 — Short100..Short25)
//  [12] hold_count (exposure_idx 4 — Flat)
//  [13] unique_action_count (popcount of 9-bit mask — how many exposure levels 0..8 were used)

extern "C" __global__ void compute_backtest_metrics(
    const float* __restrict__ step_returns,   // [n_windows * max_len]
    const float* __restrict__ portfolio_state, // [n_windows * 8]
    const int*   __restrict__ window_lens,    // [n_windows]
    const int*   __restrict__ actions_history, // [n_windows * max_len] for trade counting
    float*       metrics_out,                 // [n_windows * 14]
    int          n_windows,
    int          max_len,
    float        annualization_factor          // sqrt(252) for daily
) {
    int w = blockIdx.x;
    if (w >= n_windows) return;

    int wlen = window_lens[w];
    int tid = threadIdx.x;
    int stride = blockDim.x;
    int base = w * max_len;

    // Shared memory layout (6 reduction arrays + sort scratch):
    //   [0 .. stride)        : s_sum
    //   [stride .. 2*stride) : s_sq_sum
    //   [2*stride..3*stride) : s_down_sq
    //   [3*stride..4*stride) : s_buys
    //   [4*stride..5*stride) : s_sells
    //   [5*stride..6*stride) : s_holds
    //   [6*stride .. 6*stride + 4096) : s_sorted  (bitonic sort scratch + boundary data)
    extern __shared__ float shmem[];
    float* s_sum     = shmem;                // [blockDim.x]
    float* s_sq_sum  = shmem + stride;       // [blockDim.x]
    float* s_down_sq = shmem + 2*stride;     // [blockDim.x] (downside deviation)
    // Action distribution: exposure 0-3=sell, 4=flat, 5-8=buy
    float* s_buys    = shmem + 3*stride;     // [blockDim.x]
    float* s_sells   = shmem + 4*stride;     // [blockDim.x]
    float* s_holds   = shmem + 5*stride;     // [blockDim.x]
    float* s_sorted  = shmem + 6*stride;     // [4096] for bitonic sort + boundary stitching

    // Pass 1: per-thread local accumulators
    float local_sum = 0.0f, local_sq = 0.0f, local_down = 0.0f;
    // Multiplicative compounding for equity product reduction (total_pnl)
    float local_cum = 1.0f;
    int local_buys = 0, local_sells = 0, local_holds = 0;
    int local_action_mask = 0;  // 9-bit bitmask: bit i set if exposure_idx i was seen

    // Boundary-aware trade tracking: each thread exports metadata about
    // trades at chunk boundaries. Thread 0 stitches these to produce
    // exact trade counts and win rates across the full window.
    int bnd_first_action = -1;     // -1 signals empty chunk
    int bnd_last_action = -1;
    float bnd_prefix_return = 0.0f;
    float bnd_suffix_return = 0.0f;
    int bnd_complete_trades = 0;
    int bnd_complete_wins = 0;
    int bnd_num_changes = 0;
    int bnd_cur_action = -1;
    float bnd_cur_return = 0.0f;

    // Process consecutive chunks (not strided) for correct drawdown tracking.
    // Strided access (bars 0, 256, 512...) skips intermediate peaks/troughs,
    // producing wrong drawdown estimates. Chunked access ensures each thread
    // sees consecutive bars within its chunk.
    int chunk_size = (wlen + stride - 1) / stride;
    int chunk_start = tid * chunk_size;
    int chunk_end = chunk_start + chunk_size;
    if (chunk_end > wlen) chunk_end = wlen;

    for (int i = chunk_start; i < chunk_end; i++) {
        float r = step_returns[base + i];
        // Sharpe/Sortino use arithmetic mean of per-bar returns
        local_sum += r;
        local_sq += r * r;
        if (r < 0.0f) local_down += r * r;

        // Equity product for total_pnl (drawdown computed exactly by thread 0 later)
        local_cum *= (1.0f + r);

        // Action distribution: decode exposure from factored action.
        // Factored action = exposure * (b1*b2) + order * b2 + urgency.
        // exposure_idx: 0-3=sell, 4=flat, 5-8=buy (for 9-action space)
        int act = actions_history[base + i];
        int exp_idx = act / (DQN_ORDER_ACTIONS * DQN_URGENCY_ACTIONS);
        if (exp_idx < DQN_NUM_ACTIONS / 2)       local_sells++;
        else if (exp_idx == DQN_NUM_ACTIONS / 2)  local_holds++;
        else                                      local_buys++;

        // Unique action tracking: set bit per exposure index (0..8)
        if (exp_idx >= 0 && exp_idx < DQN_NUM_ACTIONS)
            local_action_mask |= (1 << exp_idx);

        // Boundary-aware trade tracking
        if (i == chunk_start) {
            bnd_first_action = act;
            bnd_cur_action = act;
        }

        if (act != bnd_cur_action) {
            bnd_num_changes++;
            if (bnd_num_changes == 1) {
                // First change: everything accumulated so far is the prefix
                bnd_prefix_return = bnd_cur_return;
            } else {
                // Interior change: segment since last change is a complete trade
                bnd_complete_trades++;
                if (bnd_cur_return > 0.0f) bnd_complete_wins++;
            }
            bnd_cur_return = 0.0f;
            bnd_cur_action = act;
        }

        bnd_cur_return += r;
        bnd_last_action = act;
    }

    // Post-loop: assign suffix and handle no-change edge case
    bnd_suffix_return = bnd_cur_return;
    if (bnd_num_changes == 0) {
        // Entire chunk is one partial trade segment (prefix == suffix)
        bnd_prefix_return = bnd_cur_return;
    }

    // Store local values to shared memory for parallel reduction
    s_sum[tid] = local_sum;
    s_sq_sum[tid] = local_sq;
    s_down_sq[tid] = local_down;
    s_buys[tid] = (float)local_buys;
    s_sells[tid] = (float)local_sells;
    s_holds[tid] = (float)local_holds;
    __syncthreads();

    // Block-level parallel reduction for 6 arrays
    for (int s = stride / 2; s > 0; s >>= 1) {
        if (tid < s) {
            s_sum[tid] += s_sum[tid + s];
            s_sq_sum[tid] += s_sq_sum[tid + s];
            s_down_sq[tid] += s_down_sq[tid + s];
            s_buys[tid] += s_buys[tid + s];
            s_sells[tid] += s_sells[tid + s];
            s_holds[tid] += s_holds[tid + s];
        }
        __syncthreads();
    }

    // ── Equity product reduction (multiplicative) ──────────────────────────
    // Each thread computed local_cum = product((1+r_i)) for its consecutive
    // chunk. The full-window equity = product of all thread equities.
    // Threads with no data (chunk_start >= chunk_end) have local_cum = 1.0
    // (identity for multiplication), so they don't affect the product.
    // We use s_sorted[0..stride] for this reduction (safe: bitonic sort is later).
    s_sorted[tid] = local_cum;
    __syncthreads();

    for (int s = stride / 2; s > 0; s >>= 1) {
        if (tid < s) {
            s_sorted[tid] *= s_sorted[tid + s];  // multiplicative reduction
        }
        __syncthreads();
    }
    // s_sorted[0] now holds the full-window equity = product of all (1+r_i)

    // ── Export boundary trade data to shared memory ──────────────────────
    // Reuse s_sorted[stride..8*stride] for 7 boundary arrays (1792 floats).
    // This region is free between the equity product reduction (which only
    // uses s_sorted[0]) and the bitonic sort (which starts later).
    // Int values stored via __int_as_float for bit-exact round-trip.
    {
        float* s_bnd = &s_sorted[stride];
        s_bnd[0 * stride + tid] = __int_as_float(bnd_first_action);
        s_bnd[1 * stride + tid] = __int_as_float(bnd_last_action);
        s_bnd[2 * stride + tid] = bnd_prefix_return;
        s_bnd[3 * stride + tid] = bnd_suffix_return;
        s_bnd[4 * stride + tid] = __int_as_float(bnd_complete_trades);
        s_bnd[5 * stride + tid] = __int_as_float(bnd_complete_wins);
        s_bnd[6 * stride + tid] = __int_as_float(bnd_num_changes);
    }
    __syncthreads();

    // Thread 0 computes final metrics from fully reduced values
    if (tid == 0) {
        float n = (float)wlen;
        float mean = s_sum[0] / n;
        float var = s_sq_sum[0] / n - mean * mean;
        float std = sqrtf(fmaxf(var, 1e-10f));
        float down_std = sqrtf(fmaxf(s_down_sq[0] / n, 1e-10f));

        // ── Exact max drawdown via sequential scan ──────────────────────
        // The per-thread parallel reduction misses cross-chunk drawdowns
        // because each thread tracks equity from local_peak=1.0, unaware
        // of the global peak. Sequential scan over step_returns gives the
        // TRUE maximum peak-to-trough drawdown.
        // Cost: O(wlen) ≈ 10K iterations ≈ microseconds on GPU.
        float exact_equity = 1.0f;
        float exact_peak = 1.0f;
        float exact_max_dd = 0.0f;
        for (int i = 0; i < wlen; i++) {
            float r = step_returns[base + i];
            exact_equity *= (1.0f + r);
            exact_peak = fmaxf(exact_peak, exact_equity);
            float dd = (exact_peak - exact_equity) / fmaxf(exact_peak, 1e-8f);
            exact_max_dd = fmaxf(exact_max_dd, dd);
        }

        // ── Boundary stitching: exact trade count + win rate ─────────────
        // Sum interior complete trades from all threads, then stitch
        // boundary trades by walking the per-thread metadata left-to-right.
        float* s_bnd = &s_sorted[stride];
        int total_trades = 0;
        int total_wins = 0;

        // Sum interior trades (parallel-counted, always correct)
        for (int t = 0; t < stride; t++) {
            total_trades += __float_as_int(s_bnd[4 * stride + t]);
            total_wins   += __float_as_int(s_bnd[5 * stride + t]);
        }

        // Stitch boundary trades
        float open_return = 0.0f;
        int open_action = -1;

        for (int t = 0; t < stride; t++) {
            int fa = __float_as_int(s_bnd[0 * stride + t]);
            int la = __float_as_int(s_bnd[1 * stride + t]);
            float pr = s_bnd[2 * stride + t];
            float sr = s_bnd[3 * stride + t];
            int nc = __float_as_int(s_bnd[6 * stride + t]);

            if (fa < 0) continue;  // empty chunk

            // Step 1: connect this chunk's start to the open trade
            if (open_action < 0) {
                // First valid chunk — open the first trade (no closure yet)
                open_action = fa;
                open_return = 0.0f;
            } else if (fa != open_action) {
                // Action changed at chunk boundary — close the open trade
                total_trades++;
                if (open_return > 0.0f) total_wins++;
                open_action = fa;
                open_return = 0.0f;
            }

            // Step 2: process this chunk's contribution
            if (nc == 0) {
                // No action changes — entire chunk extends the open trade
                open_return += pr;
            } else {
                // Prefix extends and closes the open trade
                open_return += pr;
                total_trades++;
                if (open_return > 0.0f) total_wins++;
                // Interior trades already counted above
                // Suffix becomes the new open trade
                open_action = la;
                open_return = sr;
            }
        }

        // Close the final open trade (window end)
        if (open_action >= 0) {
            total_trades++;
            if (open_return > 0.0f) total_wins++;
        }

        int out_base = w * 14;
        metrics_out[out_base + 0] = (mean / std) * annualization_factor;  // Sharpe
        // Total return: (final_equity / initial_equity) - 1.0 as fraction
        metrics_out[out_base + 1] = s_sorted[0] - 1.0f;
        metrics_out[out_base + 2] = exact_max_dd;  // max drawdown (exact sequential scan)
        metrics_out[out_base + 3] = (mean / down_std) * annualization_factor;  // Sortino

        float total_trades_f = (float)total_trades;
        metrics_out[out_base + 4] = (total_trades_f > 0.0f) ? (float)total_wins / total_trades_f : 0.0f;
        metrics_out[out_base + 5] = total_trades_f;

        metrics_out[out_base + 10] = s_buys[0];
        metrics_out[out_base + 11] = s_sells[0];
        metrics_out[out_base + 12] = s_holds[0];
    }

    // ── Unique action count via bitmask OR-reduction ─────────────────────────
    // Reuses s_sorted memory (safe: equity reduction is complete, OR-reduction
    // completes before bitonic sort).
    // Each thread accumulated a 9-bit mask of observed exposure indices. OR-reduce
    // across the block, then popcount gives unique_action_count (1-9).
    {
        int* s_action_mask = (int*)s_sorted;  // alias — sequential use, no conflict
        s_action_mask[tid] = local_action_mask;
        __syncthreads();

        for (int s = stride / 2; s > 0; s >>= 1) {
            if (tid < s) {
                s_action_mask[tid] |= s_action_mask[tid + s];
            }
            __syncthreads();
        }

        if (tid == 0) {
            int out_base = w * 14;
            metrics_out[out_base + 13] = (float)__popc(s_action_mask[0]);
        }
    }

    // ── Extended metrics: VaR, CVaR, Calmar, Omega via bitonic sort ──────────
    //
    // All threads cooperate to load and sort up to 4096 step_returns for this
    // window into s_sorted (ascending order). Thread 0 then scans the sorted
    // array to derive the tail-risk metrics.

    int sort_len = wlen < 4096 ? wlen : 4096;

    // Load returns into sort scratch.
    // For windows > 4096 bars, use strided sampling to cover the full window
    // evenly instead of only the first 4096 bars (which biases VaR/CVaR
    // toward early-window volatility characteristics).
    int sample_step = (wlen + sort_len - 1) / sort_len;  // ceil(wlen / sort_len)
    for (int i = tid; i < sort_len; i += stride) {
        int src = i * sample_step;
        s_sorted[i] = (src < wlen) ? step_returns[base + src] : 0.0f;
    }

    // Compute next power-of-two for bitonic sort
    int padded_len = 1;
    while (padded_len < sort_len) padded_len <<= 1;

    // Pad with +inf so sentinel values sort to the end (ascending)
    for (int i = sort_len + tid; i < padded_len; i += stride) {
        s_sorted[i] = 1e30f;
    }
    __syncthreads();

    // Bitonic sort — ascending
    for (int k = 2; k <= padded_len; k <<= 1) {
        for (int j = k >> 1; j > 0; j >>= 1) {
            for (int i = tid; i < padded_len; i += stride) {
                int ixj = i ^ j;
                if (ixj > i) {
                    bool ascending = ((i & k) == 0);
                    if ((ascending  && s_sorted[i] > s_sorted[ixj]) ||
                        (!ascending && s_sorted[i] < s_sorted[ixj])) {
                        float tmp   = s_sorted[i];
                        s_sorted[i] = s_sorted[ixj];
                        s_sorted[ixj] = tmp;
                    }
                }
            }
            __syncthreads();
        }
    }

    // Thread 0 computes extended metrics from the sorted array
    if (tid == 0) {
        int out_base = w * 14;

        // VaR at 95% confidence = 5th-percentile of sorted returns
        int var_idx = (int)(0.05f * (float)sort_len);
        if (var_idx < 0) var_idx = 0;
        if (var_idx >= sort_len) var_idx = sort_len - 1;
        float var_95 = s_sorted[var_idx];

        // CVaR (Expected Shortfall): mean of returns strictly below VaR index
        int cvar_count = (var_idx > 0) ? var_idx : 1;
        float cvar_sum = 0.0f;
        for (int i = 0; i < cvar_count; i++) {
            cvar_sum += s_sorted[i];
        }
        float cvar_95 = cvar_sum / (float)cvar_count;

        // Calmar ratio: annualised mean return / max drawdown
        // Compute raw mean directly from sorted array (avoids reversing annualisation).
        float total_return = 0.0f;
        for (int i = 0; i < sort_len; i++) {
            total_return += s_sorted[i];
        }
        float daily_mean = total_return / (float)sort_len;
        float max_dd = metrics_out[out_base + 2];
        float calmar = (max_dd > 1e-8f)
            ? (daily_mean * annualization_factor * annualization_factor) / max_dd
            : 0.0f;

        // Omega ratio: sum of positive returns / sum of |negative returns|
        float gain_sum = 0.0f, loss_sum = 0.0f;
        for (int i = 0; i < sort_len; i++) {
            if (s_sorted[i] > 0.0f) gain_sum += s_sorted[i];
            else                     loss_sum -= s_sorted[i];
        }
        float omega = (loss_sum > 1e-10f) ? gain_sum / loss_sum : 0.0f;

        metrics_out[out_base + 6] = var_95;
        metrics_out[out_base + 7] = cvar_95;
        metrics_out[out_base + 8] = calmar;
        metrics_out[out_base + 9] = omega;
    }
}
