grid1d 0.5.0

//! Classification of a coordinate slice by grid interval.
//!
//! Given a set of coordinates and a 1D grid (any type implementing
//! [`FindIntervalIdOfPoint`] + [`HasIntervalIdRange`]), every coordinate is mapped to
//! either the grid interval that contains it or to an *outside-domain* set when no
//! interval contains it.  The result is captured in a
//! [`CoordsByIntervalClassification`] value that provides ergonomic accessors for
//! downstream processing (finite-element assembly, interpolation, adaptive
//! refinement, etc.).
//!
//! # Accepted coordinate inputs
//!
//! The slice-based constructors accept `coords: &[T]`.  Because [`Coords1D<T>`]
//! implements `Deref<Target = [T]>`, a `&Coords1D<T>` coerces automatically and
//! can be passed anywhere a `&[T]` is expected — no explicit `.as_ref()` or
//! `.deref()` call is needed:
//!
//! ```ignore
//! let classification =
//!     CoordsByIntervalClassification::try_new_dense_from_coords(&my_coords1d, &grid)?;
//! //                                                             ^^^^^^^^^^^^^
//! //                   &Coords1D<f64> coerces to &[f64] via Deref
//! ```
//!
//! # Storage strategies
//!
//! Two internal representations are offered for the inside-interval mapping:
//!
//! | Strategy | Constructor | Best for |
//! |----------|-------------|----------|
//! | **Dense** | [`CoordsByIntervalClassification::try_new_dense_from_coords`] | Most intervals are occupied |
//! | **Sparse** | [`CoordsByIntervalClassification::try_new_sparse_from_coords`] | Few intervals are occupied |
//!
//! Both strategies guarantee that the coordinate-index sets stored in each bucket
//! are **sorted in ascending order**, so callers can rely on that property without
//! an extra sorting step.
//!
//! # Parallel variants
//!
//! Every sequential constructor has a `_par` counterpart
//! ([`try_new_dense_from_coords_par`](CoordsByIntervalClassification::try_new_dense_from_coords_par),
//! [`try_new_sparse_from_coords_par`](CoordsByIntervalClassification::try_new_sparse_from_coords_par))
//! that classifies coordinates in parallel using [Rayon].  The grid type must
//! implement `Sync`; the coordinate type must implement `Send + Sync`.
//!
//! # Sorted-input fast path
//!
//! When the input coordinates are already sorted (wrapped in a [`Coords1D`]),
//! [`CoordsByIntervalClassification::new_dense_from_coords_sorted`] performs a
//! single linear sweep through both the coordinate slice and the interval list,
//! avoiding a per-coordinate binary-search lookup and returning `Self` directly
//! (infallible, because [`Coords1D`] is guaranteed non-empty).
//!
//! [Rayon]: https://docs.rs/rayon

#![deny(rustdoc::broken_intra_doc_links)]

use crate::{
    PositiveNumPoints1D,
    coords::Coords1D,
    grids::traits::{FindIntervalIdOfPoint, HasIntervalIdRange, HasIntervals},
    intervals::{Contains, IntervalBoundsRuntime},
    scalars::{IntervalId, NumIntervals},
};
use num_valid::core::errors::capture_backtrace;
use rayon::iter::{IntoParallelIterator, IntoParallelRefIterator, ParallelIterator};
use serde::{Deserialize, Serialize};
use sorted_vec::SortedSet;
use std::{backtrace::Backtrace, collections::BTreeMap};
use thiserror::Error;
use try_create::TryNew;

/// Errors that can occur when classifying coordinates by grid interval.
#[derive(Debug, Error)]
pub enum ErrorsCoordsByIntervalClassification {
    #[error("Cannot classify points by interval: no coordinates provided")]
    /// No coordinates were supplied to the classification function.
    EmptyCoords {
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },

    #[error("Cannot classify point: interval with id {interval_id} not present in grid partition")]
    /// The requested interval ID does not exist in the grid partition.
    IntervalNotPresent {
        /// The interval ID that was not found.
        interval_id: IntervalId,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },
}

/// Stores the mapping from each grid interval to the coordinate indices that fall inside it.
///
/// Two storage strategies are available:
/// - [`Dense`](CoordIdsByInterval::Dense): a contiguous `Vec` indexed by [`IntervalId`];
///   slots for empty intervals are `None`. Efficient when most intervals are occupied.
/// - [`Sparse`](CoordIdsByInterval::Sparse): a [`BTreeMap`] keyed by [`IntervalId`];
///   absent keys mean empty intervals. Efficient when only a few intervals are occupied.
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
pub enum CoordIdsByInterval {
    /// Dense storage: `buckets[i]` holds the sorted coordinate indices for interval `i`,
    /// or `None` if no coordinates fall in that interval.
    Dense {
        /// Per-interval buckets; index equals the [`IntervalId`] value.
        buckets: Vec<Option<SortedSet<usize>>>,
    },
    /// Sparse storage: only intervals that contain at least one coordinate appear as keys.
    Sparse {
        /// Map from [`IntervalId`] to sorted coordinate indices.
        buckets: BTreeMap<IntervalId, SortedSet<usize>>,
    },
}

impl CoordIdsByInterval {
    /// Returns the sorted set of coordinate indices inside the given interval, if any.
    ///
    /// Returns `None` if no coordinates were classified into `interval_id`
    /// (either the interval is empty or the id is out of range for a dense storage).
    #[inline]
    pub fn point_ids_in_interval(&self, interval_id: &IntervalId) -> Option<&SortedSet<usize>> {
        match self {
            CoordIdsByInterval::Dense { buckets } => {
                let bucket_index = *interval_id.as_ref();
                buckets.get(bucket_index)?.as_ref()
            }
            CoordIdsByInterval::Sparse { buckets } => buckets.get(interval_id),
        }
    }

    /// Returns all non-empty intervals together with their coordinate indices, in ascending [`IntervalId`] order.
    ///
    /// Each element of the returned `Vec` is a pair `(interval_id, coord_ids)` where
    /// `coord_ids` is the sorted set of coordinate indices that fall inside that interval.
    /// Intervals with no coordinates assigned are omitted.
    #[must_use]
    pub fn point_ids(&self) -> Vec<(IntervalId, &SortedSet<usize>)> {
        match self {
            CoordIdsByInterval::Dense { buckets } => buckets
                .iter()
                .enumerate()
                .filter_map(|(i, opt_bucket)| {
                    opt_bucket
                        .as_ref()
                        .map(|bucket| (IntervalId::new(i), bucket))
                })
                .collect(),
            CoordIdsByInterval::Sparse { buckets } => {
                buckets.iter().map(|(id, bucket)| (*id, bucket)).collect()
            }
        }
    }
}

/// Classification of a coordinate slice by grid interval.
///
/// Given a slice of coordinates and a grid, every coordinate is assigned either to
/// the grid interval that contains it, or to an *outside-domain* set if no interval
/// contains it.
///
/// The inside mapping can be stored in two ways, selected at construction time:
/// - **Dense** ([`try_new_dense_from_coords`](CoordsByIntervalClassification::try_new_dense_from_coords),
///   [`new_dense_from_coords_sorted`](CoordsByIntervalClassification::new_dense_from_coords_sorted)):
///   a `Vec` indexed by [`IntervalId`]; best when most intervals are occupied.
/// - **Sparse** ([`try_new_sparse_from_coords`](CoordsByIntervalClassification::try_new_sparse_from_coords)):
///   a [`BTreeMap`] keyed by [`IntervalId`]; best when only few intervals are occupied.
///
/// Parallel variants (`_par` suffix) classify coordinates using Rayon and require
/// the grid type to implement `Sync`.
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
pub struct CoordsByIntervalClassification {
    /// Mapping from each occupied grid interval to the sorted set of coordinate indices
    /// that fall inside it. The representation (dense or sparse) is chosen at
    /// construction time and is fixed for the lifetime of this value.
    inside: CoordIdsByInterval,
    /// Sorted set of coordinate indices that did not belong to any interval in the
    /// partition, i.e. they lie strictly outside the grid domain.
    outside_domain: SortedSet<usize>,
}

/// Mutable accumulator used while building a [`CoordIdsByInterval::Dense`] result.
///
/// `buckets[i]` is `None` until the first coordinate is classified into interval `i`;
/// after that it is a `Vec<usize>` grown in ascending `coord_id` order so that the
/// final [`SortedSet`] can be constructed without sorting.
struct DenseData {
    /// Per-interval buckets, indexed directly by [`IntervalId`] as a `usize`.
    buckets: Vec<Option<Vec<usize>>>,
    /// Coordinate indices that fell outside the domain, appended in ascending order.
    outside_domain: Vec<usize>,
    /// Pre-computed per-bucket capacity hint derived from `num_points / num_intervals`.
    hint_num_points_per_bucket: Option<PositiveNumPoints1D>,
}

