rustsim_spaces/

density.rs

//! Spatial density grid and Level-of-Service (LoS) analysis for transportation.
//!
//! LoS is not a single metric -- the Highway Capacity Manual (HCM) defines
//! separate criteria for different facility types and modes:
//!
//! | Criteria type | Input metric | Typical use |
//! |---------------|-------------|-------------|
//! | [`PedestrianWalkway`] | Density (pax/m2) | Sidewalks, corridors (Fruin) |
//! | [`PedestrianStairway`] | Density (pax/m2) | Stairs, escalator approaches |
//! | [`PedestrianQueuing`] | Density (pax/m2) | Waiting areas, platforms |
//! | [`VehicularFreeway`] | Density (pc/km/ln) | Freeway segments |
//! | [`VehicularUrbanStreet`] | Delay (s/veh) | Signalized intersections, arterials |
//! | [`BicycleFacility`] | Events per min | Bike lanes, shared paths |
//! | [`TransitCapacity`] | Load factor | Bus, rail passenger loading |
//! | [`CustomLosCriteria`] | User-defined thresholds | Any metric |
//!
//! # Architecture
//!
//! The [`LosCriteria`] trait defines how a raw measurement is classified
//! into a [`LosGrade`] (A through F). Built-in implementations cover
//! the standard HCM facility types. Users can implement the trait for
//! custom metrics.
//!
//! [`DensityGrid`] is a spatial grid that counts agents per cell and
//! can classify each cell using any [`LosCriteria`] implementation.
//!
//! # Usage
//!
//! ```ignore
//! use rustsim_spaces::density::*;
//!
//! let mut grid = DensityGrid::new(100.0, 50.0, 2.0);
//!
//! for agent in model.agents() {
//!     grid.add_position(agent.x, agent.y);
//! }
//!
//! // Pedestrian walkway LoS (Fruin / HCM)
//! let los = PedestrianWalkway.classify(grid.density_at(10.0, 20.0));
//!
//! // Pedestrian stairway LoS
//! let los = PedestrianStairway.classify(grid.density_at(5.0, 3.0));
//!
//! // Full statistics with a specific criteria
//! let stats = grid.statistics(&PedestrianWalkway);
//!
//! // Vehicular delay-based LoS (not grid-based -- classify directly)
//! let los = VehicularUrbanStreet.classify(35.0); // 35 s/veh delay
//!
//! grid.clear();
//! ```

// ---------------------------------------------------------------------------
// LoS grade
// ---------------------------------------------------------------------------

use thiserror::Error;

/// Level-of-Service grade A through F (HCM convention).
///
/// Ordered from best (A) to worst (F). The meaning of each grade
/// depends on the facility type and the [`LosCriteria`] used.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum LosGrade {
    /// Best operating conditions.
    A,
    /// Stable flow / minor conflicts.
    B,
    /// Stable flow with some restrictions.
    C,
    /// Approaching unstable flow.
    D,
    /// Unstable flow / significant congestion.
    E,
    /// Breakdown / forced flow.
    F,
}

impl std::fmt::Display for LosGrade {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::A => write!(f, "A"),
            Self::B => write!(f, "B"),
            Self::C => write!(f, "C"),
            Self::D => write!(f, "D"),
            Self::E => write!(f, "E"),
            Self::F => write!(f, "F"),
        }
    }
}

// ---------------------------------------------------------------------------
// LosCriteria trait
// ---------------------------------------------------------------------------

/// Trait for classifying a raw measurement into a [`LosGrade`].
///
/// Each implementation encodes the HCM thresholds for a specific
/// facility type and mode. The `value` parameter is the measurement
/// (density, delay, events/min, load factor, etc.).
pub trait LosCriteria {
    /// Classify a measurement value into a LoS grade.
    fn classify(&self, value: f64) -> LosGrade;

    /// Human-readable name of this criteria type.
    fn name(&self) -> &'static str;

    /// Unit of the input measurement (e.g. "pax/m2", "s/veh").
    fn unit(&self) -> &'static str;
}

// ---------------------------------------------------------------------------
// Pedestrian walkway (Fruin / HCM Exhibit 24-1)
// ---------------------------------------------------------------------------

