rustsim-crowd 0.0.1

//! CUDA implementation of [Generalized Centrifugal Force](crate::generalized_centrifugal_force).
//!
//! Closes the GCF half of the P0-1 "CUDA kernels for CFS, AVM, GCF,
//! OSM" item from `docs/rustsim-crowd.md`. Mirrors the
//! [`crate::cuda::social_force`] design exactly:
//!
//! - One CUDA thread per pedestrian. O(n²) pair loop inside the
//!   kernel; same double-buffered `*_in` / `*_out` column layout as
//!   the SFM kernel.
//! - Internal `f32` arithmetic (GCF is force-based and similarly
//!   tolerant to single precision; the CPU reference runs in `f64`
//!   and the test harness allows `1e-2 m` divergence).
//! - NVRTC at runtime, no CUDA toolkit needed at build time.
//! - `CudaState::step` is the stateless per-tick entry point;
//!   [`CudaResident`] keeps device buffers alive across ticks;
//!   [`step_with_fallback`] is the production wrapper that retries on
//!   CPU when CUDA init or the kernel launch fails.
//!
//! See `src/cuda/social_force.rs` for the full design rationale.

use cudarc::driver::{
    CudaContext, CudaFunction, CudaSlice, CudaStream, LaunchConfig, PushKernelArg,
};
use cudarc::nvrtc::compile_ptx;
use std::sync::Arc;

use crate::common::{Pedestrian, WallSegment};
use crate::generalized_centrifugal_force::Params;