/// Mutable accumulator used while building a [`CoordIdsByInterval::Sparse`] result.
///
/// Only intervals that receive at least one coordinate are inserted as keys, so the
/// map stays compact for grids with many empty intervals.
struct SparseData {
    /// Per-interval buckets; only occupied intervals appear as keys.
    buckets: BTreeMap<IntervalId, Vec<usize>>,
    /// Coordinate indices that fell outside the domain, appended in ascending order.
    outside_domain: Vec<usize>,
    /// Pre-computed per-bucket capacity hint derived from `num_points / num_intervals`.
    hint_num_points_per_bucket: Option<PositiveNumPoints1D>,
}

/// Sealed interface shared by [`DenseData`] and [`SparseData`].
///
/// A single classification pass calls only `push_inside` / `push_outside` and
/// then `finalize` (sequential) or `finalize_par` (parallel finalisation),
/// keeping the pass generic over the chosen storage strategy.
trait BucketsTrait {
    /// Allocate an empty accumulator pre-sized for `num_intervals` intervals and
    /// `hint_num_points` total coordinates.
    fn create_empty(num_intervals: NumIntervals, hint_num_points: PositiveNumPoints1D) -> Self;

    /// Record that coordinate `coord_id` falls inside `interval_id`.
    ///
    /// Implementations must append `coord_id` in non-decreasing order within each
    /// bucket so that [`finalize`](BucketsTrait::finalize) can call
    /// `SortedSet::from_sorted` without extra sorting.
    fn push_inside(&mut self, interval_id: IntervalId, coord_id: usize);

    /// Record that coordinate `coord_id` lies outside the grid domain.
    fn push_outside(&mut self, coord_id: usize);

    /// Consume the accumulator and produce the final [`CoordsByIntervalClassification`],
    /// converting the inner `Vec` buckets to [`SortedSet`]s sequentially.
    fn finalize(self) -> CoordsByIntervalClassification;

    /// Like [`finalize`](BucketsTrait::finalize) but converts the buckets in parallel
    /// using Rayon. Use after the two-phase parallel classification path.
    fn finalize_par(self) -> CoordsByIntervalClassification;

    /// Helper: allocate the `outside_domain` vec with capacity `hint_num_points` and
    /// derive the per-bucket point hint. Both outputs are needed by every `create_empty`
    /// implementation, so the logic lives here to avoid duplication.
    #[inline(always)]
    fn allocate_outside_domain_vec(
        num_intervals: NumIntervals,
        hint_num_points: PositiveNumPoints1D,
    ) -> (Vec<usize>, Option<PositiveNumPoints1D>) {
        (
            Vec::with_capacity(*hint_num_points.as_ref()),
            compute_hint_num_points_per_interval(hint_num_points, num_intervals),
        )
    }
}

/// Compute a per-interval capacity hint as `ceil(num_points / num_intervals)`.
///
/// The `+1` guards against zero-capacity allocations when `num_points < num_intervals`.
/// Returns `Some` so the result can be stored directly in the `Option` field of the
/// accumulator structs.
#[inline(always)]
fn compute_hint_num_points_per_interval(
    num_points: PositiveNumPoints1D,
    num_intervals: NumIntervals,
) -> Option<PositiveNumPoints1D> {
    Some(
        PositiveNumPoints1D::try_new((num_points.into_inner() / num_intervals.into_inner()) + 1)
            .unwrap(),
    )
}

impl BucketsTrait for DenseData {
    fn push_inside(&mut self, interval_id: IntervalId, coord_id: usize) {
        let bucket_index = *interval_id.as_ref();
        if bucket_index >= self.buckets.len() {
            self.buckets.resize(bucket_index + 1, None);
        }
        self.buckets[bucket_index]
            .get_or_insert_with(|| match self.hint_num_points_per_bucket {
                Some(hint) => Vec::with_capacity(*hint.as_ref()),
                None => Vec::new(),
            })
            .push(coord_id);
    }

    fn push_outside(&mut self, coord_id: usize) {
        self.outside_domain.push(coord_id);
    }

    fn create_empty(num_intervals: NumIntervals, hint_num_points: PositiveNumPoints1D) -> Self {
        let buckets = Vec::<Option<Vec<usize>>>::with_capacity(*num_intervals.as_ref());

        let (outside_domain, hint_num_points_per_bucket) =
            Self::allocate_outside_domain_vec(num_intervals, hint_num_points);

        Self {
            buckets,
            outside_domain,
            hint_num_points_per_bucket,
        }
    }

    fn finalize(self) -> CoordsByIntervalClassification {
        let buckets = self
            .buckets
            .into_iter()
            // SAFETY: coord_ids are added via enumerate() in ascending order, so each bucket is sorted.
            .map(|opt| opt.map(|bucket| unsafe { SortedSet::from_sorted(bucket) }))
            .collect::<Vec<Option<SortedSet<usize>>>>();

        let outside_domain = unsafe { SortedSet::from_sorted(self.outside_domain) };

        CoordsByIntervalClassification {
            inside: CoordIdsByInterval::Dense { buckets },
            outside_domain,
        }
    }

    fn finalize_par(self) -> CoordsByIntervalClassification {
        let buckets = self
            .buckets
            .into_par_iter()
            // SAFETY: coord_ids are added via enumerate() in ascending order, so each bucket is sorted.
            .map(|opt| opt.map(|bucket| unsafe { SortedSet::from_sorted(bucket) }))
            .collect::<Vec<Option<SortedSet<usize>>>>();

        let outside_domain = unsafe { SortedSet::from_sorted(self.outside_domain) };

        CoordsByIntervalClassification {
            inside: CoordIdsByInterval::Dense { buckets },
            outside_domain,
        }
    }
}

impl BucketsTrait for SparseData {
    fn create_empty(num_intervals: NumIntervals, hint_num_points: PositiveNumPoints1D) -> Self {
        let (outside_domain, hint_num_points_per_bucket) =
            Self::allocate_outside_domain_vec(num_intervals, hint_num_points);

        Self {
            buckets: BTreeMap::new(),
            outside_domain,
            hint_num_points_per_bucket,
        }
    }

    fn push_inside(&mut self, interval_id: IntervalId, coord_id: usize) {
        self.buckets
            .entry(interval_id)
            .or_insert_with(|| match self.hint_num_points_per_bucket {
                Some(hint) => Vec::with_capacity(*hint.as_ref()),
                None => Vec::new(),
            })
            .push(coord_id);
    }

    fn push_outside(&mut self, coord_id: usize) {
        self.outside_domain.push(coord_id);
    }

    fn finalize(self) -> CoordsByIntervalClassification {
        let buckets = self
            .buckets
            .into_iter()
            // SAFETY: We ensure that each bucket is sorted by construction, since we push coord_ids in order.
            .map(|(id, bucket)| (id, unsafe { SortedSet::from_sorted(bucket) }))
            .collect::<BTreeMap<IntervalId, SortedSet<usize>>>();

        let outside_domain = unsafe { SortedSet::from_sorted(self.outside_domain) };

        CoordsByIntervalClassification {
            inside: CoordIdsByInterval::Sparse { buckets },
            outside_domain,
        }
    }

    fn finalize_par(self) -> CoordsByIntervalClassification {
        let buckets = self
            .buckets
            .into_par_iter()
            // SAFETY: We ensure that each bucket is sorted by construction, since we push coord_ids in order.
            .map(|(id, bucket)| (id, unsafe { SortedSet::from_sorted(bucket) }))
            .collect::<BTreeMap<IntervalId, SortedSet<usize>>>();

        let outside_domain = unsafe { SortedSet::from_sorted(self.outside_domain) };

        CoordsByIntervalClassification {
            inside: CoordIdsByInterval::Sparse { buckets },
            outside_domain,
        }
    }
}

impl CoordsByIntervalClassification {
    /// Returns the sorted set of coordinate indices that fall outside the grid domain.
    ///
    /// These are indices of coordinates that did not belong to any interval
    /// in the partition (i.e. they lie outside the domain).
    #[inline(always)]
    pub fn point_ids_outside(&self) -> &SortedSet<usize> {
        &self.outside_domain
    }