/// Pedestrian walkway LoS based on density (Fruin / HCM Chapter 24).
///
/// | Grade | Density (pax/m2) | Space (m2/pax) | Description |
/// |-------|-----------------|----------------|-------------|
/// | A     | <= 0.31         | >= 3.24        | Free flow, no conflicts |
/// | B     | <= 0.43         | >= 2.32        | Minor conflicts |
/// | C     | <= 0.72         | >= 1.39        | Restricted speed |
/// | D     | <= 1.08         | >= 0.93        | Severely restricted |
/// | E     | <= 2.17         | >= 0.46        | Shuffling gait |
/// | F     | > 2.17          | < 0.46         | Breakdown |
#[derive(Debug, Clone, Copy)]
pub struct PedestrianWalkway;

impl LosCriteria for PedestrianWalkway {
    fn classify(&self, density: f64) -> LosGrade {
        if density <= 0.31 {
            LosGrade::A
        } else if density <= 0.43 {
            LosGrade::B
        } else if density <= 0.72 {
            LosGrade::C
        } else if density <= 1.08 {
            LosGrade::D
        } else if density <= 2.17 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Pedestrian Walkway"
    }
    fn unit(&self) -> &'static str {
        "pax/m2"
    }
}

// ---------------------------------------------------------------------------
// Pedestrian stairway (HCM Exhibit 24-2)
// ---------------------------------------------------------------------------

/// Pedestrian stairway LoS based on density (HCM Chapter 24).
///
/// Thresholds are tighter than walkways because stairways have
/// reduced maneuvering space and slower speeds.
///
/// | Grade | Density (pax/m2) | Space (m2/pax) |
/// |-------|-----------------|----------------|
/// | A     | <= 0.54         | >= 1.85        |
/// | B     | <= 0.72         | >= 1.39        |
/// | C     | <= 1.08         | >= 0.93        |
/// | D     | <= 1.54         | >= 0.65        |
/// | E     | <= 2.70         | >= 0.37        |
/// | F     | > 2.70          | < 0.37         |
#[derive(Debug, Clone, Copy)]
pub struct PedestrianStairway;

impl LosCriteria for PedestrianStairway {
    fn classify(&self, density: f64) -> LosGrade {
        if density <= 0.54 {
            LosGrade::A
        } else if density <= 0.72 {
            LosGrade::B
        } else if density <= 1.08 {
            LosGrade::C
        } else if density <= 1.54 {
            LosGrade::D
        } else if density <= 2.70 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Pedestrian Stairway"
    }
    fn unit(&self) -> &'static str {
        "pax/m2"
    }
}

// ---------------------------------------------------------------------------
// Pedestrian queuing (HCM Exhibit 24-3)
// ---------------------------------------------------------------------------

/// Pedestrian queuing area LoS based on density (HCM Chapter 24).
///
/// For waiting areas, platforms, ticket halls, etc.
///
/// | Grade | Density (pax/m2) | Space (m2/pax) |
/// |-------|-----------------|----------------|
/// | A     | <= 0.83         | >= 1.21        |
/// | B     | <= 1.11         | >= 0.90        |
/// | C     | <= 1.43         | >= 0.70        |
/// | D     | <= 3.33         | >= 0.30        |
/// | E     | <= 5.00         | >= 0.20        |
/// | F     | > 5.00          | < 0.20         |
#[derive(Debug, Clone, Copy)]
pub struct PedestrianQueuing;

impl LosCriteria for PedestrianQueuing {
    fn classify(&self, density: f64) -> LosGrade {
        if density <= 0.83 {
            LosGrade::A
        } else if density <= 1.11 {
            LosGrade::B
        } else if density <= 1.43 {
            LosGrade::C
        } else if density <= 3.33 {
            LosGrade::D
        } else if density <= 5.00 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Pedestrian Queuing"
    }
    fn unit(&self) -> &'static str {
        "pax/m2"
    }
}

// ---------------------------------------------------------------------------
// Vehicular freeway (HCM Chapter 12, density-based)
// ---------------------------------------------------------------------------

/// Vehicular freeway LoS based on density (HCM Chapter 12).
///
/// Input: density in passenger cars per km per lane (pc/km/ln).
///
/// | Grade | Density (pc/km/ln) | Description |
/// |-------|-------------------|-------------|
/// | A     | <= 7              | Free flow   |
/// | B     | <= 11             | Reasonably free flow |
/// | C     | <= 16             | Stable flow |
/// | D     | <= 22             | Approaching unstable |
/// | E     | <= 28             | Unstable flow |
/// | F     | > 28              | Forced / breakdown |
#[derive(Debug, Clone, Copy)]
pub struct VehicularFreeway;