/// CUDA source for the GCF step kernel.
///
/// Force law mirrors `generalized_centrifugal_force::ped_force` /
/// `wall_force` exactly, modulo `f32` precision:
/// `F_ij = m * (v0 + approach)^2 / clearance * e_ji`,
/// `clearance = max(d - r_i - r_j, min_clearance)`,
/// `r(v) = a + b * |v|`,
/// `approach = max(0, -dot(v_i - v_j, e_ji))`,
/// `F_wall = wall_gain * m * v0^2 / clearance_wall * e_wall` against
/// the closest point of every wall segment.
///
/// Integrator is semi-implicit Euler (`v ← clamp(v + a·dt, max_speed); p ← p + v·dt`)
/// with `a` clamped to `max_accel` first — exactly as on the CPU.
const GCF_CUDA_SRC: &str = r#"
extern "C" __global__ void gcf_step(
    const float* __restrict__ pos_x_in,
    const float* __restrict__ pos_y_in,
    const float* __restrict__ vel_x_in,
    const float* __restrict__ vel_y_in,
    const float* __restrict__ desired_speed,
    const float* __restrict__ dest_x,
    const float* __restrict__ dest_y,
    float* __restrict__ pos_x_out,
    float* __restrict__ pos_y_out,
    float* __restrict__ vel_x_out,
    float* __restrict__ vel_y_out,
    const float* __restrict__ wall_ax,
    const float* __restrict__ wall_ay,
    const float* __restrict__ wall_bx,
    const float* __restrict__ wall_by,
    unsigned int n,
    unsigned int n_walls,
    float tau,
    float a_radius,
    float b_radius,
    float wall_gain,
    float mass,
    float max_speed,
    float max_accel,
    float min_clearance,
    float arrival_radius,
    float dt)
{
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= n) return;

    float pxi = pos_x_in[i];
    float pyi = pos_y_in[i];
    float vxi = vel_x_in[i];
    float vyi = vel_y_in[i];
    float v0  = desired_speed[i];
    float dxi = dest_x[i];
    float dyi = dest_y[i];

    // Body radius depends on speed: r(v) = a + b * |v|.
    float speed_i = sqrtf(vxi * vxi + vyi * vyi);
    float ri = a_radius + b_radius * speed_i;

    // Driving force with arrival taper (mirrors CPU `driving_force`).
    float to_dest_x = dxi - pxi;
    float to_dest_y = dyi - pyi;
    float d_to_dest = sqrtf(to_dest_x * to_dest_x + to_dest_y * to_dest_y);
    float eff_v0 = v0;
    if (arrival_radius > 0.0f && d_to_dest < arrival_radius) {
        eff_v0 = v0 * (d_to_dest / arrival_radius);
    }
    float dir_x = 0.0f;
    float dir_y = 0.0f;
    if (d_to_dest > 1.0e-12f) {
        dir_x = to_dest_x / d_to_dest;
        dir_y = to_dest_y / d_to_dest;
    }
    float fx = mass * (eff_v0 * dir_x - vxi) / tau;
    float fy = mass * (eff_v0 * dir_y - vyi) / tau;

    // Pedestrian centrifugal repulsion — O(n²) pair loop.
    for (unsigned int j = 0; j < n; ++j) {
        if (j == i) continue;
        float ddx = pxi - pos_x_in[j];
        float ddy = pyi - pos_y_in[j];
        float dist = sqrtf(ddx * ddx + ddy * ddy);
        if (dist < 1.0e-9f) continue;
        float inv_d = 1.0f / dist;
        float ex = ddx * inv_d;
        float ey = ddy * inv_d;

        // Approach speed along e_ji = diff/d (positive part).
        float vrx = vxi - vel_x_in[j];
        float vry = vyi - vel_y_in[j];
        float approach = -(vrx * ex + vry * ey);
        if (approach < 0.0f) approach = 0.0f;

        float speed_j = sqrtf(
            vel_x_in[j] * vel_x_in[j] + vel_y_in[j] * vel_y_in[j]);
        float rj = a_radius + b_radius * speed_j;
        float clearance = dist - ri - rj;
        if (clearance < min_clearance) clearance = min_clearance;

        float vsum = v0 + approach;
        float mag = mass * vsum * vsum / clearance;
        fx += mag * ex;
        fy += mag * ey;
    }

    // Wall centrifugal repulsion against the closest segment point.
    for (unsigned int k = 0; k < n_walls; ++k) {
        float ax = wall_ax[k];
        float ay = wall_ay[k];
        float bx = wall_bx[k];
        float by = wall_by[k];
        float abx = bx - ax;
        float aby = by - ay;
        float denom = abx * abx + aby * aby;
        float t = 0.0f;
        if (denom > 1.0e-18f) {
            t = ((pxi - ax) * abx + (pyi - ay) * aby) / denom;
            if (t < 0.0f) t = 0.0f;
            if (t > 1.0f) t = 1.0f;
        }
        float cpx = ax + t * abx;
        float cpy = ay + t * aby;
        float ddx = pxi - cpx;
        float ddy = pyi - cpy;
        float dist = sqrtf(ddx * ddx + ddy * ddy);
        if (dist < 1.0e-9f) continue;
        float inv_d = 1.0f / dist;
        float ex = ddx * inv_d;
        float ey = ddy * inv_d;
        float clearance = dist - ri;
        if (clearance < min_clearance) clearance = min_clearance;
        float mag = wall_gain * mass * v0 * v0 / clearance;
        fx += mag * ex;
        fy += mag * ey;
    }

    // a = F / m, clamp to max_accel for explicit-Euler stability.
    float ax_acc = fx / mass;
    float ay_acc = fy / mass;
    float a_mag = sqrtf(ax_acc * ax_acc + ay_acc * ay_acc);
    if (a_mag > max_accel && a_mag > 1.0e-12f) {
        float s = max_accel / a_mag;
        ax_acc *= s;
        ay_acc *= s;
    }

    // Semi-implicit Euler: v ← clamp(v + a·dt); p ← p + v·dt.
    float vx_new = vxi + ax_acc * dt;
    float vy_new = vyi + ay_acc * dt;
    float v_mag = sqrtf(vx_new * vx_new + vy_new * vy_new);
    if (v_mag > max_speed && v_mag > 1.0e-12f) {
        float s = max_speed / v_mag;
        vx_new *= s;
        vy_new *= s;
    }
    float px_new = pxi + vx_new * dt;
    float py_new = pyi + vy_new * dt;

    pos_x_out[i] = px_new;
    pos_y_out[i] = py_new;
    vel_x_out[i] = vx_new;
    vel_y_out[i] = vy_new;
}
"#;

