rustsim-crowd 0.0.1

//! CUDA implementation of [Collision-Free Speed](crate::collision_free_speed).
//!
//! Closes the CFS slot of the P0-1 "CUDA kernels for CFS, AVM, OSM"
//! item from `docs/rustsim-crowd.md`. Mirrors the SFM/GCF design:
//!
//! - One CUDA thread per pedestrian. Two O(n²) pair loops inside the
//!   kernel (one for the chosen-direction blend, one for the free
//!   headroom along that direction); same double-buffered `*_in` /
//!   `*_out` column layout.
//! - Internal `f32` arithmetic (CFS is first-order kinematic, the
//!   speed cap dominates over numerical noise; the test harness
//!   tolerates `5e-2 m` divergence after 5 ticks).
//! - 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.
//!
//! Unlike SFM/GCF (force-based, second-order semi-implicit Euler),
//! CFS is **kinematic**: each agent's new velocity is a direction
//! blend (alpha-weighted desired + neighbour/wall repulsions) scaled
//! by a speed clamped to `(s - radius) / time_gap`, where `s` is the
//! free headroom along the chosen direction. The kernel writes both
//! `vel_out` and `pos_out` directly without going through an
//! acceleration column.
//!
//! 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::collision_free_speed::Params;
use crate::common::{Pedestrian, WallSegment};

/// CUDA source for the CFS step kernel.
///
/// Mirrors `collision_free_speed::chosen_direction` +
/// `free_headroom` exactly (modulo `f32` precision):
///
/// 1. `acc = alpha * e_dest + Σ_j w_ij * e_ij + Σ_walls w_iw * e_iw`,
///    with `w = strength * exp(-max(0, clearance) / range)`.
/// 2. `e = normalize(acc)`; if `acc ≈ 0`, `e = 0`.
/// 3. `s = min over forward neighbours of |q - p|`, with the same
///    "in front and laterally within `r_i + r_j`" filter as the CPU.
/// 4. `v = min(eff_v0, max(0, (s - radius_i) / time_gap))`, where
///    `eff_v0` applies the arrival-radius linear taper.
/// 5. `new_vel = e * v`; `new_pos = p + new_vel * dt`.
const CFS_CUDA_SRC: &str = r#"
extern "C" __global__ void cfs_step(
    const float* __restrict__ pos_x_in,
    const float* __restrict__ pos_y_in,
    const float* __restrict__ radius,
    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 time_gap,
    float alpha,
    float interaction_range,
    float interaction_strength,
    float wall_range,
    float wall_strength,
    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 ri  = radius[i];
    float v0  = desired_speed[i];
    float dxi = dest_x[i];
    float dyi = dest_y[i];

    // Desired direction toward destination + arrival taper on v0.
    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 edx = 0.0f;
    float edy = 0.0f;
    if (d_to_dest > 1.0e-12f) {
        edx = to_dest_x / d_to_dest;
        edy = to_dest_y / d_to_dest;
    }
    float eff_v0 = v0;
    if (arrival_radius > 0.0f && d_to_dest < arrival_radius) {
        eff_v0 = v0 * (d_to_dest / arrival_radius);
    }

    // Direction-blend accumulator: alpha * e_dest + Σ repulsions.
    float ax_dir = alpha * edx;
    float ay_dir = alpha * edy;

    // Pedestrian repulsion contribution to direction.
    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 d = sqrtf(ddx * ddx + ddy * ddy);
        if (d < 1.0e-9f) continue;
        float inv_d = 1.0f / d;
        float ex = ddx * inv_d;
        float ey = ddy * inv_d;
        float clearance = d - (ri + radius[j]);
        if (clearance < 0.0f) clearance = 0.0f;
        float weight =
            interaction_strength * expf(-clearance / interaction_range);
        ax_dir += weight * ex;
        ay_dir += weight * ey;
    }

    // Wall repulsion contribution to direction.
    for (unsigned int k = 0; k < n_walls; ++k) {
        float wax = wall_ax[k];
        float way = wall_ay[k];
        float wbx = wall_bx[k];
        float wby = wall_by[k];
        float abx = wbx - wax;
        float aby = wby - way;
        float denom = abx * abx + aby * aby;
        float t = 0.0f;
        if (denom > 1.0e-18f) {
            t = ((pxi - wax) * abx + (pyi - way) * aby) / denom;
            if (t < 0.0f) t = 0.0f;
            if (t > 1.0f) t = 1.0f;
        }
        float cpx = wax + t * abx;
        float cpy = way + t * aby;
        float ddx = pxi - cpx;
        float ddy = pyi - cpy;
        float d = sqrtf(ddx * ddx + ddy * ddy);
        if (d < 1.0e-9f) continue;
        float inv_d = 1.0f / d;
        float ex = ddx * inv_d;
        float ey = ddy * inv_d;
        float clearance = d - ri;
        if (clearance < 0.0f) clearance = 0.0f;
        float weight = wall_strength * expf(-clearance / wall_range);
        ax_dir += weight * ex;
        ay_dir += weight * ey;
    }

    // Normalize the chosen direction; zero accumulator → zero direction.
    float a_mag = sqrtf(ax_dir * ax_dir + ay_dir * ay_dir);
    float ex = 0.0f;
    float ey = 0.0f;
    if (a_mag > 1.0e-12f) {
        ex = ax_dir / a_mag;
        ey = ay_dir / a_mag;
    }

    // Free headroom along (ex, ey): nearest forward neighbour whose
    // lateral offset is within ri + r_j of the line.
    float s = 1.0e30f;
    for (unsigned int j = 0; j < n; ++j) {
        if (j == i) continue;
        float relx = pos_x_in[j] - pxi;
        float rely = pos_y_in[j] - pyi;
        float forward = relx * ex + rely * ey;
        if (forward <= 0.0f) continue;
        float projx = ex * forward;
        float projy = ey * forward;
        float latx = relx - projx;
        float laty = rely - projy;
        float lat = sqrtf(latx * latx + laty * laty);
        float r_sum = ri + radius[j];
        if (lat > r_sum) continue;
        float d = sqrtf(relx * relx + rely * rely);
        if (d < s) s = d;
    }

    // Speed cap from headroom; clamp eff_v0 to it.
    float cap = (s - ri) / time_gap;
    if (cap < 0.0f) cap = 0.0f;
    float v = eff_v0 < cap ? eff_v0 : cap;

    float vx_new = ex * v;
    float vy_new = ey * v;
    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 CFS kernel.