impl LosCriteria for VehicularFreeway {
    fn classify(&self, density_pc_km_ln: f64) -> LosGrade {
        if density_pc_km_ln <= 7.0 {
            LosGrade::A
        } else if density_pc_km_ln <= 11.0 {
            LosGrade::B
        } else if density_pc_km_ln <= 16.0 {
            LosGrade::C
        } else if density_pc_km_ln <= 22.0 {
            LosGrade::D
        } else if density_pc_km_ln <= 28.0 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Vehicular Freeway"
    }
    fn unit(&self) -> &'static str {
        "pc/km/ln"
    }
}

// ---------------------------------------------------------------------------
// Vehicular urban street / signalized intersection (HCM Ch 16/19, delay-based)
// ---------------------------------------------------------------------------

/// Vehicular urban street and signalized intersection LoS based on
/// control delay (HCM Chapters 16 and 19).
///
/// Input: average control delay in seconds per vehicle.
///
/// | Grade | Delay (s/veh) | Description |
/// |-------|--------------|-------------|
/// | A     | <= 10        | Free flow   |
/// | B     | <= 20        | Stable flow |
/// | C     | <= 35        | Stable, acceptable delay |
/// | D     | <= 55        | Approaching unstable |
/// | E     | <= 80        | Unstable, significant delay |
/// | F     | > 80         | Forced flow / gridlock |
#[derive(Debug, Clone, Copy)]
pub struct VehicularUrbanStreet;

impl LosCriteria for VehicularUrbanStreet {
    fn classify(&self, delay_s: f64) -> LosGrade {
        if delay_s <= 10.0 {
            LosGrade::A
        } else if delay_s <= 20.0 {
            LosGrade::B
        } else if delay_s <= 35.0 {
            LosGrade::C
        } else if delay_s <= 55.0 {
            LosGrade::D
        } else if delay_s <= 80.0 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Vehicular Urban Street"
    }
    fn unit(&self) -> &'static str {
        "s/veh"
    }
}

// ---------------------------------------------------------------------------
// Unsignalized intersection (HCM Ch 20/21, delay-based)
// ---------------------------------------------------------------------------

/// Unsignalized intersection LoS based on control delay
/// (HCM Chapters 20 and 21).
///
/// Input: average control delay in seconds per vehicle.
///
/// | Grade | Delay (s/veh) | Description |
/// |-------|--------------|-------------|
/// | A     | <= 10        | Little or no delay |
/// | B     | <= 15        | Short delays |
/// | C     | <= 25        | Average delays |
/// | D     | <= 35        | Long delays |
/// | E     | <= 50        | Very long delays |
/// | F     | > 50         | Extreme delay / failure |
#[derive(Debug, Clone, Copy)]
pub struct VehicularUnsignalized;

impl LosCriteria for VehicularUnsignalized {
    fn classify(&self, delay_s: f64) -> LosGrade {
        if delay_s <= 10.0 {
            LosGrade::A
        } else if delay_s <= 15.0 {
            LosGrade::B
        } else if delay_s <= 25.0 {
            LosGrade::C
        } else if delay_s <= 35.0 {
            LosGrade::D
        } else if delay_s <= 50.0 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Vehicular Unsignalized Intersection"
    }
    fn unit(&self) -> &'static str {
        "s/veh"
    }
}

// ---------------------------------------------------------------------------
// Bicycle facility (HCM Ch 24, event-based)
// ---------------------------------------------------------------------------

/// Bicycle facility LoS based on hindrance events per minute
/// (HCM Chapter 24).
///
/// Events include meetings, passings, and active conflicts with
/// other bicyclists or pedestrians on shared-use paths.
///
/// | Grade | Events/min | Description |
/// |-------|-----------|-------------|
/// | A     | <= 10     | Few conflicts |
/// | B     | <= 20     | Occasional conflicts |
/// | C     | <= 30     | Frequent conflicts |
/// | D     | <= 40     | Significant conflicts |
/// | E     | <= 60     | Serious conflicts |
/// | F     | > 60      | Very crowded, breakdown |
#[derive(Debug, Clone, Copy)]
pub struct BicycleFacility;