/// Persistent CUDA state for the GCF kernel.
///
/// Construct once per simulation; dropping it releases every CUDA
/// resource. Mirrors [`crate::cuda::social_force::CudaState`].
pub struct CudaState {
    ctx: Arc<CudaContext>,
    stream: Arc<CudaStream>,
    func: CudaFunction,
    block_size: u32,
}

impl CudaState {
    /// Initialise CUDA on device 0, compile the GCF kernel via NVRTC,
    /// and load it. Returns `Err` if no device is available, the
    /// driver fails, or NVRTC rejects the source.
    pub fn new() -> Result<Self, String> {
        Self::with_block_size(256)
    }

    /// Like [`CudaState::new`] but with a configurable CUDA block size.
    pub fn with_block_size(block_size: u32) -> Result<Self, String> {
        if block_size == 0 {
            return Err("block_size must be positive".to_string());
        }
        let ctx = super::new_context(0)?;
        let stream = ctx.default_stream();
        let ptx = compile_ptx(GCF_CUDA_SRC).map_err(|e| format!("NVRTC compile failed: {e}"))?;
        let module = ctx
            .load_module(ptx)
            .map_err(|e| format!("module load failed: {e}"))?;
        let func = module
            .load_function("gcf_step")
            .map_err(|e| format!("kernel lookup failed: {e}"))?;
        Ok(Self {
            ctx,
            stream,
            func,
            block_size,
        })
    }