    /// Returns the sorted set of coordinate indices inside the given interval, or an error
    /// if no coordinates were classified into `interval_id`.
    ///
    /// This is the fallible counterpart of
    /// [`CoordIdsByInterval::point_ids_in_interval`]: instead of returning `None` when
    /// the interval is absent or empty, it returns
    /// [`ErrorsCoordsByIntervalClassification::IntervalNotPresent`], which carries the
    /// offending [`IntervalId`] and a backtrace.
    ///
    /// Prefer this method when the caller considers an empty interval to be a logic
    /// error and wants to propagate a descriptive error rather than handle `None`
    /// explicitly.
    ///
    /// # Errors
    /// Returns [`ErrorsCoordsByIntervalClassification::IntervalNotPresent`] when:
    /// - **Dense storage**: the bucket slot for `interval_id` is `None` (no coordinates
    ///   were classified there) or `interval_id` is beyond the allocated bucket range.
    /// - **Sparse storage**: `interval_id` has no entry in the map.
    #[inline(always)]
    pub fn try_point_ids_in_interval(
        &self,
        interval_id: &IntervalId,
    ) -> Result<&SortedSet<usize>, ErrorsCoordsByIntervalClassification> {
        self.inside
            .point_ids_in_interval(interval_id)
            .ok_or_else(
                || ErrorsCoordsByIntervalClassification::IntervalNotPresent {
                    interval_id: *interval_id,
                    backtrace: capture_backtrace(),
                },
            )
    }

    /// Returns all non-empty intervals together with their coordinate indices, in ascending [`IntervalId`] order.
    ///
    /// Each element of the returned `Vec` is a pair `(interval_id, coord_ids)` where
    /// `coord_ids` is the sorted set of coordinate indices that fall inside that interval.
    /// Intervals with no coordinates assigned are omitted.
    #[inline(always)]
    pub fn point_ids_inside(&self) -> Vec<(IntervalId, &SortedSet<usize>)> {
        self.inside.point_ids()
    }

    /// Generic sequential classification kernel.
    ///
    /// Iterates `coords` once, calling [`FindIntervalIdOfPoint::find_interval_id_of_point`]
    /// for each element and routing the result to either
    /// [`BucketsTrait::push_inside`] or [`BucketsTrait::push_outside`]. The
    /// storage strategy (dense or sparse) is selected by the `DataType` type
    /// parameter.
    fn try_new_from_coords<G, DataType>(
        coords: &[G::Point1DType],
        grid: &G,
    ) -> Result<Self, ErrorsCoordsByIntervalClassification>
    where
        G: FindIntervalIdOfPoint + HasIntervalIdRange,
        DataType: BucketsTrait,
    {
        if coords.is_empty() {
            return Err(ErrorsCoordsByIntervalClassification::EmptyCoords {
                backtrace: capture_backtrace(),
            });
        }

        let num_coords = PositiveNumPoints1D::try_new(coords.len()).unwrap();
        let mut data = DataType::create_empty(grid.num_intervals(), num_coords);

        coords.iter().enumerate().for_each(|(coord_id, coord)| {
            match grid.find_interval_id_of_point(coord) {
                Some(interval_id) => {
                    data.push_inside(interval_id, coord_id);
                }
                None => {
                    data.push_outside(coord_id);
                }
            }
        });

        Ok(data.finalize())
    }

    /// Classifies `coords` against `grid` using **dense** storage.
    ///
    /// Each coordinate is located in the grid and its index placed into the bucket
    /// for the containing interval, or into the outside-domain set if it lies outside
    /// the grid domain. Coordinates are processed sequentially.
    ///
    /// `coords` accepts any `&[T]`, including a `&Coords1D<T>` which coerces
    /// automatically via `Deref`.
    ///
    /// # Arguments
    /// - `coords` – slice of points to classify (or a `&Coords1D<T>` via deref coercion).
    /// - `grid`   – the partition used to locate each point.
    ///
    /// # Errors
    /// Returns [`ErrorsCoordsByIntervalClassification::EmptyCoords`] when `coords` is empty.
    pub fn try_new_dense_from_coords<G>(
        coords: &[G::Point1DType],
        grid: &G,
    ) -> Result<Self, ErrorsCoordsByIntervalClassification>
    where
        G: FindIntervalIdOfPoint + HasIntervalIdRange,
    {
        Self::try_new_from_coords::<G, DenseData>(coords, grid)
    }

    /// Classifies `coords` against `grid` using **sparse** storage.
    ///
    /// Like [`try_new_dense_from_coords`](Self::try_new_dense_from_coords) but stores the
    /// inside mapping in a [`BTreeMap`]: only intervals that actually contain at least one
    /// coordinate appear as keys. Prefer this variant when the grid has many intervals
    /// but only a small subset are occupied.
    ///
    /// `coords` accepts any `&[T]`, including a `&Coords1D<T>` which coerces
    /// automatically via `Deref`.
    ///
    /// # Arguments
    /// - `coords` – slice of points to classify (or a `&Coords1D<T>` via deref coercion).
    /// - `grid`   – the partition used to locate each point.
    ///
    /// # Errors
    /// Returns [`ErrorsCoordsByIntervalClassification::EmptyCoords`] when `coords` is empty.
    pub fn try_new_sparse_from_coords<G>(
        coords: &[G::Point1DType],
        grid: &G,
    ) -> Result<Self, ErrorsCoordsByIntervalClassification>
    where
        G: FindIntervalIdOfPoint + HasIntervalIdRange,
    {
        Self::try_new_from_coords::<G, SparseData>(coords, grid)
    }

    /// Generic two-phase parallel classification kernel.
    ///
    /// **Phase 1 (parallel):** each coordinate is mapped to `Option<IntervalId>` in
    /// parallel via Rayon. Because `par_iter()` on a slice is an
    /// `IndexedParallelIterator`, `collect()` preserves input order:
    /// `classified[i]` always corresponds to `coords[i]`.
    ///
    /// **Phase 2 (sequential):** the ordered `classified` vector is iterated once
    /// to fill the accumulator in ascending `coord_id` order, satisfying the
    /// sorted-bucket invariant required by `SortedSet::from_sorted`.
    ///
    /// The storage strategy (dense or sparse) is selected by the `DataType` type
    /// parameter.
    fn try_new_from_coords_par<G, DataType>(
        coords: &[G::Point1DType],
        grid: &G,
    ) -> Result<Self, ErrorsCoordsByIntervalClassification>
    where
        G: FindIntervalIdOfPoint + HasIntervalIdRange + Sync,
        G::Point1DType: Send + Sync,
        DataType: BucketsTrait + Send,
    {
        if coords.is_empty() {
            return Err(ErrorsCoordsByIntervalClassification::EmptyCoords {
                backtrace: capture_backtrace(),
            });
        }

        // Phase 1 (parallel): classify each coord. par_iter() on a slice is an IndexedParallelIterator,
        // so collect() preserves input order: classified[i] corresponds to coords[i].
        let classified: Vec<Option<IntervalId>> = coords
            .par_iter()
            .map(|coord| grid.find_interval_id_of_point(coord))
            .collect();

        // Phase 2 (sequential): iterate in coord_id order so each bucket stays sorted.
        let num_coords = PositiveNumPoints1D::try_new(coords.len()).unwrap();
        let mut data = DataType::create_empty(grid.num_intervals(), num_coords);

        for (coord_id, maybe_interval) in classified.into_iter().enumerate() {
            match maybe_interval {
                Some(interval_id) => {
                    data.push_inside(interval_id, coord_id);
                }
                None => data.push_outside(coord_id),
            }
        }

        Ok(data.finalize_par())
    }

    /// Classifies `coords` against `grid` using **dense** storage and **parallel** execution.
    ///
    /// Phase 1 maps each coordinate to its interval in parallel via Rayon, preserving
    /// input order. Phase 2 fills the dense buckets sequentially so that each bucket
    /// remains sorted by coordinate index.
    ///
    /// Prefer [`try_new_dense_from_coords`](Self::try_new_dense_from_coords) for small
    /// inputs where parallelism overhead outweighs the benefit.
    ///
    /// # Errors
    /// Returns [`ErrorsCoordsByIntervalClassification::EmptyCoords`] when `coords` is empty.
    pub fn try_new_dense_from_coords_par<G>(
        coords: &[G::Point1DType],
        grid: &G,
    ) -> Result<Self, ErrorsCoordsByIntervalClassification>
    where
        G: FindIntervalIdOfPoint + HasIntervalIdRange + Sync,
        G::Point1DType: Send + Sync,
    {
        Self::try_new_from_coords_par::<G, DenseData>(coords, grid)
    }

    /// Classifies `coords` against `grid` using **sparse** storage and **parallel** execution.
    ///
    /// Phase 1 maps each coordinate to its interval in parallel via Rayon, preserving
    /// input order. Phase 2 fills the sparse [`BTreeMap`] sequentially so that each
    /// bucket remains sorted by coordinate index.
    ///
    /// Prefer [`try_new_sparse_from_coords`](Self::try_new_sparse_from_coords) for small
    /// inputs where parallelism overhead outweighs the benefit.
    ///
    /// # Errors
    /// Returns [`ErrorsCoordsByIntervalClassification::EmptyCoords`] when `coords` is empty.
    pub fn try_new_sparse_from_coords_par<G>(
        coords: &[G::Point1DType],
        grid: &G,
    ) -> Result<Self, ErrorsCoordsByIntervalClassification>
    where
        G: FindIntervalIdOfPoint + HasIntervalIdRange + Sync,
        G::Point1DType: Send + Sync,
    {
        Self::try_new_from_coords_par::<G, SparseData>(coords, grid)
    }