impl LosCriteria for BicycleFacility {
    fn classify(&self, events_per_min: f64) -> LosGrade {
        if events_per_min <= 10.0 {
            LosGrade::A
        } else if events_per_min <= 20.0 {
            LosGrade::B
        } else if events_per_min <= 30.0 {
            LosGrade::C
        } else if events_per_min <= 40.0 {
            LosGrade::D
        } else if events_per_min <= 60.0 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Bicycle Facility"
    }
    fn unit(&self) -> &'static str {
        "events/min"
    }
}

// ---------------------------------------------------------------------------
// Transit capacity (TCQSM, load-factor-based)
// ---------------------------------------------------------------------------

/// Transit vehicle LoS based on passenger load factor
/// (Transit Capacity and Quality of Service Manual).
///
/// Load factor = passengers / seats. Values > 1.0 mean standees.
///
/// | Grade | Load factor | Description |
/// |-------|------------|-------------|
/// | A     | <= 0.50    | Many empty seats |
/// | B     | <= 0.75    | Some empty seats |
/// | C     | <= 1.00    | All seats occupied |
/// | D     | <= 1.25    | Comfortable standee load |
/// | E     | <= 1.50    | Maximum schedule load |
/// | F     | > 1.50     | Crush loading |
#[derive(Debug, Clone, Copy)]
pub struct TransitCapacity;

impl LosCriteria for TransitCapacity {
    fn classify(&self, load_factor: f64) -> LosGrade {
        if load_factor <= 0.50 {
            LosGrade::A
        } else if load_factor <= 0.75 {
            LosGrade::B
        } else if load_factor <= 1.00 {
            LosGrade::C
        } else if load_factor <= 1.25 {
            LosGrade::D
        } else if load_factor <= 1.50 {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        "Transit Capacity"
    }
    fn unit(&self) -> &'static str {
        "load factor"
    }
}

// ---------------------------------------------------------------------------
// Custom LoS criteria
// ---------------------------------------------------------------------------

/// Errors returned by custom LoS criteria validation.
#[derive(Debug, Clone, Copy, PartialEq, Error)]
pub enum LosCriteriaConfigError {
    #[error("thresholds must be in strictly ascending order")]
    NonAscendingThresholds,
}

/// User-defined LoS criteria with custom thresholds.
#[derive(Debug, Clone)]
pub struct CustomLosCriteria {
    name: &'static str,
    unit: &'static str,
    thresholds: [f64; 5],
}

impl CustomLosCriteria {
    /// Create custom LoS criteria.
    pub fn new(
        name: &'static str,
        unit: &'static str,
        thresholds: [f64; 5],
    ) -> Result<Self, LosCriteriaConfigError> {
        for i in 1..5 {
            if thresholds[i] <= thresholds[i - 1] {
                return Err(LosCriteriaConfigError::NonAscendingThresholds);
            }
        }
        Ok(Self {
            name,
            unit,
            thresholds,
        })
    }
}

impl LosCriteria for CustomLosCriteria {
    fn classify(&self, value: f64) -> LosGrade {
        if value <= self.thresholds[0] {
            LosGrade::A
        } else if value <= self.thresholds[1] {
            LosGrade::B
        } else if value <= self.thresholds[2] {
            LosGrade::C
        } else if value <= self.thresholds[3] {
            LosGrade::D
        } else if value <= self.thresholds[4] {
            LosGrade::E
        } else {
            LosGrade::F
        }
    }
    fn name(&self) -> &'static str {
        self.name
    }
    fn unit(&self) -> &'static str {
        self.unit
    }
}

// ---------------------------------------------------------------------------
// DensityStatistics
// ---------------------------------------------------------------------------

/// Summary statistics from a density grid.
#[derive(Debug, Clone)]
pub struct DensityStatistics {
    /// Maximum density across all cells (agents per cell area).
    pub max_density: f64,
    /// Mean density across occupied cells.
    pub mean_density: f64,
    /// Number of cells with at least one agent.
    pub occupied_cells: usize,
    /// Total number of agents counted.
    pub total_agents: usize,
    /// Worst LoS grade observed.
    pub worst_los: LosGrade,
    /// Count of cells at each LoS grade (A=0, B=1, ..., F=5).
    pub los_distribution: [usize; 6],
}