    /// Advance `peds` by one tick of `dt` seconds on the GPU.
    ///
    /// Stateless w.r.t. device memory: every column is uploaded and
    /// downloaded per call. Returns the kernel execution time in
    /// microseconds.
    pub fn step(
        &self,
        peds: &mut [Pedestrian],
        walls: &[WallSegment],
        params: &Params,
        dt: f64,
    ) -> Result<u128, String> {
        let n = peds.len();
        if n == 0 {
            return Ok(0);
        }

        let stream = &self.stream;
        let mut pos_x: Vec<f32> = Vec::with_capacity(n);
        let mut pos_y: Vec<f32> = Vec::with_capacity(n);
        let mut vel_x: Vec<f32> = Vec::with_capacity(n);
        let mut vel_y: Vec<f32> = Vec::with_capacity(n);
        let mut desired_speed: Vec<f32> = Vec::with_capacity(n);
        let mut dest_x: Vec<f32> = Vec::with_capacity(n);
        let mut dest_y: Vec<f32> = Vec::with_capacity(n);
        for p in peds.iter() {
            pos_x.push(p.pos[0] as f32);
            pos_y.push(p.pos[1] as f32);
            vel_x.push(p.vel[0] as f32);
            vel_y.push(p.vel[1] as f32);
            desired_speed.push(p.desired_speed as f32);
            dest_x.push(p.destination[0] as f32);
            dest_y.push(p.destination[1] as f32);
        }

        let n_walls = walls.len();
        let mut wall_ax: Vec<f32> = Vec::with_capacity(n_walls.max(1));
        let mut wall_ay: Vec<f32> = Vec::with_capacity(n_walls.max(1));
        let mut wall_bx: Vec<f32> = Vec::with_capacity(n_walls.max(1));
        let mut wall_by: Vec<f32> = Vec::with_capacity(n_walls.max(1));
        if n_walls == 0 {
            // CUDA disallows zero-sized device allocations; pad with a
            // single degenerate segment that the kernel skips via
            // `n_walls == 0`.
            wall_ax.push(0.0);
            wall_ay.push(0.0);
            wall_bx.push(0.0);
            wall_by.push(0.0);
        } else {
            for w in walls {
                wall_ax.push(w.a[0] as f32);
                wall_ay.push(w.a[1] as f32);
                wall_bx.push(w.b[0] as f32);
                wall_by.push(w.b[1] as f32);
            }
        }

        let d_pos_x: CudaSlice<f32> = stream
            .clone_htod(&pos_x)
            .map_err(|e| format!("htod pos_x failed: {e}"))?;
        let d_pos_y: CudaSlice<f32> = stream
            .clone_htod(&pos_y)
            .map_err(|e| format!("htod pos_y failed: {e}"))?;
        let d_vel_x: CudaSlice<f32> = stream
            .clone_htod(&vel_x)
            .map_err(|e| format!("htod vel_x failed: {e}"))?;
        let d_vel_y: CudaSlice<f32> = stream
            .clone_htod(&vel_y)
            .map_err(|e| format!("htod vel_y failed: {e}"))?;
        let d_desired_speed: CudaSlice<f32> = stream
            .clone_htod(&desired_speed)
            .map_err(|e| format!("htod desired_speed failed: {e}"))?;
        let d_dest_x: CudaSlice<f32> = stream
            .clone_htod(&dest_x)
            .map_err(|e| format!("htod dest_x failed: {e}"))?;
        let d_dest_y: CudaSlice<f32> = stream
            .clone_htod(&dest_y)
            .map_err(|e| format!("htod dest_y failed: {e}"))?;
        let mut d_pos_x_out: CudaSlice<f32> = stream
            .alloc_zeros(n)
            .map_err(|e| format!("alloc pos_x_out failed: {e}"))?;
        let mut d_pos_y_out: CudaSlice<f32> = stream
            .alloc_zeros(n)
            .map_err(|e| format!("alloc pos_y_out failed: {e}"))?;
        let mut d_vel_x_out: CudaSlice<f32> = stream
            .alloc_zeros(n)
            .map_err(|e| format!("alloc vel_x_out failed: {e}"))?;
        let mut d_vel_y_out: CudaSlice<f32> = stream
            .alloc_zeros(n)
            .map_err(|e| format!("alloc vel_y_out failed: {e}"))?;
        let d_wall_ax: CudaSlice<f32> = stream
            .clone_htod(&wall_ax)
            .map_err(|e| format!("htod wall_ax failed: {e}"))?;
        let d_wall_ay: CudaSlice<f32> = stream
            .clone_htod(&wall_ay)
            .map_err(|e| format!("htod wall_ay failed: {e}"))?;
        let d_wall_bx: CudaSlice<f32> = stream
            .clone_htod(&wall_bx)
            .map_err(|e| format!("htod wall_bx failed: {e}"))?;
        let d_wall_by: CudaSlice<f32> = stream
            .clone_htod(&wall_by)
            .map_err(|e| format!("htod wall_by failed: {e}"))?;

        let n_u32 = n as u32;
        let n_walls_u32 = n_walls as u32;
        let tau = params.tau as f32;
        let a_radius = params.a as f32;
        let b_radius = params.b as f32;
        let wall_gain = params.wall_gain as f32;
        let mass = params.mass as f32;
        let max_speed = params.max_speed as f32;
        let max_accel = params.max_accel as f32;
        let min_clearance = params.min_clearance as f32;
        let arrival_radius = params.arrival_radius as f32;
        let dt_f32 = dt as f32;

        let grid = n.div_ceil(self.block_size as usize) as u32;
        let cfg = LaunchConfig {
            grid_dim: (grid.max(1), 1, 1),
            block_dim: (self.block_size, 1, 1),
            shared_mem_bytes: 0,
        };

        let t0 = std::time::Instant::now();

        // SAFETY:
        // - every device slice was allocated on `stream` and outlives
        //   the launch;
        // - the kernel bounds-checks every access against `n` and
        //   `n_walls`;
        // - argument order matches `gcf_step` exactly;
        // - `stream.synchronize()` runs before any host-visible read.
        unsafe {
            let mut b = stream.launch_builder(&self.func);
            b.arg(&d_pos_x);
            b.arg(&d_pos_y);
            b.arg(&d_vel_x);
            b.arg(&d_vel_y);
            b.arg(&d_desired_speed);
            b.arg(&d_dest_x);
            b.arg(&d_dest_y);
            b.arg(&mut d_pos_x_out);
            b.arg(&mut d_pos_y_out);
            b.arg(&mut d_vel_x_out);
            b.arg(&mut d_vel_y_out);
            b.arg(&d_wall_ax);
            b.arg(&d_wall_ay);
            b.arg(&d_wall_bx);
            b.arg(&d_wall_by);
            b.arg(&n_u32);
            b.arg(&n_walls_u32);
            b.arg(&tau);
            b.arg(&a_radius);
            b.arg(&b_radius);
            b.arg(&wall_gain);
            b.arg(&mass);
            b.arg(&max_speed);
            b.arg(&max_accel);
            b.arg(&min_clearance);
            b.arg(&arrival_radius);
            b.arg(&dt_f32);
            b.launch(cfg)
                .map_err(|e| format!("kernel launch failed: {e}"))?;
        }

        stream
            .synchronize()
            .map_err(|e| format!("stream sync failed: {e}"))?;
        let kernel_us = t0.elapsed().as_micros();

        stream
            .memcpy_dtoh(&d_pos_x_out, &mut pos_x)
            .map_err(|e| format!("dtoh pos_x failed: {e}"))?;
        stream
            .memcpy_dtoh(&d_pos_y_out, &mut pos_y)
            .map_err(|e| format!("dtoh pos_y failed: {e}"))?;
        stream
            .memcpy_dtoh(&d_vel_x_out, &mut vel_x)
            .map_err(|e| format!("dtoh vel_x failed: {e}"))?;
        stream
            .memcpy_dtoh(&d_vel_y_out, &mut vel_y)
            .map_err(|e| format!("dtoh vel_y failed: {e}"))?;

        for (i, p) in peds.iter_mut().enumerate() {
            p.pos = [pos_x[i] as f64, pos_y[i] as f64];
            p.vel = [vel_x[i] as f64, vel_y[i] as f64];
        }

        let _ = &self.ctx; // anchor ctx lifetime
        Ok(kernel_us)
    }
}