    /// Shared implementation for [`new_dense_from_coords_sorted`](Self::new_dense_from_coords_sorted)
    /// and [`new_sparse_from_coords_sorted`](Self::new_sparse_from_coords_sorted).
    ///
    /// Performs a single linear sweep through the sorted `coords` slice, advancing a
    /// single `interval_id` cursor forward only when the current coordinate no longer
    /// belongs to the current interval. The storage strategy (dense vs. sparse) is
    /// chosen at compile time via the `DataType` type parameter.
    fn new_from_coords_sorted<G, DataType>(
        coords: &Coords1D<<<G as HasIntervals>::IntervalType as IntervalBoundsRuntime>::RealType>,
        grid: &G,
    ) -> Self
    where
        G: FindIntervalIdOfPoint<Point1DType = <<G as HasIntervals>::IntervalType as IntervalBoundsRuntime>::RealType>
            + HasIntervals,
        DataType: BucketsTrait,
    {
        let num_coords = coords.num_points();
        let coords = coords.as_ref();

        let first_interval_id = grid.first_interval_id();
        let last_interval_id = grid.last_interval_id();

        let num_intervals = grid.num_intervals();

        let mut data = DataType::create_empty(num_intervals, num_coords);

        // look for the first coordinate that is inside the domain, to determine where to start iterating.
        // This is needed to ensure that we iterate in coord_id order, so that each bucket is sorted by construction.
        let mut coord_id = 0; // dummy initialization to satisfy the borrow checker - will be overwritten when first interval is found
        let mut interval_id = first_interval_id; // dummy initialization to satisfy the borrow checker - will be overwritten when first interval is found
        for coord in coords.iter() {
            if let Some(id) = grid.find_interval_id_of_point(coord) {
                interval_id = id;
                break;
            }
            data.push_outside(coord_id); // this coord is outside, but we need to keep track of its coord_id for the final result
            coord_id += 1; // keep looking for the first inside coordinate
        }

        let num_coords = num_coords.into_inner();
        'outer: while coord_id < num_coords {
            // At this point, we have found the first interval that contains at least one coordinate (if any), and we know its interval_id.
            // SAFETY: We ensure that we iterate in coord_id order, so each bucket is sorted by construction.
            let coord = &coords[coord_id]; // safe because we check coord_id < num_coords in the loop condition

            // We know that coords[0..coord_id] are all outside the domain, and coords[coord_id..] are all inside the domain of interval_id or a subsequent interval (if any).
            // We need to check if coords[coord_id] is inside the current interval. If it is, we classify it as inside and move on to the next coord.
            // If it's not, we need to check the next interval(s) until we find one that contains this coord (in which case it's inside the domain) or we exhaust all intervals (in which case it's outside the domain).
            while !grid.interval(&interval_id).contains_point(coord) {
                // This coord is outside the current interval, but we don't know yet if it's outside the whole domain or just this interval.
                // We need to check the next interval(s) until we find one that contains this coord (in which case it's inside the domain) or we exhaust all intervals (in which case it's outside the domain).
                interval_id = IntervalId::new(interval_id.as_ref() + 1);
                if interval_id > last_interval_id {
                    // We've exhausted all intervals without finding one that contains this coord, so it must be outside the domain.
                    //outside_domain.push(coord_id);
                    break 'outer; // we can stop iterating because all subsequent coords will also be outside the domain (since coords are sorted).
                }
            }
            // we found an interval that contains this coord, so we can classify it as inside and move on to the next coord.

            data.push_inside(interval_id, coord_id);
            coord_id += 1;
        }

        for remaining_coord_id in coord_id..num_coords {
            // All remaining coords are outside the domain, since coords are sorted and we've exhausted all intervals.
            data.push_outside(remaining_coord_id);
        }