// ---------------------------------------------------------------------------
// DensityGrid
// ---------------------------------------------------------------------------

/// Errors returned by density-grid validation.
#[derive(Debug, Clone, Copy, PartialEq, Error)]
pub enum DensityGridError {
    #[error("extent must be positive")]
    InvalidExtent,
    #[error("cell_size must be positive")]
    InvalidCellSize,
}

/// Spatial density grid for computing agents-per-square-meter.
#[derive(Debug, Clone)]
pub struct DensityGrid {
    extent_x: f64,
    extent_y: f64,
    cell_size: f64,
    cell_area: f64,
    grid_w: usize,
    grid_h: usize,
    counts: Vec<u32>,
    total_agents: usize,
}

impl DensityGrid {
    /// Create a new density grid.
    pub fn new(extent_x: f64, extent_y: f64, cell_size: f64) -> Result<Self, DensityGridError> {
        if extent_x <= 0.0 || extent_y <= 0.0 {
            return Err(DensityGridError::InvalidExtent);
        }
        if cell_size <= 0.0 {
            return Err(DensityGridError::InvalidCellSize);
        }
        let grid_w = (extent_x / cell_size).ceil() as usize;
        let grid_h = (extent_y / cell_size).ceil() as usize;
        Ok(Self {
            extent_x,
            extent_y,
            cell_size,
            cell_area: cell_size * cell_size,
            grid_w,
            grid_h,
            counts: vec![0; grid_w * grid_h],
            total_agents: 0,
        })
    }

    /// Grid dimensions as `(width_cells, height_cells)`.
    pub fn grid_dimensions(&self) -> (usize, usize) {
        (self.grid_w, self.grid_h)
    }

    /// Cell side length in meters.
    pub fn cell_size(&self) -> f64 {
        self.cell_size
    }

    /// Space extent as `(extent_x, extent_y)`.
    pub fn extent(&self) -> (f64, f64) {
        (self.extent_x, self.extent_y)
    }

    /// Register an agent's position.
    ///
    /// Positions outside the extent are clamped to the nearest edge cell.
    pub fn add_position(&mut self, x: f64, y: f64) {
        let cx = ((x / self.cell_size).floor() as usize).min(self.grid_w.saturating_sub(1));
        let cy = ((y / self.cell_size).floor() as usize).min(self.grid_h.saturating_sub(1));
        let idx = cy * self.grid_w + cx;
        self.counts[idx] += 1;
        self.total_agents += 1;
    }

    /// Register multiple agent positions from an iterator of `(x, y)` pairs.
    pub fn add_positions(&mut self, positions: impl IntoIterator<Item = (f64, f64)>) {
        for (x, y) in positions {
            self.add_position(x, y);
        }
    }

    /// Agent count in the cell containing `(x, y)`.
    pub fn count_at(&self, x: f64, y: f64) -> u32 {
        let cx = ((x / self.cell_size).floor() as usize).min(self.grid_w.saturating_sub(1));
        let cy = ((y / self.cell_size).floor() as usize).min(self.grid_h.saturating_sub(1));
        self.counts[cy * self.grid_w + cx]
    }

    /// Density (agents per cell area) in the cell containing `(x, y)`.
    pub fn density_at(&self, x: f64, y: f64) -> f64 {
        self.count_at(x, y) as f64 / self.cell_area
    }

    /// LoS grade for the cell containing `(x, y)` using the given criteria.
    pub fn los_at(&self, x: f64, y: f64, criteria: &dyn LosCriteria) -> LosGrade {
        criteria.classify(self.density_at(x, y))
    }

    /// Density for a specific cell index `(cx, cy)`.
    pub fn density_at_cell(&self, cx: usize, cy: usize) -> f64 {
        if cx >= self.grid_w || cy >= self.grid_h {
            return 0.0;
        }
        self.counts[cy * self.grid_w + cx] as f64 / self.cell_area
    }

    /// Agent count for a specific cell index `(cx, cy)`.
    pub fn count_at_cell(&self, cx: usize, cy: usize) -> u32 {
        if cx >= self.grid_w || cy >= self.grid_h {
            return 0;
        }
        self.counts[cy * self.grid_w + cx]
    }

    /// Maximum density across all cells.
    pub fn max_density(&self) -> f64 {
        self.counts.iter().copied().max().unwrap_or(0) as f64 / self.cell_area
    }