/// CUDA state with persistent device-resident buffers for GCF.
///
/// Uploads the initial pedestrian/wall columns once, then advances the
/// simulation with zero host-device traffic until [`download`](Self::download)
/// is called. Population, desired speeds, destinations, and wall segments are
/// fixed for the lifetime of the resident store.
pub struct CudaResident {
    ctx: Arc<CudaContext>,
    stream: Arc<CudaStream>,
    func: CudaFunction,
    block_size: u32,

    n: usize,
    n_walls: usize,

    d_pos_x_a: CudaSlice<f32>,
    d_pos_y_a: CudaSlice<f32>,
    d_vel_x_a: CudaSlice<f32>,
    d_vel_y_a: CudaSlice<f32>,
    d_pos_x_b: CudaSlice<f32>,
    d_pos_y_b: CudaSlice<f32>,
    d_vel_x_b: CudaSlice<f32>,
    d_vel_y_b: CudaSlice<f32>,

    d_radius: CudaSlice<f32>,
    d_desired_speed: CudaSlice<f32>,
    d_dest_x: CudaSlice<f32>,
    d_dest_y: CudaSlice<f32>,

    d_wall_ax: CudaSlice<f32>,
    d_wall_ay: CudaSlice<f32>,
    d_wall_bx: CudaSlice<f32>,
    d_wall_by: CudaSlice<f32>,
}

impl CudaResident {
    /// Upload an initial pedestrian and wall configuration to the device.
    pub fn upload(peds: &[Pedestrian], walls: &[WallSegment]) -> Result<Self, String> {
        Self::upload_with_block_size(peds, walls, 256)
    }