        data.finalize()
    }

    /// Classifies a **pre-sorted** coordinate slice against `grid` using **dense** storage.
    ///
    /// This is a specialised, infallible alternative to
    /// [`try_new_dense_from_coords`](Self::try_new_dense_from_coords) that exploits the
    /// sorted order of `coords` to sweep through intervals linearly rather than performing
    /// an independent point-location query for every coordinate:
    ///
    /// 1. An initial pass skips coordinates that lie before the first domain interval
    ///    and records them as outside.
    /// 2. A main loop advances a single `interval_id` cursor forward only when needed,
    ///    filling each dense bucket in order — guaranteeing that buckets are sorted by
    ///    construction without any extra sorting step.
    /// 3. A cleanup pass records all remaining coordinates (past the last interval) as
    ///    outside.
    ///
    /// Because `coords` is guaranteed non-empty by the [`Coords1D`] wrapper, this
    /// method is infallible and returns `Self` directly.
    ///
    /// # Panics
    /// Does not panic under normal use. The method relies on `coords` being truly
    /// sorted; passing an unsorted slice via an unsafe construction of [`Coords1D`]
    /// yields unspecified (but memory-safe) results.
    #[must_use]
    pub fn new_dense_from_coords_sorted<G>(
        coords: &Coords1D<<<G as HasIntervals>::IntervalType as IntervalBoundsRuntime>::RealType>,
        grid: &G,
    ) -> Self
    where
        G: FindIntervalIdOfPoint<Point1DType = <<G as HasIntervals>::IntervalType as IntervalBoundsRuntime>::RealType>
            + HasIntervals,
    {
        Self::new_from_coords_sorted::<G, DenseData>(coords, grid)
    }

    /// Classifies a **pre-sorted** coordinate slice against `grid` using **sparse** storage.
    ///
    /// This is the sparse counterpart to
    /// [`new_dense_from_coords_sorted`](Self::new_dense_from_coords_sorted): it uses the
    /// same single-pass linear sweep that exploits sorted order, but stores the inside
    /// mapping in a [`BTreeMap`] so that only intervals that contain at least one coordinate
    /// appear as keys. Prefer this variant when the grid has many intervals but coordinates
    /// are clustered in a small subset of them.
    ///
    /// Because `coords` is guaranteed non-empty by the [`Coords1D`] wrapper, this
    /// method is infallible and returns `Self` directly.
    ///
    /// # Panics
    /// Does not panic under normal use. The method relies on `coords` being truly
    /// sorted; passing an unsorted slice via an unsafe construction of [`Coords1D`]
    /// yields unspecified (but memory-safe) results.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use grid1d::coords::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 1.0),
    ///     NumIntervals::try_new(4).unwrap(),
    /// );
    ///
    /// // Pre-sorted coordinates — some outside the domain
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.3, 0.6, 1.5])).unwrap();
    ///
    /// let classification =
    ///     CoordsByIntervalClassification::new_sparse_from_coords_sorted(&coords, &grid);
    ///
    /// // Only the three intervals that actually contain a point appear
    /// assert_eq!(classification.point_ids_inside().len(), 3);
    /// assert_eq!(classification.point_ids_outside().as_slice(), &[3]); // 1.5
    /// ```
    #[must_use]
    pub fn new_sparse_from_coords_sorted<G>(
        coords: &Coords1D<<<G as HasIntervals>::IntervalType as IntervalBoundsRuntime>::RealType>,
        grid: &G,
    ) -> Self
    where
        G: FindIntervalIdOfPoint<Point1DType = <<G as HasIntervals>::IntervalType as IntervalBoundsRuntime>::RealType>
            + HasIntervals,
    {
        Self::new_from_coords_sorted::<G, SparseData>(coords, grid)
    }
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        Grid1D, Grid1DUniform,
        intervals::{IntervalClosed, bounded::IntervalFromBounds},
        scalars::NumIntervals,
    };
    use try_create::TryNew;

    fn make_grid_3_intervals() -> Grid1DUniform<IntervalClosed<f64>> {
        // [0, 3] with 3 intervals of width 1: [0,1), [1,2), [2,3]
        Grid1DUniform::new(
            IntervalClosed::new(0.0, 3.0),
            NumIntervals::try_new(3).unwrap(),
        )
    }

    mod try_new_dense_from_coords {
        use super::*;
        #[test]
        fn test_all_inside() {
            let grid = make_grid_3_intervals();
            // One point per interval
            let coords = vec![0.5, 1.5, 2.5];
            let result =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Dense { buckets } = &result.inside else {
                panic!("Expected Dense variant");
            };
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![0usize]);
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![1usize]);
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![2usize]);
        }

        #[test]
        fn test_all_outside() {
            let grid = make_grid_3_intervals();
            let coords = vec![-1.0, 5.0];
            let result =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize]);
            let CoordIdsByInterval::Dense { buckets } = &result.inside else {
                panic!("Expected Dense variant");
            };
            assert!(buckets.is_empty());
        }

        #[test]
        fn test_mixed_inside_and_outside() {
            let grid = make_grid_3_intervals();
            // coord 0: outside (below), coord 1: interval 0, coord 2: outside (above), coord 3: interval 2
            let coords = vec![-1.0, 0.5, 10.0, 2.5];
            let result =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 2usize]);
            let CoordIdsByInterval::Dense { buckets } = &result.inside else {
                panic!("Expected Dense variant");
            };
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![1usize]);
            assert!(buckets[1].is_none()); // no points in interval 1
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![3usize]);
        }

        #[test]
        fn test_multiple_points_in_same_interval() {
            let grid = make_grid_3_intervals();
            // coords 0, 1, 2 all fall in interval 1
            let coords = vec![1.1, 1.5, 1.9];
            let result =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Dense { buckets } = &result.inside else {
                panic!("Expected Dense variant");
            };
            assert_eq!(buckets.len(), 2); // only bucket 0 (empty) and bucket 1 are created
            assert!(buckets[0].is_none()); // slot 0 is None: no points in interval 0
            assert_eq!(
                buckets[1].as_ref().unwrap().to_vec(),
                vec![0usize, 1usize, 2usize]
            );
        }

        #[test]
        fn test_boundary_points_closed_domain() {
            let grid = make_grid_3_intervals();
            // For [a,b] domain, interior boundary point pₖ belongs to the RIGHT interval
            // p=1.0 → interval 1, p=2.0 → interval 2
            // Domain boundaries: p=0.0 → interval 0, p=3.0 → interval 2 (last)
            let coords = vec![0.0, 1.0, 2.0, 3.0];
            let result =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Dense { buckets } = &result.inside else {
                panic!("Expected Dense variant");
            };
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![0usize]); // p=0.0 → interval 0
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![1usize]); // p=1.0 → interval 1 (right)
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![2usize, 3usize]); // p=2.0 → interval 2 (right), p=3.0 → interval 2 (last)
        }

        #[test]
        fn test_empty_coords() {
            let grid = make_grid_3_intervals();
            let coords: Vec<f64> = vec![];
            let err = CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                .unwrap_err();
            assert!(matches!(
                err,
                ErrorsCoordsByIntervalClassification::EmptyCoords { .. }
            ));
        }

        #[test]
        fn test_hint_num_intervals_does_not_change_result() {
            let grid = make_grid_3_intervals();
            let coords = vec![0.5, 1.5, 2.5];
            let without_hint =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();
            let with_hint =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();
            assert_eq!(without_hint, with_hint);
        }
    }

    mod try_new_sparse_from_coords {
        use super::*;
        use crate::Grid1DWindow;

        #[test]
        fn test_all_inside() {
            let grid = make_grid_3_intervals();
            // One point per interval
            let coords = vec![0.5, 1.5, 2.5];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid).unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[&IntervalId::new(0)].to_vec(), vec![0usize]);
            assert_eq!(buckets[&IntervalId::new(1)].to_vec(), vec![1usize]);
            assert_eq!(buckets[&IntervalId::new(2)].to_vec(), vec![2usize]);
        }

        #[test]
        fn test_all_outside() {
            let grid = make_grid_3_intervals();
            let coords = vec![-1.0, 5.0];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid).unwrap();

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize]);
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            // No keys in the map when all points are outside
            assert!(buckets.is_empty());
        }

        #[test]
        fn test_mixed_inside_and_outside() {
            let grid = make_grid_3_intervals();
            // coord 0: outside (below), coord 1: interval 0, coord 2: outside (above), coord 3: interval 2
            let coords = vec![-1.0, 0.5, 10.0, 2.5];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid).unwrap();

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 2usize]);
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            // Only intervals that actually contain points appear as keys
            assert_eq!(buckets.len(), 2);
            assert_eq!(buckets[&IntervalId::new(0)].to_vec(), vec![1usize]);
            assert!(!buckets.contains_key(&IntervalId::new(1)));
            assert_eq!(buckets[&IntervalId::new(2)].to_vec(), vec![3usize]);
        }

        #[test]
        fn test_multiple_points_in_same_interval() {
            let grid = make_grid_3_intervals();
            // coords 0, 1, 2 all fall in interval 1
            let coords = vec![1.1, 1.5, 1.9];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid).unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            // Only interval 1 appears as a key - no empty entries for intervals 0 and 2
            assert_eq!(buckets.len(), 1);
            assert_eq!(
                buckets[&IntervalId::new(1)].to_vec(),
                vec![0usize, 1usize, 2usize]
            );
        }

        #[test]
        fn test_boundary_points_closed_domain() {
            let grid = make_grid_3_intervals();
            // For [a,b] domain, interior boundary point pₖ belongs to the RIGHT interval
            let coords = vec![0.0, 1.0, 2.0, 3.0];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid).unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[&IntervalId::new(0)].to_vec(), vec![0usize]); // p=0.0 → interval 0
            assert_eq!(buckets[&IntervalId::new(1)].to_vec(), vec![1usize]); // p=1.0 → interval 1 (right)
            assert_eq!(buckets[&IntervalId::new(2)].to_vec(), vec![2usize, 3usize]); // p=2.0, p=3.0 → interval 2
        }

        #[test]
        fn test_empty_coords_returns_error() {
            let grid = make_grid_3_intervals();
            let coords: Vec<f64> = vec![];
            let err = CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                .unwrap_err();
            assert!(matches!(
                err,
                ErrorsCoordsByIntervalClassification::EmptyCoords { .. }
            ));
        }

        #[test]
        fn test_sparse_on_interval_partition_window() {
            // Grid1DWindow covers intervals 1..=2 of a 4-interval grid.
            // The interval IDs in the result should be 1 and 2, not 0-based.
            let grid = Grid1D::uniform(
                IntervalClosed::new(0.0, 4.0),
                NumIntervals::try_new(4).unwrap(),
            );
            let window =
                Grid1DWindow::try_new(&grid, IntervalId::new(1), NumIntervals::try_new(2).unwrap())
                    .unwrap();

            // Points in interval 1 ([1,2)) and interval 2 ([2,3))
            let coords = vec![1.5, 2.5];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &window)
                    .unwrap();

            assert!(result.outside_domain.is_empty());
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            // Keys must be 1 and 2 (original grid interval IDs), not 0 and 1
            assert_eq!(buckets.len(), 2);
            assert!(buckets.contains_key(&IntervalId::new(1)));
            assert!(buckets.contains_key(&IntervalId::new(2)));
            assert!(!buckets.contains_key(&IntervalId::new(0)));
        }

        #[test]
        fn test_sparse_window_outside_points() {
            // Points outside the window but inside the full grid go to outside_domain
            let grid = Grid1D::uniform(
                IntervalClosed::new(0.0, 4.0),
                NumIntervals::try_new(4).unwrap(),
            );
            let window =
                Grid1DWindow::try_new(&grid, IntervalId::new(1), NumIntervals::try_new(2).unwrap())
                    .unwrap();

            // 0.5 is in interval 0 (outside the window), 1.5 is in interval 1 (inside)
            let coords = vec![0.5, 1.5];
            let result =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &window)
                    .unwrap();

            assert_eq!(result.outside_domain.to_vec(), vec![0usize]);
            let CoordIdsByInterval::Sparse { buckets } = &result.inside else {
                panic!("Expected Sparse variant");
            };
            assert_eq!(buckets.len(), 1);
            assert_eq!(buckets[&IntervalId::new(1)].to_vec(), vec![1usize]);
        }
    }

    mod try_new_dense_from_coords_par {
        use super::*;

        #[test]
        fn test_matches_sequential() {
            let grid = make_grid_3_intervals();
            let coords = vec![-1.0, 0.5, 1.5, 10.0, 2.5];
            let seq =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();
            let par = CoordsByIntervalClassification::try_new_dense_from_coords_par(&coords, &grid)
                .unwrap();
            assert_eq!(seq, par);
        }

        #[test]
        fn test_hint_matches_sequential() {
            let grid = make_grid_3_intervals();
            let coords = vec![0.5, 1.5, 2.5];
            let seq =
                CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid).unwrap();
            let par = CoordsByIntervalClassification::try_new_dense_from_coords_par(&coords, &grid)
                .unwrap();
            assert_eq!(seq, par);
        }

        #[test]
        fn test_empty_coords() {
            let grid = make_grid_3_intervals();
            let coords: Vec<f64> = vec![];
            let err = CoordsByIntervalClassification::try_new_dense_from_coords_par(&coords, &grid)
                .unwrap_err();
            assert!(matches!(
                err,
                ErrorsCoordsByIntervalClassification::EmptyCoords { .. }
            ));
        }
    }

    mod try_new_sparse_from_coords_par {
        use super::*;

        #[test]
        fn test_matches_sequential() {
            let grid = make_grid_3_intervals();
            let coords = vec![-1.0, 0.5, 1.5, 10.0, 2.5];
            let seq =
                CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid).unwrap();
            let par =
                CoordsByIntervalClassification::try_new_sparse_from_coords_par(&coords, &grid)
                    .unwrap();
            assert_eq!(seq, par);
        }

        #[test]
        fn test_empty_coords() {
            let grid = make_grid_3_intervals();
            let coords: Vec<f64> = vec![];
            let err =
                CoordsByIntervalClassification::try_new_sparse_from_coords_par(&coords, &grid)
                    .unwrap_err();
            assert!(matches!(
                err,
                ErrorsCoordsByIntervalClassification::EmptyCoords { .. }
            ));
        }
    }

    mod new_dense_from_coords_sorted {
        use super::*;
        use sorted_vec::partial::SortedSet;

        fn make_coords(values: Vec<f64>) -> Coords1D<f64> {
            Coords1D::try_from(SortedSet::from_unsorted(values)).unwrap()
        }

        // ── helpers ──────────────────────────────────────────────────────────

        fn call(coords: &Coords1D<f64>) -> CoordsByIntervalClassification {
            let grid = make_grid_3_intervals();
            CoordsByIntervalClassification::new_dense_from_coords_sorted(coords, &grid)
        }

        fn dense_buckets(
            result: &CoordsByIntervalClassification,
        ) -> &Vec<Option<sorted_vec::SortedSet<usize>>> {
            let CoordIdsByInterval::Dense { buckets } = &result.inside else {
                panic!("Expected Dense variant");
            };
            buckets
        }

        // ── basic correctness ────────────────────────────────────────────────

        #[test]
        fn test_all_inside() {
            // One point per interval, no outside coords.
            let coords = make_coords(vec![0.5, 1.5, 2.5]);
            let result = call(&coords);

            assert!(result.outside_domain.is_empty());
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![0usize]);
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![1usize]);
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![2usize]);
        }

        #[test]
        fn test_all_outside_left() {
            // All coords are below the domain: the initial for loop exhausts all
            // coords without finding any inside, so coord_id == num_coords and
            // the outer while must NOT be entered (this was the `||` bug scenario).
            let coords = make_coords(vec![-2.0, -1.0]);
            let result = call(&coords);

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize]);
            let buckets = dense_buckets(&result);
            assert!(buckets.is_empty());
        }

        #[test]
        fn test_all_outside_right() {
            // All coords are above the domain: the inner while exhausts all
            // intervals and triggers `break 'outer`, then the cleanup loop
            // pushes the remaining coord_ids.
            let coords = make_coords(vec![5.0, 6.0, 7.0]);
            let result = call(&coords);

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize, 2usize]);
            let buckets = dense_buckets(&result);
            assert!(buckets.is_empty());
        }

        #[test]
        fn test_outside_left_then_inside() {
            // Coords before the domain are placed in outside_domain by the
            // initial for loop; the remaining coords are classified normally.
            let coords = make_coords(vec![-1.0, 0.5, 1.5, 2.5]);
            let result = call(&coords);

            assert_eq!(result.outside_domain.to_vec(), vec![0usize]);
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![1usize]);
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![2usize]);
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![3usize]);
        }

        #[test]
        fn test_inside_then_outside_right() {
            // Coords after the domain end trigger `break 'outer`; the cleanup
            // loop must add those coord_ids to outside_domain.
            let coords = make_coords(vec![0.5, 1.5, 2.5, 5.0, 6.0]);
            let result = call(&coords);

            assert_eq!(result.outside_domain.to_vec(), vec![3usize, 4usize]);
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![0usize]);
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![1usize]);
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![2usize]);
        }

        #[test]
        fn test_all_three_regions() {
            // outside-left + inside + outside-right: exercises all code paths.
            // coord 0 (-2.0): outside left (initial for loop)
            // coord 1 (0.5):  interval 0
            // coord 2 (1.5):  interval 1
            // coord 3 (2.5):  interval 2
            // coord 4 (9.0):  outside right (cleanup loop)
            let coords = make_coords(vec![-2.0, 0.5, 1.5, 2.5, 9.0]);
            let result = call(&coords);

            assert_eq!(result.outside_domain.to_vec(), vec![0usize, 4usize]);
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![1usize]);
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![2usize]);
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![3usize]);
        }

        #[test]
        fn test_multiple_points_in_same_interval() {
            // Three coords in interval 1; intervals 0 and 2 get no points.
            let coords = make_coords(vec![1.1, 1.5, 1.9]);
            let result = call(&coords);

            assert!(result.outside_domain.is_empty());
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 2); // slots 0 and 1; slot 0 is None
            assert!(buckets[0].is_none());
            assert_eq!(
                buckets[1].as_ref().unwrap().to_vec(),
                vec![0usize, 1usize, 2usize]
            );
        }

        #[test]
        fn test_boundary_points_closed_domain() {
            // p=0.0 → interval 0 (left domain boundary)
            // p=1.0 → interval 1 (interior boundary, assigned to the right interval)
            // p=2.0 → interval 2 (interior boundary, assigned to the right interval)
            // p=3.0 → interval 2 (right domain boundary, last interval)
            let coords = make_coords(vec![0.0, 1.0, 2.0, 3.0]);
            let result = call(&coords);

            assert!(result.outside_domain.is_empty());
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 3);
            assert_eq!(buckets[0].as_ref().unwrap().to_vec(), vec![0usize]);
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![1usize]);
            assert_eq!(buckets[2].as_ref().unwrap().to_vec(), vec![2usize, 3usize]);
        }

        #[test]
        fn test_single_coord_inside() {
            let coords = make_coords(vec![1.5]);
            let result = call(&coords);

            assert!(result.outside_domain.is_empty());
            let buckets = dense_buckets(&result);
            assert_eq!(buckets.len(), 2); // slot 0 is None, slot 1 has the point
            assert!(buckets[0].is_none());
            assert_eq!(buckets[1].as_ref().unwrap().to_vec(), vec![0usize]);
        }

        #[test]
        fn test_single_coord_outside_left() {
            let coords = make_coords(vec![-1.0]);
            let result = call(&coords);

            assert_eq!(result.outside_domain.to_vec(), vec![0usize]);
            assert!(dense_buckets(&result).is_empty());
        }

        #[test]
        fn test_matches_try_new_dense_from_coords() {
            // For sorted inputs the two methods must produce identical results.
            let grid = make_grid_3_intervals();
            let values = vec![-1.0, 0.5, 1.5, 2.5, 9.0];
            let coords_obj = make_coords(values.clone());

            let via_sorted =
                CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords_obj, &grid);
            let via_generic =
                CoordsByIntervalClassification::try_new_dense_from_coords(&values, &grid).unwrap();

            assert_eq!(via_sorted, via_generic);
        }

        #[test]
        fn test_all_outside_right_matches_try_new() {
            // Regression for the `||` bug: previously panicked with OOB access
            // when all coords are outside the right boundary.
            let grid = make_grid_3_intervals();
            let values = vec![5.0, 6.0, 7.0];
            let coords_obj = make_coords(values.clone());

            let via_sorted =
                CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords_obj, &grid);
            let via_generic =
                CoordsByIntervalClassification::try_new_dense_from_coords(&values, &grid).unwrap();

            assert_eq!(via_sorted, via_generic);
        }

        #[test]
        fn test_all_outside_left_matches_try_new() {
            // Regression for the `||` bug: previously panicked with OOB access
            // when all coords are outside the left boundary.
            let grid = make_grid_3_intervals();
            let values = vec![-2.0, -1.0];
            let coords_obj = make_coords(values.clone());

            let via_sorted =
                CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords_obj, &grid);
            let via_generic =
                CoordsByIntervalClassification::try_new_dense_from_coords(&values, &grid).unwrap();

            assert_eq!(via_sorted, via_generic);
        }

        // ── window tests ─────────────────────────────────────────────────────

        mod window {
            use super::*;
            use crate::{Grid1D, Grid1DNonUniform, Grid1DUniform, Grid1DWindow};

            // 5-interval uniform grid [0.0, 5.0]:
            //   0:[0.0,1.0)  1:[1.0,2.0)  2:[2.0,3.0)  3:[3.0,4.0)  4:[4.0,5.0]
            fn make_grid_5() -> Grid1DUniform<IntervalClosed<f64>> {
                Grid1DUniform::new(
                    IntervalClosed::new(0.0, 5.0),
                    NumIntervals::try_new(5).unwrap(),
                )
            }

            // Interior window: absolute intervals 1..=3, covers [1.0, 4.0)
            fn interior_window(
                grid: &Grid1DUniform<IntervalClosed<f64>>,
            ) -> Grid1DWindow<'_, Grid1DUniform<IntervalClosed<f64>>> {
                Grid1DWindow::try_new(grid, IntervalId::new(1), NumIntervals::try_new(3).unwrap())
                    .unwrap()
            }

            // End window: absolute intervals 3..=4, covers [3.0, 5.0] (right boundary closed)
            fn end_window(
                grid: &Grid1DUniform<IntervalClosed<f64>>,
            ) -> Grid1DWindow<'_, Grid1DUniform<IntervalClosed<f64>>> {
                Grid1DWindow::try_new(grid, IntervalId::new(3), NumIntervals::try_new(2).unwrap())
                    .unwrap()
            }

            // ── interior-window tests ─────────────────────────────────────────

            #[test]
            fn test_all_inside_interior_window() {
                // One point per window interval.
                // Absolute interval IDs 1, 2, 3 → buckets[1], [2], [3]; buckets[0] is None.
                let grid = make_grid_5();
                let w = interior_window(&grid);
                let coords = super::make_coords(vec![1.5, 2.5, 3.5]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert!(result.outside_domain.is_empty());
                let b = super::dense_buckets(&result);
                assert_eq!(b.len(), 4); // indices 0..=3
                assert!(b[0].is_none());
                assert_eq!(b[1].as_ref().unwrap().to_vec(), vec![0usize]);
                assert_eq!(b[2].as_ref().unwrap().to_vec(), vec![1usize]);
                assert_eq!(b[3].as_ref().unwrap().to_vec(), vec![2usize]);
            }

            #[test]
            fn test_before_window_is_outside() {
                // Coords in [0.0, 1.0): inside parent grid but before the window start.
                // The initial for-loop exhausts all coords without finding any inside.
                let grid = make_grid_5();
                let w = interior_window(&grid);
                let coords = super::make_coords(vec![0.2, 0.7]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize]);
                assert!(super::dense_buckets(&result).is_empty());
            }

            #[test]
            fn test_after_window_is_outside() {
                // Coords in (4.0, 5.0]: inside parent grid but past the open window end.
                // The 'outer cleanup loop pushes them to outside_domain.
                let grid = make_grid_5();
                let w = interior_window(&grid);
                let coords = super::make_coords(vec![4.1, 4.9]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize]);
                assert!(super::dense_buckets(&result).is_empty());
            }

            #[test]
            fn test_all_three_regions() {
                // coord 0 (0.5) → outside (before window, initial loop)
                // coord 1 (1.5) → interval 1
                // coord 2 (2.5) → interval 2
                // coord 3 (3.5) → interval 3
                // coord 4 (4.5) → outside (after window, cleanup loop)
                let grid = make_grid_5();
                let w = interior_window(&grid);
                let coords = super::make_coords(vec![0.5, 1.5, 2.5, 3.5, 4.5]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert_eq!(result.outside_domain.to_vec(), vec![0usize, 4usize]);
                let b = super::dense_buckets(&result);
                assert_eq!(b.len(), 4);
                assert!(b[0].is_none());
                assert_eq!(b[1].as_ref().unwrap().to_vec(), vec![1usize]);
                assert_eq!(b[2].as_ref().unwrap().to_vec(), vec![2usize]);
                assert_eq!(b[3].as_ref().unwrap().to_vec(), vec![3usize]);
            }

            #[test]
            fn test_window_boundary_points() {
                // x = 1.0: left edge of window → interval 1 (left-closed [1.0, 2.0))
                // x = 2.0: interior boundary → interval 2 (left-closed, so right interval)
                // x = 3.0: interior boundary → interval 3
                // x = 4.0: open right edge of interior window → outside
                let grid = make_grid_5();
                let w = interior_window(&grid);
                let coords = super::make_coords(vec![1.0, 2.0, 3.0, 4.0]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert_eq!(result.outside_domain.to_vec(), vec![3usize]); // x=4.0
                let b = super::dense_buckets(&result);
                assert_eq!(b.len(), 4);
                assert!(b[0].is_none());
                assert_eq!(b[1].as_ref().unwrap().to_vec(), vec![0usize]); // x=1.0
                assert_eq!(b[2].as_ref().unwrap().to_vec(), vec![1usize]); // x=2.0
                assert_eq!(b[3].as_ref().unwrap().to_vec(), vec![2usize]); // x=3.0
            }

            // ── end-window test ───────────────────────────────────────────────

            #[test]
            fn test_end_window_right_boundary_is_closed() {
                // End window covers [3.0, 5.0]: intervals 3 and 4.
                //   interval 3: [3.0, 4.0), interval 4: [4.0, 5.0]
                // Since the window ends at the last partition interval, 5.0 is inside.
                let grid = make_grid_5();
                let w = end_window(&grid);
                let coords = super::make_coords(vec![3.0, 3.5, 4.0, 5.0]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert!(result.outside_domain.is_empty());
                let b = super::dense_buckets(&result);
                assert_eq!(b.len(), 5); // indices 0..=4; first 3 are None
                assert!(b[0].is_none());
                assert!(b[1].is_none());
                assert!(b[2].is_none());
                assert_eq!(b[3].as_ref().unwrap().to_vec(), vec![0usize, 1usize]); // 3.0, 3.5
                assert_eq!(b[4].as_ref().unwrap().to_vec(), vec![2usize, 3usize]); // 4.0, 5.0
            }

            // ── cross-check test ──────────────────────────────────────────────

            #[test]
            fn test_matches_try_new_dense_from_coords() {
                // For sorted input, new_dense_from_coords_sorted and
                // try_new_dense_from_coords must produce identical results.
                let grid = make_grid_5();
                let w = interior_window(&grid);
                let values = vec![0.5, 1.5, 2.5, 3.5, 4.5];
                let coords_obj = super::make_coords(values.clone());

                let via_sorted =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords_obj, &w);
                let via_generic =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&values, &w).unwrap();

                assert_eq!(via_sorted, via_generic);
            }

            // ── non-uniform grid window test ──────────────────────────────────

            #[test]
            fn test_non_uniform_interior_window() {
                // Non-uniform grid from coords [0.0, 1.0, 3.0, 5.0, 6.0]:
                //   interval 0: [0.0, 1.0), 1: [1.0, 3.0), 2: [3.0, 5.0), 3: [5.0, 6.0]
                // Interior window: intervals 1..=2, covers [1.0, 5.0) (right-open)
                //
                // coord 0 (0.5) → outside (before window)
                // coord 1 (2.0) → interval 1 ([1.0, 3.0))
                // coord 2 (3.0) → interval 2 ([3.0, 5.0), interior boundary → right interval)
                // coord 3 (4.0) → interval 2
                // coord 4 (5.0) → outside (open right edge of interior window)
                // coord 5 (6.0) → outside
                let grid =
                    Grid1D::<IntervalClosed<f64>>::try_from_sorted(SortedSet::from_unsorted(vec![
                        0.0, 1.0, 3.0, 5.0, 6.0,
                    ]))
                    .unwrap();
                let w = Grid1DWindow::try_new(
                    &grid,
                    IntervalId::new(1),
                    NumIntervals::try_new(2).unwrap(),
                )
                .unwrap();
                let coords = super::make_coords(vec![0.5, 2.0, 3.0, 4.0, 5.0, 6.0]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert_eq!(result.outside_domain.to_vec(), vec![0usize, 4usize, 5usize]);
                let b = super::dense_buckets(&result);
                assert_eq!(b.len(), 3); // indices 0..=2; 0 is None
                assert!(b[0].is_none());
                assert_eq!(b[1].as_ref().unwrap().to_vec(), vec![1usize]); // x=2.0
                assert_eq!(b[2].as_ref().unwrap().to_vec(), vec![2usize, 3usize]); // x=3.0, x=4.0
            }

            #[test]
            fn test_non_uniform_direct_grid1d_non_uniform() {
                // Same grid as test_non_uniform_interior_window but using Grid1DNonUniform
                // directly, which has its own FindIntervalIdOfPoint impl (binary-search based).
                // Grid coords [0.0, 1.0, 3.0, 5.0, 6.0]:
                //   interval 0: [0.0, 1.0), 1: [1.0, 3.0), 2: [3.0, 5.0), 3: [5.0, 6.0]
                // End window: intervals 2..=3, covers [3.0, 6.0] (right boundary closed)
                //
                // coord 0 (0.5) → outside (before window)
                // coord 1 (1.5) → outside (before window)
                // coord 2 (3.5) → interval 2
                // coord 3 (4.5) → interval 2
                // coord 4 (5.0) → interval 3 ([5.0, 6.0], interior boundary → right interval)
                // coord 5 (6.0) → interval 3 (closed right end)
                let grid_coords = super::make_coords(vec![0.0, 1.0, 3.0, 5.0, 6.0]);
                let grid =
                    Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(grid_coords)
                        .unwrap();
                let w = Grid1DWindow::try_new(
                    &grid,
                    IntervalId::new(2),
                    NumIntervals::try_new(2).unwrap(),
                )
                .unwrap();
                let coords = super::make_coords(vec![0.5, 1.5, 3.5, 4.5, 5.0, 6.0]);
                let result =
                    CoordsByIntervalClassification::new_dense_from_coords_sorted(&coords, &w);

                assert_eq!(result.outside_domain.to_vec(), vec![0usize, 1usize]);
                let b = super::dense_buckets(&result);
                assert_eq!(b.len(), 4); // indices 0..=3; 0 and 1 are None
                assert!(b[0].is_none());
                assert!(b[1].is_none());
                assert_eq!(b[2].as_ref().unwrap().to_vec(), vec![2usize, 3usize]); // 3.5, 4.5
                assert_eq!(b[3].as_ref().unwrap().to_vec(), vec![4usize, 5usize]); // 5.0, 6.0
            }
        }
    }

    mod accessors {
        use super::*;

        // ── point_ids_outside ─────────────────────────────────────────────────

        mod point_ids_outside {
            use super::*;

            #[test]
            fn dense_all_inside() {
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 1.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert!(result.point_ids_outside().is_empty());
            }

            #[test]
            fn dense_all_outside() {
                let grid = make_grid_3_intervals();
                let coords = vec![-1.0, 5.0];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(result.point_ids_outside().to_vec(), vec![0usize, 1usize]);
            }

            #[test]
            fn dense_mixed() {
                let grid = make_grid_3_intervals();
                // coord 0 (-1.0): outside, coord 1 (0.5): inside, coord 2 (10.0): outside, coord 3 (2.5): inside
                let coords = vec![-1.0, 0.5, 10.0, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(result.point_ids_outside().to_vec(), vec![0usize, 2usize]);
            }

            #[test]
            fn sparse_all_inside() {
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 1.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert!(result.point_ids_outside().is_empty());
            }

            #[test]
            fn sparse_all_outside() {
                let grid = make_grid_3_intervals();
                let coords = vec![-1.0, 5.0];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(result.point_ids_outside().to_vec(), vec![0usize, 1usize]);
            }

            #[test]
            fn sparse_mixed() {
                let grid = make_grid_3_intervals();
                let coords = vec![-1.0, 0.5, 10.0, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(result.point_ids_outside().to_vec(), vec![0usize, 2usize]);
            }
        }

        // ── point_ids_inside ──────────────────────────────────────────────────

        mod point_ids_inside {
            use super::*;

            // Helper: collapses Vec<(IntervalId, &SortedSet<usize>)> into
            // Vec<(usize, Vec<usize>)> for easy assertions.
            fn collect_pairs(result: &CoordsByIntervalClassification) -> Vec<(usize, Vec<usize>)> {
                result
                    .point_ids_inside()
                    .into_iter()
                    .map(|(id, set)| (*id.as_ref(), set.to_vec()))
                    .collect()
            }

            #[test]
            fn dense_all_inside() {
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 1.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(0, vec![0usize]), (1, vec![1usize]), (2, vec![2usize])]
                );
            }

            #[test]
            fn dense_all_outside() {
                let grid = make_grid_3_intervals();
                let coords = vec![-1.0, 5.0];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert!(collect_pairs(&result).is_empty());
            }

            #[test]
            fn dense_mixed() {
                let grid = make_grid_3_intervals();
                // coord 0 (-1.0): outside, coord 1 (0.5): interval 0,
                // coord 2 (10.0): outside, coord 3 (2.5): interval 2
                // Interval 1 has no points → must not appear in the result.
                let coords = vec![-1.0, 0.5, 10.0, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(0, vec![1usize]), (2, vec![3usize])]
                );
            }

            #[test]
            fn dense_skips_empty_intervals() {
                // All coords in interval 1: intervals 0 and 2 must not appear.
                let grid = make_grid_3_intervals();
                let coords = vec![1.1, 1.5, 1.9];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(1, vec![0usize, 1usize, 2usize])]
                );
            }

            #[test]
            fn dense_pairs_in_interval_id_ascending_order() {
                // coord 0 (1.5) → interval 1, coord 1 (0.5) → interval 0, coord 2 (2.5) → interval 2.
                // Dense iterates buckets in index order → pairs must be ordered by IntervalId.
                let grid = make_grid_3_intervals();
                let coords = vec![1.5, 0.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(0, vec![1usize]), (1, vec![0usize]), (2, vec![2usize])]
                );
            }

            #[test]
            fn sparse_all_inside() {
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 1.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(0, vec![0usize]), (1, vec![1usize]), (2, vec![2usize])]
                );
            }

            #[test]
            fn sparse_all_outside() {
                let grid = make_grid_3_intervals();
                let coords = vec![-1.0, 5.0];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert!(collect_pairs(&result).is_empty());
            }

            #[test]
            fn sparse_mixed() {
                let grid = make_grid_3_intervals();
                let coords = vec![-1.0, 0.5, 10.0, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(0, vec![1usize]), (2, vec![3usize])]
                );
            }

            #[test]
            fn sparse_pairs_in_interval_id_ascending_order() {
                // coord 0 (1.5) → interval 1, coord 1 (0.5) → interval 0, coord 2 (2.5) → interval 2.
                // BTreeMap iterates by key → pairs must be ordered by IntervalId.
                let grid = make_grid_3_intervals();
                let coords = vec![1.5, 0.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    collect_pairs(&result),
                    vec![(0, vec![1usize]), (1, vec![0usize]), (2, vec![2usize])]
                );
            }
        }

        // ── point_ids_in_interval ─────────────────────────────────────────────

        mod point_ids_in_interval {
            use super::*;

            #[test]
            fn dense_returns_ids_for_occupied_interval() {
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 1.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(0))
                        .unwrap()
                        .to_vec(),
                    vec![0usize]
                );
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(1))
                        .unwrap()
                        .to_vec(),
                    vec![1usize]
                );
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(2))
                        .unwrap()
                        .to_vec(),
                    vec![2usize]
                );
            }

            #[test]
            fn dense_returns_err_for_empty_slot() {
                // All coords in interval 1 → Dense slot 0 is None.
                let grid = make_grid_3_intervals();
                let coords = vec![1.1, 1.5, 1.9];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                let err = result
                    .try_point_ids_in_interval(&IntervalId::new(0))
                    .unwrap_err();
                assert!(
                    matches!(err, ErrorsCoordsByIntervalClassification::IntervalNotPresent { interval_id, .. } if interval_id == IntervalId::new(0))
                );
            }

            #[test]
            fn dense_returns_none_for_out_of_range_id() {
                // Single coord in interval 1 → buckets has len 2 (indices 0 and 1).
                // IntervalId::new(2) is beyond buckets.len() - 1.
                let grid = make_grid_3_intervals();
                let coords = vec![1.5];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                let err = result
                    .try_point_ids_in_interval(&IntervalId::new(2))
                    .unwrap_err();
                assert!(
                    matches!(err, ErrorsCoordsByIntervalClassification::IntervalNotPresent { interval_id, .. } if interval_id == IntervalId::new(2))
                );
            }

            #[test]
            fn dense_multiple_points_same_interval() {
                let grid = make_grid_3_intervals();
                let coords = vec![1.1, 1.5, 1.9];
                let result =
                    CoordsByIntervalClassification::try_new_dense_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(1))
                        .unwrap()
                        .to_vec(),
                    vec![0usize, 1usize, 2usize]
                );
            }

            #[test]
            fn sparse_returns_ids_for_occupied_interval() {
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 1.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(0))
                        .unwrap()
                        .to_vec(),
                    vec![0usize]
                );
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(2))
                        .unwrap()
                        .to_vec(),
                    vec![2usize]
                );
            }

            #[test]
            fn sparse_returns_none_for_absent_interval() {
                // coord 0 (0.5) → interval 0, coord 1 (2.5) → interval 2.
                // Interval 1 has no key in the BTreeMap.
                let grid = make_grid_3_intervals();
                let coords = vec![0.5, 2.5];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                let err = result
                    .try_point_ids_in_interval(&IntervalId::new(1))
                    .unwrap_err();
                assert!(
                    matches!(err, ErrorsCoordsByIntervalClassification::IntervalNotPresent { interval_id, .. } if interval_id == IntervalId::new(1))
                );
            }

            #[test]
            fn sparse_multiple_points_same_interval() {
                let grid = make_grid_3_intervals();
                let coords = vec![1.1, 1.5, 1.9];
                let result =
                    CoordsByIntervalClassification::try_new_sparse_from_coords(&coords, &grid)
                        .unwrap();
                assert_eq!(
                    result
                        .try_point_ids_in_interval(&IntervalId::new(1))
                        .unwrap()
                        .to_vec(),
                    vec![0usize, 1usize, 2usize]
                );
            }
        }
    }
}