    /// Mean density across occupied cells only.
    ///
    /// Returns 0.0 if no cells are occupied.
    pub fn mean_density_occupied(&self) -> f64 {
        let occupied: Vec<u32> = self.counts.iter().copied().filter(|&c| c > 0).collect();
        if occupied.is_empty() {
            return 0.0;
        }
        let sum: u32 = occupied.iter().sum();
        (sum as f64 / occupied.len() as f64) / self.cell_area
    }

    /// Mean density across all cells.
    pub fn mean_density_all(&self) -> f64 {
        let total_area = self.grid_w as f64 * self.grid_h as f64 * self.cell_area;
        if total_area == 0.0 {
            return 0.0;
        }
        self.total_agents as f64 / total_area
    }

    /// Compute full statistics using the given LoS criteria.
    pub fn statistics(&self, criteria: &dyn LosCriteria) -> DensityStatistics {
        let mut max_count: u32 = 0;
        let mut occupied_cells = 0usize;
        let mut los_distribution = [0usize; 6];

        for &count in &self.counts {
            if count > 0 {
                occupied_cells += 1;
                if count > max_count {
                    max_count = count;
                }
                let density = count as f64 / self.cell_area;
                let los = criteria.classify(density);
                let idx = match los {
                    LosGrade::A => 0,
                    LosGrade::B => 1,
                    LosGrade::C => 2,
                    LosGrade::D => 3,
                    LosGrade::E => 4,
                    LosGrade::F => 5,
                };
                los_distribution[idx] += 1;
            }
        }

        let max_density = max_count as f64 / self.cell_area;
        let worst_los = criteria.classify(max_density);

        let mean_density = if occupied_cells > 0 {
            (self.total_agents as f64 / occupied_cells as f64) / self.cell_area
        } else {
            0.0
        };

        DensityStatistics {
            max_density,
            mean_density,
            occupied_cells,
            total_agents: self.total_agents,
            worst_los,
            los_distribution,
        }
    }

    /// Raw counts slice (row-major: `counts[cy * grid_w + cx]`).
    pub fn counts(&self) -> &[u32] {
        &self.counts
    }

    /// Total number of agents registered.
    pub fn total_agents(&self) -> usize {
        self.total_agents
    }