///
/// Construct once per simulation; dropping it releases every CUDA
/// resource. Mirrors [`crate::cuda::generalized_centrifugal_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 CFS 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(CFS_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("cfs_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![0.0; n];
        let mut vel_y: Vec<f32> = vec![0.0; n];
        let mut radius: 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);
            radius.push(p.radius 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_radius: CudaSlice<f32> = stream
            .clone_htod(&radius)
            .map_err(|e| format!("htod radius 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 time_gap = params.time_gap as f32;
        let alpha = params.alpha as f32;
        let interaction_range = params.interaction_range as f32;
        let interaction_strength = params.interaction_strength as f32;
        let wall_range = params.wall_range as f32;
        let wall_strength = params.wall_strength 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 `cfs_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_radius);
            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(&time_gap);
            b.arg(&alpha);
            b.arg(&interaction_range);
            b.arg(&interaction_strength);
            b.arg(&wall_range);
            b.arg(&wall_strength);
            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 CFS.
///
/// Uploads the initial pedestrian/wall columns once, then advances the
/// simulation with zero host-device traffic until [`download`](Self::download)
/// is called. The population, radii, desired speeds, destinations, and wall
/// segments are fixed for the lifetime of the resident store; create a new
/// resident store when any of those static columns change.
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(CFS_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("cfs_step")
            .map_err(|e| format!("kernel lookup failed: {e}"))?;

        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 p in peds {
            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);
            radius.push(p.radius 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 (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 ax = Vec::with_capacity(n_walls);
            let mut ay = Vec::with_capacity(n_walls);
            let mut bx = Vec::with_capacity(n_walls);
            let mut by = Vec::with_capacity(n_walls);
            for w in walls {
                ax.push(w.a[0] as f32);
                ay.push(w.a[1] as f32);
                bx.push(w.b[0] as f32);
                by.push(w.b[1] as f32);
            }
            (ax, ay, bx, by)
        };

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

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

        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 time_gap = params.time_gap as f32;
        let alpha = params.alpha as f32;
        let interaction_range = params.interaction_range as f32;
        let interaction_strength = params.interaction_strength as f32;
        let wall_range = params.wall_range as f32;
        let wall_strength = params.wall_strength 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_radius,
            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 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 `cfs_step` exactly;
        // - `_a` position buffers are read-only inputs and `_b` buffers are
        //   distinct output allocations.
        unsafe {
            let mut b = stream.launch_builder(func);
            b.arg(d_pos_x_a);
            b.arg(d_pos_y_a);
            b.arg(d_radius);
            b.arg(d_desired_speed);
            b.arg(d_dest_x);
            b.arg(d_dest_y);
            b.arg(d_pos_x_b);
            b.arg(d_pos_y_b);
            b.arg(d_vel_x_b);
            b.arg(d_vel_y_b);
            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(&time_gap);
            b.arg(&alpha);
            b.arg(&interaction_range);
            b.arg(&interaction_strength);
            b.arg(&wall_range);
            b.arg(&wall_strength);
            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();

        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(|e| format!("dtoh pos_x failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_pos_y_a, &mut pos_y)
            .map_err(|e| format!("dtoh pos_y failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_vel_x_a, &mut vel_x)
            .map_err(|e| format!("dtoh vel_x failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_vel_y_a, &mut vel_y)
            .map_err(|e| format!("dtoh vel_y failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_radius, &mut radius)
            .map_err(|e| format!("dtoh radius failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_desired_speed, &mut desired_speed)
            .map_err(|e| format!("dtoh desired_speed failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_dest_x, &mut dest_x)
            .map_err(|e| format!("dtoh dest_x failed: {e}"))?;
        self.stream
            .memcpy_dtoh(&self.d_dest_y, &mut dest_y)
            .map_err(|e| format!("dtoh dest_y failed: {e}"))?;

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

        Ok(())
    }
}

/// One-shot wrapper: initialise CUDA, run one CFS 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::collision_free_speed::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 CFS path");
                #[allow(deprecated)]
                crate::collision_free_speed::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 CFS step failed ({e}), falling back to CPU");
            #[allow(deprecated)]
            crate::collision_free_speed::step(peds, walls, params, dt);
            false
        }
    }
}