// ═══════════════════════════════════════════════════════════════════════════
// CUDA kernels for GPU-resident Integrated Gradients (IG).
//
// Sundararajan et al. (2017). Attributions are obtained by integrating the
// model's gradient along the straight-line path baseline -> input:
//
//     attribution_i = (x_i - baseline_i) * (1/M) * sum_{k=0..M-1} d F/d x_i
//                                                  evaluated at the k-th
//                                                  interpolation point.
//
// These kernels handle the two per-step GPU primitives:
//
//   1. `interpolate_input`   : builds x(alpha) = baseline + alpha * (input - baseline)
//   2. `perturb_dimension`   : copies a buffer and perturbs a single coordinate
//                              (used to build x_plus / x_minus for central
//                              finite differences).
//
// All kernels use a 1-D grid-stride layout with 256 threads per block.
// Launch config: blocks = ceil(n / 256), block_dim = 256, shared_mem = 0.
// Deterministic — no atomicAdd, no warp-level reductions, no randomness.
// ═══════════════════════════════════════════════════════════════════════════

// Build: output[i] = baseline[i] + alpha * diff[i]    (element-wise)
//
// Launch: grid_dim = (ceil(n / 256), 1, 1), block_dim = (256, 1, 1).
extern "C" __global__
void interpolate_input(
    const float* __restrict__ baseline,
    const float* __restrict__ diff,
    float*       __restrict__ output,
    float alpha,
    int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        output[i] = baseline[i] + alpha * diff[i];
    }
}

// Copy src to dst, then apply dst[dim] += delta.
//
// The full copy is performed on every thread (rather than a targeted single
// write) so dst is always a complete, coherent snapshot of src with exactly
// one perturbed coordinate. This avoids needing a separate cudaMemcpy before
// the perturbation and keeps the operation a single kernel launch per probe.
//
// Launch: grid_dim = (ceil(n / 256), 1, 1), block_dim = (256, 1, 1).
extern "C" __global__
void perturb_dimension(
    const float* __restrict__ src,
    float*       __restrict__ dst,
    int dim,
    float delta,
    int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        float val = src[i];
        if (i == dim) {
            val += delta;
        }
        dst[i] = val;
    }
}