    /// Clear all counts for the next time step.
    pub fn clear(&mut self) {
        self.counts.fill(0);
        self.total_agents = 0;
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn pedestrian_walkway_thresholds() {
        let c = PedestrianWalkway;
        assert_eq!(c.classify(0.0), LosGrade::A);
        assert_eq!(c.classify(0.31), LosGrade::A);
        assert_eq!(c.classify(0.32), LosGrade::B);
        assert_eq!(c.classify(0.43), LosGrade::B);
        assert_eq!(c.classify(0.44), LosGrade::C);
        assert_eq!(c.classify(0.72), LosGrade::C);
        assert_eq!(c.classify(0.73), LosGrade::D);
        assert_eq!(c.classify(1.08), LosGrade::D);
        assert_eq!(c.classify(1.09), LosGrade::E);
        assert_eq!(c.classify(2.17), LosGrade::E);
        assert_eq!(c.classify(2.18), LosGrade::F);
        assert_eq!(c.classify(10.0), LosGrade::F);
    }

    #[test]
    fn pedestrian_stairway_thresholds() {
        let c = PedestrianStairway;
        assert_eq!(c.classify(0.5), LosGrade::A);
        assert_eq!(c.classify(0.55), LosGrade::B);
        assert_eq!(c.classify(0.73), LosGrade::C);
        assert_eq!(c.classify(1.1), LosGrade::D);
        assert_eq!(c.classify(1.55), LosGrade::E);
        assert_eq!(c.classify(3.0), LosGrade::F);
    }

    #[test]
    fn pedestrian_queuing_thresholds() {
        let c = PedestrianQueuing;
        assert_eq!(c.classify(0.5), LosGrade::A);
        assert_eq!(c.classify(0.9), LosGrade::B);
        assert_eq!(c.classify(1.2), LosGrade::C);
        assert_eq!(c.classify(2.0), LosGrade::D);
        assert_eq!(c.classify(4.0), LosGrade::E);
        assert_eq!(c.classify(6.0), LosGrade::F);
    }

    #[test]
    fn vehicular_freeway_thresholds() {
        let c = VehicularFreeway;
        assert_eq!(c.classify(5.0), LosGrade::A);
        assert_eq!(c.classify(9.0), LosGrade::B);
        assert_eq!(c.classify(14.0), LosGrade::C);
        assert_eq!(c.classify(20.0), LosGrade::D);
        assert_eq!(c.classify(26.0), LosGrade::E);
        assert_eq!(c.classify(35.0), LosGrade::F);
    }

    #[test]
    fn vehicular_urban_street_thresholds() {
        let c = VehicularUrbanStreet;
        assert_eq!(c.classify(5.0), LosGrade::A);
        assert_eq!(c.classify(15.0), LosGrade::B);
        assert_eq!(c.classify(30.0), LosGrade::C);
        assert_eq!(c.classify(45.0), LosGrade::D);
        assert_eq!(c.classify(70.0), LosGrade::E);
        assert_eq!(c.classify(100.0), LosGrade::F);
    }

    #[test]
    fn vehicular_unsignalized_thresholds() {
        let c = VehicularUnsignalized;
        assert_eq!(c.classify(5.0), LosGrade::A);
        assert_eq!(c.classify(12.0), LosGrade::B);
        assert_eq!(c.classify(20.0), LosGrade::C);
        assert_eq!(c.classify(30.0), LosGrade::D);
        assert_eq!(c.classify(45.0), LosGrade::E);
        assert_eq!(c.classify(60.0), LosGrade::F);
    }

    #[test]
    fn bicycle_facility_thresholds() {
        let c = BicycleFacility;
        assert_eq!(c.classify(5.0), LosGrade::A);
        assert_eq!(c.classify(15.0), LosGrade::B);
        assert_eq!(c.classify(25.0), LosGrade::C);
        assert_eq!(c.classify(35.0), LosGrade::D);
        assert_eq!(c.classify(55.0), LosGrade::E);
        assert_eq!(c.classify(70.0), LosGrade::F);
    }

    #[test]
    fn transit_capacity_thresholds() {
        let c = TransitCapacity;
        assert_eq!(c.classify(0.3), LosGrade::A);
        assert_eq!(c.classify(0.6), LosGrade::B);
        assert_eq!(c.classify(0.9), LosGrade::C);
        assert_eq!(c.classify(1.1), LosGrade::D);
        assert_eq!(c.classify(1.4), LosGrade::E);
        assert_eq!(c.classify(2.0), LosGrade::F);
    }

    #[test]
    fn custom_criteria() {
        let c = CustomLosCriteria::new("Test", "units", [1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        assert_eq!(c.classify(0.5), LosGrade::A);
        assert_eq!(c.classify(1.5), LosGrade::B);
        assert_eq!(c.classify(2.5), LosGrade::C);
        assert_eq!(c.classify(3.5), LosGrade::D);
        assert_eq!(c.classify(4.5), LosGrade::E);
        assert_eq!(c.classify(5.5), LosGrade::F);
        assert_eq!(c.name(), "Test");
        assert_eq!(c.unit(), "units");
    }

    #[test]
    fn custom_criteria_bad_order() {
        let err = CustomLosCriteria::new("Bad", "x", [1.0, 3.0, 2.0, 4.0, 5.0]).unwrap_err();
        assert_eq!(err, LosCriteriaConfigError::NonAscendingThresholds);
    }

    #[test]
    fn basic_density() {
        let mut grid = DensityGrid::new(10.0, 10.0, 1.0).unwrap();
        grid.add_position(0.5, 0.5);
        grid.add_position(0.1, 0.9);
        grid.add_position(0.3, 0.3);
        grid.add_position(0.9, 0.1);

        assert_eq!(grid.count_at(0.5, 0.5), 4);
        assert!((grid.density_at(0.5, 0.5) - 4.0).abs() < 1e-9);
        assert_eq!(grid.los_at(0.5, 0.5, &PedestrianWalkway), LosGrade::F);
    }

    #[test]
    fn empty_cell_density() {
        let grid = DensityGrid::new(10.0, 10.0, 1.0).unwrap();
        assert_eq!(grid.count_at(5.0, 5.0), 0);
        assert!((grid.density_at(5.0, 5.0)).abs() < 1e-9);
        assert_eq!(grid.los_at(5.0, 5.0, &PedestrianWalkway), LosGrade::A);
    }

    #[test]
    fn statistics_with_walkway_criteria() {
        let mut grid = DensityGrid::new(10.0, 10.0, 2.0).unwrap();
        // 2 agents in one cell: density = 2/4 = 0.5 pax/m2 -> LoS C (walkway)
        grid.add_position(1.0, 1.0);
        grid.add_position(1.5, 1.5);
        // 1 agent in another cell: density = 1/4 = 0.25 pax/m2 -> LoS A
        grid.add_position(5.0, 5.0);

        let stats = grid.statistics(&PedestrianWalkway);
        assert_eq!(stats.total_agents, 3);
        assert_eq!(stats.occupied_cells, 2);
        assert!((stats.max_density - 0.5).abs() < 1e-9);
        assert_eq!(stats.worst_los, LosGrade::C);
        assert_eq!(stats.los_distribution[0], 1); // A
        assert_eq!(stats.los_distribution[2], 1); // C
    }

    #[test]
    fn statistics_with_queuing_criteria() {
        let mut grid = DensityGrid::new(10.0, 10.0, 2.0).unwrap();
        // Same data as above, but with queuing criteria:
        // 0.5 pax/m2 -> LoS A for queuing (threshold is 0.83)
        grid.add_position(1.0, 1.0);
        grid.add_position(1.5, 1.5);
        grid.add_position(5.0, 5.0);

        let stats = grid.statistics(&PedestrianQueuing);
        assert_eq!(stats.worst_los, LosGrade::A); // 0.5 < 0.83
        assert_eq!(stats.los_distribution[0], 2); // both cells are A
    }

    #[test]
    fn different_criteria_same_density() {
        // 0.6 pax/m2 is:
        //   Walkway: C (0.43 < 0.6 <= 0.72)
        //   Stairway: B (0.54 < 0.6 <= 0.72)
        //   Queuing: A (0.6 <= 0.83)
        assert_eq!(PedestrianWalkway.classify(0.6), LosGrade::C);
        assert_eq!(PedestrianStairway.classify(0.6), LosGrade::B);
        assert_eq!(PedestrianQueuing.classify(0.6), LosGrade::A);
    }

    #[test]
    fn clear_resets() {
        let mut grid = DensityGrid::new(10.0, 10.0, 1.0).unwrap();
        grid.add_position(0.5, 0.5);
        assert_eq!(grid.total_agents(), 1);
        grid.clear();
        assert_eq!(grid.total_agents(), 0);
        assert_eq!(grid.count_at(0.5, 0.5), 0);
    }

    #[test]
    fn add_positions_batch() {
        let mut grid = DensityGrid::new(10.0, 10.0, 1.0).unwrap();
        let positions = vec![(0.5, 0.5), (0.1, 0.1), (5.0, 5.0)];
        grid.add_positions(positions);
        assert_eq!(grid.total_agents(), 3);
        assert_eq!(grid.count_at(0.5, 0.5), 2);
        assert_eq!(grid.count_at(5.0, 5.0), 1);
    }

    #[test]
    fn clamping_out_of_bounds() {
        let mut grid = DensityGrid::new(10.0, 10.0, 1.0).unwrap();
        grid.add_position(-5.0, -5.0);
        assert_eq!(grid.count_at(0.0, 0.0), 1);
        grid.add_position(100.0, 100.0);
        assert_eq!(grid.count_at(9.9, 9.9), 1);
    }

    #[test]
    fn los_grade_ordering() {
        assert!(LosGrade::A < LosGrade::B);
        assert!(LosGrade::B < LosGrade::C);
        assert!(LosGrade::E < LosGrade::F);
    }

    #[test]
    fn los_grade_display() {
        assert_eq!(format!("{}", LosGrade::A), "A");
        assert_eq!(format!("{}", LosGrade::F), "F");
    }

    #[test]
    fn criteria_metadata() {
        assert_eq!(PedestrianWalkway.name(), "Pedestrian Walkway");
        assert_eq!(PedestrianWalkway.unit(), "pax/m2");
        assert_eq!(VehicularUrbanStreet.name(), "Vehicular Urban Street");
        assert_eq!(VehicularUrbanStreet.unit(), "s/veh");
        assert_eq!(TransitCapacity.name(), "Transit Capacity");
        assert_eq!(TransitCapacity.unit(), "load factor");
    }
}