    /// Like [`CudaResident::upload`] but with a configurable CUDA block size.
    pub fn upload_with_block_size(
        peds: &[Pedestrian],
        walls: &[WallSegment],
        block_size: u32,
    ) -> Result<Self, String> {
        if block_size == 0 {
            return Err("block_size must be positive".to_string());
        }
        let n = peds.len();
        if n == 0 {
            return Err("CudaResident::upload requires at least one pedestrian".to_string());
        }

        let ctx = super::new_context(0)?;
        let stream = ctx.default_stream();
        let ptx =
            compile_ptx(GCF_CUDA_SRC).map_err(|err| format!("NVRTC compile failed: {err}"))?;
        let module = ctx
            .load_module(ptx)
            .map_err(|err| format!("module load failed: {err}"))?;
        let func = module
            .load_function("gcf_step")
            .map_err(|err| format!("kernel lookup failed: {err}"))?;

        let mut pos_x = Vec::with_capacity(n);
        let mut pos_y = Vec::with_capacity(n);
        let mut vel_x = Vec::with_capacity(n);
        let mut vel_y = Vec::with_capacity(n);
        let mut radius = Vec::with_capacity(n);
        let mut desired_speed = Vec::with_capacity(n);
        let mut dest_x = Vec::with_capacity(n);
        let mut dest_y = Vec::with_capacity(n);
        for pedestrian in peds {
            pos_x.push(pedestrian.pos[0] as f32);
            pos_y.push(pedestrian.pos[1] as f32);
            vel_x.push(pedestrian.vel[0] as f32);
            vel_y.push(pedestrian.vel[1] as f32);
            radius.push(pedestrian.radius as f32);
            desired_speed.push(pedestrian.desired_speed as f32);
            dest_x.push(pedestrian.destination[0] as f32);
            dest_y.push(pedestrian.destination[1] as f32);
        }

        let n_walls = walls.len();
        let (wall_ax, wall_ay, wall_bx, wall_by) = if n_walls == 0 {
            (vec![0.0f32], vec![0.0f32], vec![0.0f32], vec![0.0f32])
        } else {
            let mut wall_ax = Vec::with_capacity(n_walls);
            let mut wall_ay = Vec::with_capacity(n_walls);
            let mut wall_bx = Vec::with_capacity(n_walls);
            let mut wall_by = Vec::with_capacity(n_walls);
            for wall in walls {
                wall_ax.push(wall.a[0] as f32);
                wall_ay.push(wall.a[1] as f32);
                wall_bx.push(wall.b[0] as f32);
                wall_by.push(wall.b[1] as f32);
            }
            (wall_ax, wall_ay, wall_bx, wall_by)
        };

        let d_pos_x_a = stream
            .clone_htod(&pos_x)
            .map_err(|err| format!("htod pos_x failed: {err}"))?;
        let d_pos_y_a = stream
            .clone_htod(&pos_y)
            .map_err(|err| format!("htod pos_y failed: {err}"))?;
        let d_vel_x_a = stream
            .clone_htod(&vel_x)
            .map_err(|err| format!("htod vel_x failed: {err}"))?;
        let d_vel_y_a = stream
            .clone_htod(&vel_y)
            .map_err(|err| format!("htod vel_y failed: {err}"))?;
        let d_pos_x_b = stream
            .alloc_zeros(n)
            .map_err(|err| format!("alloc pos_x_b failed: {err}"))?;
        let d_pos_y_b = stream
            .alloc_zeros(n)
            .map_err(|err| format!("alloc pos_y_b failed: {err}"))?;
        let d_vel_x_b = stream
            .alloc_zeros(n)
            .map_err(|err| format!("alloc vel_x_b failed: {err}"))?;
        let d_vel_y_b = stream
            .alloc_zeros(n)
            .map_err(|err| format!("alloc vel_y_b failed: {err}"))?;
        let d_radius = stream
            .clone_htod(&radius)
            .map_err(|err| format!("htod radius failed: {err}"))?;
        let d_desired_speed = stream
            .clone_htod(&desired_speed)
            .map_err(|err| format!("htod desired_speed failed: {err}"))?;
        let d_dest_x = stream
            .clone_htod(&dest_x)
            .map_err(|err| format!("htod dest_x failed: {err}"))?;
        let d_dest_y = stream
            .clone_htod(&dest_y)
            .map_err(|err| format!("htod dest_y failed: {err}"))?;
        let d_wall_ax = stream
            .clone_htod(&wall_ax)
            .map_err(|err| format!("htod wall_ax failed: {err}"))?;
        let d_wall_ay = stream
            .clone_htod(&wall_ay)
            .map_err(|err| format!("htod wall_ay failed: {err}"))?;
        let d_wall_bx = stream
            .clone_htod(&wall_bx)
            .map_err(|err| format!("htod wall_bx failed: {err}"))?;
        let d_wall_by = stream
            .clone_htod(&wall_by)
            .map_err(|err| format!("htod wall_by failed: {err}"))?;

        stream
            .synchronize()
            .map_err(|err| format!("initial sync failed: {err}"))?;

        Ok(Self {
            ctx,
            stream,
            func,
            block_size,
            n,
            n_walls,
            d_pos_x_a,
            d_pos_y_a,
            d_vel_x_a,
            d_vel_y_a,
            d_pos_x_b,
            d_pos_y_b,
            d_vel_x_b,
            d_vel_y_b,
            d_radius,
            d_desired_speed,
            d_dest_x,
            d_dest_y,
            d_wall_ax,
            d_wall_ay,
            d_wall_bx,
            d_wall_by,
        })
    }

    /// Number of resident pedestrians.
    #[inline]
    pub fn len(&self) -> usize {
        self.n
    }

    /// Whether any pedestrians are resident.
    #[inline]
    pub fn is_empty(&self) -> bool {
        self.n == 0
    }

    /// Advance one tick of `dt` seconds fully on the GPU.
    pub fn step(&mut self, params: &Params, dt: f64) -> Result<u128, String> {
        let n_u32 = self.n as u32;
        let n_walls_u32 = self.n_walls as u32;
        let tau = params.tau as f32;
        let a_radius = params.a as f32;
        let b_radius = params.b as f32;
        let wall_gain = params.wall_gain as f32;
        let mass = params.mass as f32;
        let max_speed = params.max_speed as f32;
        let max_accel = params.max_accel as f32;
        let min_clearance = params.min_clearance as f32;
        let arrival_radius = params.arrival_radius as f32;
        let dt_f32 = dt as f32;

        let grid = self.n.div_ceil(self.block_size as usize) as u32;
        let cfg = LaunchConfig {
            grid_dim: (grid.max(1), 1, 1),
            block_dim: (self.block_size, 1, 1),
            shared_mem_bytes: 0,
        };

        let Self {
            ref stream,
            ref func,
            ref d_desired_speed,
            ref d_dest_x,
            ref d_dest_y,
            ref d_wall_ax,
            ref d_wall_ay,
            ref d_wall_bx,
            ref d_wall_by,
            ref d_pos_x_a,
            ref d_pos_y_a,
            ref d_vel_x_a,
            ref d_vel_y_a,
            ref mut d_pos_x_b,
            ref mut d_pos_y_b,
            ref mut d_vel_x_b,
            ref mut d_vel_y_b,
            ..
        } = *self;

        let t0 = std::time::Instant::now();

        // SAFETY:
        // - every device slice is owned by `self` and outlives the launch;
        // - the kernel bounds-checks accesses against `n` and `n_walls`;
        // - argument order matches `gcf_step` exactly;
        // - `_a` buffers are read-only inputs and `_b` buffers are distinct
        //   output allocations.
        unsafe {
            let mut builder = stream.launch_builder(func);
            builder.arg(d_pos_x_a);
            builder.arg(d_pos_y_a);
            builder.arg(d_vel_x_a);
            builder.arg(d_vel_y_a);
            builder.arg(d_desired_speed);
            builder.arg(d_dest_x);
            builder.arg(d_dest_y);
            builder.arg(d_pos_x_b);
            builder.arg(d_pos_y_b);
            builder.arg(d_vel_x_b);
            builder.arg(d_vel_y_b);
            builder.arg(d_wall_ax);
            builder.arg(d_wall_ay);
            builder.arg(d_wall_bx);
            builder.arg(d_wall_by);
            builder.arg(&n_u32);
            builder.arg(&n_walls_u32);
            builder.arg(&tau);
            builder.arg(&a_radius);
            builder.arg(&b_radius);
            builder.arg(&wall_gain);
            builder.arg(&mass);
            builder.arg(&max_speed);
            builder.arg(&max_accel);
            builder.arg(&min_clearance);
            builder.arg(&arrival_radius);
            builder.arg(&dt_f32);
            builder
                .launch(cfg)
                .map_err(|err| format!("kernel launch failed: {err}"))?;
        }

        stream
            .synchronize()
            .map_err(|err| format!("stream sync failed: {err}"))?;
        let kernel_us = t0.elapsed().as_micros();

        std::mem::swap(&mut self.d_pos_x_a, &mut self.d_pos_x_b);
        std::mem::swap(&mut self.d_pos_y_a, &mut self.d_pos_y_b);
        std::mem::swap(&mut self.d_vel_x_a, &mut self.d_vel_x_b);
        std::mem::swap(&mut self.d_vel_y_a, &mut self.d_vel_y_b);

        let _ = &self.ctx;
        Ok(kernel_us)
    }

    /// Bulk download of the current resident state back into `peds`.
    pub fn download(&self, peds: &mut Vec<Pedestrian>) -> Result<(), String> {
        let n = self.n;
        peds.resize(
            n,
            Pedestrian {
                pos: [0.0, 0.0],
                vel: [0.0, 0.0],
                radius: 0.0,
                desired_speed: 0.0,
                destination: [0.0, 0.0],
            },
        );

        let mut pos_x = vec![0.0f32; n];
        let mut pos_y = vec![0.0f32; n];
        let mut vel_x = vec![0.0f32; n];
        let mut vel_y = vec![0.0f32; n];
        let mut radius = vec![0.0f32; n];
        let mut desired_speed = vec![0.0f32; n];
        let mut dest_x = vec![0.0f32; n];
        let mut dest_y = vec![0.0f32; n];

        self.stream
            .memcpy_dtoh(&self.d_pos_x_a, &mut pos_x)
            .map_err(|err| format!("dtoh pos_x failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_pos_y_a, &mut pos_y)
            .map_err(|err| format!("dtoh pos_y failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_vel_x_a, &mut vel_x)
            .map_err(|err| format!("dtoh vel_x failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_vel_y_a, &mut vel_y)
            .map_err(|err| format!("dtoh vel_y failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_radius, &mut radius)
            .map_err(|err| format!("dtoh radius failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_desired_speed, &mut desired_speed)
            .map_err(|err| format!("dtoh desired_speed failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_dest_x, &mut dest_x)
            .map_err(|err| format!("dtoh dest_x failed: {err}"))?;
        self.stream
            .memcpy_dtoh(&self.d_dest_y, &mut dest_y)
            .map_err(|err| format!("dtoh dest_y failed: {err}"))?;

        for (index, pedestrian) in peds.iter_mut().enumerate().take(n) {
            pedestrian.pos = [pos_x[index] as f64, pos_y[index] as f64];
            pedestrian.vel = [vel_x[index] as f64, vel_y[index] as f64];
            pedestrian.radius = radius[index] as f64;
            pedestrian.desired_speed = desired_speed[index] as f64;
            pedestrian.destination = [dest_x[index] as f64, dest_y[index] as f64];
        }

        Ok(())
    }
}

/// One-shot wrapper: initialise CUDA, run one GCF step, tear down.
///
/// For any production tick loop, hold a [`CudaState`] across ticks
/// instead — otherwise NVRTC compilation dominates the runtime.
pub fn step(
    peds: &mut [Pedestrian],
    walls: &[WallSegment],
    params: &Params,
    dt: f64,
) -> Result<u128, String> {
    let state = CudaState::new()?;
    state.step(peds, walls, params, dt)
}

/// CPU-fallback wrapper: try CUDA, fall back to
/// [`crate::generalized_centrifugal_force::step`] on any error.
///
/// Returns `true` iff the CUDA path succeeded.
pub fn step_with_fallback(
    state: &mut Option<CudaState>,
    peds: &mut [Pedestrian],
    walls: &[WallSegment],
    params: &Params,
    dt: f64,
) -> bool {
    if state.is_none() {
        match CudaState::new() {
            Ok(s) => *state = Some(s),
            Err(e) => {
                eprintln!("rustsim-crowd CUDA init failed ({e}), using CPU GCF path");
                #[allow(deprecated)]
                crate::generalized_centrifugal_force::step(peds, walls, params, dt);
                return false;
            }
        }
    }
    match state.as_ref().unwrap().step(peds, walls, params, dt) {
        Ok(_) => true,
        Err(e) => {
            eprintln!("rustsim-crowd CUDA GCF step failed ({e}), falling back to CPU");
            #[allow(deprecated)]
            crate::generalized_centrifugal_force::step(peds, walls, params, dt);
            false
        }
    }
}