grid1d 0.5.1

#![deny(rustdoc::broken_intra_doc_links)]

//! Core traits for grid1d.
//!
//! This module contains all the fundamental trait definitions used throughout the `grid1d` library.
//! These traits define the core abstractions for domains, coordinates, interval partitions,
//! and point transformations.
//!
//! # Trait Overview
//!
//! | Trait | Purpose | Key Methods |
//! |-------|---------|-------------|
//! | [`HasDomain1D`] | Domain access | `domain()` |
//! | [`HasCoords1D`] | Coordinate access | `coords()`, `num_points()` |
//! | [`HasIntervalIdRange`] | Interval index range | `first_interval_id()`, `last_interval_id()`, `num_intervals()` |
//! | [`HasCoordIdRange`] | Coordinate index range | `first_coord_id()`, `last_coord_id()`, `num_coord_ids()` |
//! | [`HasIntervals`] | Interval lookup and iteration | `interval()`, `iter_intervals()` |
//! | [`Grid1DIntervalBuilder`] | Sub-interval construction rules | `build_single_interval()`, `build_first_interval()`, `build_middle_interval()`, `build_last_interval()` |
//! | [`Grid1DTrait`] | Main partition interface | `find_interval_id_of_point()`, `refine()` |
//!
//! # Trait Hierarchy
//!
//! ```text
//! Grid1DTrait
//!     ├── HasCoords1D      (coordinate access: coords(), num_points())
//!     ├── HasIntervals     (interval lookup: interval(), iter_intervals())
//!     │       └── HasIntervalIdRange  (index range: first_interval_id(), last_interval_id(), num_intervals())
//!     ├── FindIntervalIdOfPoint  (point location: find_interval_id_of_point())
//!     └── HasDomain1D      (domain access: domain())
//!             └── Domain1D: Grid1DIntervalBuilder
//!
//! HasIntervals
//!     Implementors: Grid1DUniform, Grid1DNonUniform, Grid1D, Grid1DUnion,
//!                   Grid1DWindow
//!
//! HasIntervalIdRange
//!     Implementors: all HasIntervals implementors + IntervalIndexSpace1D
//!
//! HasCoordIdRange  (coordinate index range, independent of HasIntervalIdRange)
//!     Implementors: Grid1DUniform, Grid1DNonUniform, Grid1D + CoordIndexSpace1D
//!
//! Grid1DIntervalBuilder
//!     └── IntervalFinitePositiveLengthTrait
//!     Methods: build_single_interval(self), build_first_interval(&self, coords),
//!              build_middle_interval(&self, coords, interval_id), build_last_interval(&self, coords)
//! ```
//!
//! # Usage
//!
//! Most users will interact with these traits through the concrete types like
//! [`Grid1D`](crate::Grid1D), [`Grid1DUniform`](crate::Grid1DUniform), and
//! [`Grid1DNonUniform`](crate::Grid1DNonUniform). However, understanding these traits
//! is essential for writing generic code that works with any grid type.
//!
//! ## Example: Generic Function Over Any Partition
//!
//! ```rust
//! use grid1d::{Grid1DTrait, HasIntervalIdRange, HasCoords1D, HasDomain1D, scalars::IntervalId};
//! use grid1d::intervals::IntervalTrait;
//!
//! fn analyze_partition<P>(partition: &P) -> (usize, usize)
//! where
//!     P: Grid1DTrait,
//! {
//!     let num_points = *partition.num_points().as_ref();
//!     let num_intervals = *partition.num_intervals().as_ref();
//!     (num_points, num_intervals)
//! }
//! ```
//!
//! ## Example: Generic Function Over Any `HasIntervalIdRange`
//!
//! ```rust
//! use grid1d::{HasIntervalIdRange, scalars::{IntervalId, NumIntervals}};
//!
//! fn describe_range<R: HasIntervalIdRange>(range: &R) -> (IntervalId, IntervalId, NumIntervals) {
//!     (range.first_interval_id(), range.last_interval_id(), range.num_intervals())
//! }
//! ```
//!
//! ## Example: Generic Function Over Any `HasCoordIdRange`
//!
//! ```rust
//! use grid1d::{HasCoordIdRange, scalars::{CoordId, PositiveNumPoints1D}};
//!
//! fn describe_coord_range<R: HasCoordIdRange>(range: &R) -> (CoordId, CoordId, PositiveNumPoints1D) {
//!     (range.first_coord_id(), range.last_coord_id(), range.num_coord_ids())
//! }
//! ```

use crate::{
    Grid1DIndexSpaces, Side, SubIntervalInPartition,
    coords::Coords1D,
    intervals::{
        GetLowerBoundValue, GetUpperBoundValue, Interval, IntervalClosed, IntervalFiniteLength,
        IntervalFinitePositiveLength, IntervalFinitePositiveLengthTrait,
        IntervalLowerClosedUpperOpen, IntervalLowerOpenUpperClosed, IntervalOpen, IntervalTrait,
        bounded::IntervalFromBounds, operations::IntervalOperations,
    },
    operations::refinement::{Grid1DNonUniformRefinement, Grid1DUniformRefinement},
    scalars::{CoordId, IntervalId, NumIntervals, PositiveNumPoints1D},
    topology::Adjacency1D,
};
use duplicate::duplicate_item;
use more_asserts::debug_assert_lt;
use num_valid::{RealScalar, scalars::PositiveRealScalar};
use rayon::prelude::*;
use serde::{Deserialize, Serialize};
use std::collections::BTreeMap;
use try_create::TryNew;

// ============================================================================
// HasDomain1D
// ============================================================================

/// A trait for types that have an associated one-dimensional domain.
///
/// This trait provides access to the domain (interval) upon which a structure is defined.
/// It is one of the fundamental building blocks for the [`Grid1DTrait`] trait.
///
/// # Associated Types
///
/// - `Domain1D`: The interval type representing the domain. Must implement [`IntervalTrait`].
///
/// # Example
///
/// ```rust
/// use grid1d::{Grid1D, HasDomain1D, intervals::*, scalars::NumIntervals};
/// use grid1d::intervals::{GetLowerBoundValue, GetUpperBoundValue};
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(
///     IntervalClosed::new(0.0, 10.0),
///     NumIntervals::try_new(5).unwrap()
/// );
///
/// let domain = grid.domain();
/// assert_eq!(domain.lower_bound_value(), &0.0);
/// assert_eq!(domain.upper_bound_value(), &10.0);
/// ```
pub trait HasDomain1D: Sized {
    /// The type of interval upon which the current one-dimensional domain is defined.
    type Domain1D: IntervalTrait;

    /// Returns a reference to the one-dimensional domain upon which the current structure is defined.
    fn domain(&self) -> &Self::Domain1D;
}

// ============================================================================
// HasCoords1D
// ============================================================================

/// A trait for types that contain one-dimensional coordinate data.
///
/// This trait provides access to the underlying [`Coords1D`] container, which guarantees
/// that coordinates are unique, sorted in ascending order, and non-empty.
///
/// # Associated Types
///
/// - `Point1DType`: The scalar type for the coordinates. Must implement [`RealScalar`].
///
/// # Provided Methods
///
/// - [`num_points()`](HasCoords1D::num_points): Returns the number of points (delegates to `coords().num_points()`)
///
/// # Example
///
/// ```rust
/// use grid1d::{Grid1D, HasCoords1D, intervals::*, scalars::NumIntervals};
/// use std::ops::Deref;
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(4).unwrap()
/// );
///
/// // Access coordinates via the trait method
/// let coords = grid.coords();
/// assert_eq!(coords.deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
/// assert_eq!(grid.num_points().as_ref(), &5);
/// ```
pub trait HasCoords1D {
    /// The type of the 1D coordinates.
    type CoordType: RealScalar;

    /// Returns a reference to the coordinate corresponding to the given `CoordId`.
    #[inline(always)]
    fn coord(&self, coord_id: CoordId) -> &Self::CoordType {
        &self.coords()[coord_id]
    }

    /// Returns the coordinates of the 1D points.
    ///
    /// This method provides access to the underlying [`Coords1D`] container, which guarantees
    /// that the coordinates are unique, sorted in ascending order, and non-empty.
    ///
    /// # Returns
    ///
    /// A reference to the [`Coords1D<Self::CoordType>`](Coords1D) object.
    fn coords(&self) -> &Coords1D<Self::CoordType>;

    /// Returns the number of points.
    ///
    /// This is a convenience method that delegates to `self.coords().num_points()`.
    /// It provides a quick way to get the total number of discrete points defining the object.
    ///
    /// # Returns
    ///
    /// The number of points as a [`PositiveNumPoints1D`] object, which guarantees the count is at least 1.
    #[inline]
    fn num_points(&self) -> PositiveNumPoints1D {
        self.coords().num_points()
    }
}

// ============================================================================
// Grid1DIntervalBuilder
// ============================================================================

/// A trait for building sub-intervals within a partition based on the domain type.
///
/// The `Grid1DIntervalBuilder` trait defines how to construct the first, middle, and last
/// sub-intervals of a partition, ensuring that the semantics of the domain boundaries
/// are correctly preserved in the resulting partition structure.
///
/// This trait is automatically implemented for all interval types that implement
/// [`IntervalFinitePositiveLengthTrait`], and it works in conjunction with the
/// [`Grid1DTrait`] trait to create mathematically correct partitions.
///
/// # Partition Strategy
///
/// The trait implements a specific strategy for constructing sub-intervals:
///
/// - **First interval**: Constructed using [`build_first_interval`](Grid1DIntervalBuilder::build_first_interval)
/// - **Middle intervals**: Constructed using [`build_middle_interval`](Grid1DIntervalBuilder::build_middle_interval)
/// - **Last interval**: Constructed using [`build_last_interval`](Grid1DIntervalBuilder::build_last_interval)
///
/// # Domain-Specific Behavior
///
/// The exact type of each sub-interval depends on the domain:
///
/// | Domain Type | First Interval | Middle Intervals | Last Interval |
/// |-------------|----------------|------------------|---------------|
/// | `[a, b]`    | `[p_0, p_1)`   | `[p_k, p_{k+1})` | `[p_{n-1}, p_n]` |
/// | `(a, b)`    | `(p_0, p_1]`   | `(p_k, p_{k+1}]` | `(p_{n-1}, p_n)` |
/// | `[a, b)`    | `[p_0, p_1)`   | `[p_k, p_{k+1})` | `[p_{n-1}, p_n)` |
/// | `(a, b]`    | `(p_0, p_1]`   | `(p_k, p_{k+1}]` | `(p_{n-1}, p_n]` |
///
/// This ensures that:
/// 1. The union of all sub-intervals exactly reconstructs the original domain
/// 2. No points are "lost" at the boundaries
/// 3. The partition is mathematically sound regardless of the domain type
///
/// # Mathematical Rationale
///
/// The choice of interval types ensures that every point in the original domain
/// belongs to exactly one sub-interval in the partition:
///
/// - For closed domains `[a, b]`, sub-intervals are left-closed right-open `[pᵢ, pᵢ₊₁)`,
///   except for the last which is fully closed `[pₙ₋₁, pₙ]` to capture the right boundary.
///   Partition points belong to the right sub-interval.
/// - For open domains `(a, b)`, sub-intervals are left-open right-closed `(pᵢ, pᵢ₊₁]`,
///   except for the last which is fully open `(pₙ₋₁, pₙ)` to respect the domain boundary.
///   Partition points belong to the left sub-interval.
/// - For lower-closed upper-open domains `[a, b)`, sub-intervals are left-closed right-open `[pᵢ, pᵢ₊₁)`,
///   except for the last which is left-closed right-open `[pₙ₋₁, pₙ)` to capture the left boundary and respect the open right boundary.
///   Partition points belong to the right sub-interval.
/// - For lower-open upper-closed domains `(a, b]`, sub-intervals are left-open right-closed `(pᵢ, pᵢ₊₁]`,
///   except for the last which is left-open right-closed `(pₙ₋₁, pₙ]` to respect the open left boundary and capture the right boundary.
///   Partition points belong to the left sub-interval.
///
/// # See Also
///
/// - [`crate::Grid1DTrait`]: The main trait for working with interval partitions
/// - [`SubIntervalInPartition`]: The enum representing different types of sub-intervals
pub trait Grid1DIntervalBuilder: IntervalFinitePositiveLengthTrait {
    /// The interval type used for the first sub-interval in a partition.
    ///
    /// This type must be convertible to both [`Interval`] and [`IntervalFinitePositiveLength`],
    /// ensuring it can represent the first interval with appropriate boundary semantics.
    type FirstIntervalInPartition: IntervalFinitePositiveLengthTrait<RealType = Self::RealType>
        + Into<Interval<Self::RealType>>
        + Into<IntervalFinitePositiveLength<Self::RealType>>;

    /// The interval type used for middle sub-intervals in a partition.
    ///
    /// A "middle" interval is any sub-interval that is neither the first nor the last
    /// in a multi-interval partition. The concrete type encodes the boundary semantics
    /// appropriate for interior intervals of the given domain type.
    ///
    /// The middle type mirrors the first type:
    /// - Left-closed domains (`[a, b]` and `[a, b)`) use `[pᵢ, pᵢ₊₁)` for middle intervals.
    /// - Left-open domains (`(a, b)` and `(a, b]`) use `(pᵢ, pᵢ₊₁]` for middle intervals.
    /// - Enum-domain `IntervalFinitePositiveLength` uses itself and dispatches at runtime.
    type MiddleIntervalInPartition: IntervalFinitePositiveLengthTrait<RealType = Self::RealType>
        + Into<Interval<Self::RealType>>
        + Into<IntervalFinitePositiveLength<Self::RealType>>;

    /// The interval type used for the last sub-interval in a partition.
    ///
    /// This type must be convertible to both [`Interval`] and [`IntervalFinitePositiveLength`],
    /// ensuring it can represent the last interval with appropriate boundary semantics.
    type LastIntervalInPartition: IntervalFinitePositiveLengthTrait<RealType = Self::RealType>
        + Into<Interval<Self::RealType>>
        + Into<IntervalFinitePositiveLength<Self::RealType>>;

    /// Wraps the domain itself into a [`SubIntervalInPartition::Single`] variant.
    ///
    /// This is the correct construction path when the partition has exactly one interval
    /// covering the entire domain. Consumes `self`.
    fn build_single_interval(self) -> SubIntervalInPartition<Self> {
        SubIntervalInPartition::Single(self)
    }

    /// Constructs the first sub-interval from the partition coordinates.
    ///
    /// Reads `coords[0]` and `coords[1]` as the lower and upper bounds.
    /// The concrete boundary semantics (open/closed) depend on the domain type;
    /// see the [`Grid1DIntervalBuilder`] table for details.
    ///
    /// # Arguments
    ///
    /// - `coords`: The full coordinate slice of the partition (at least 2 points).
    fn build_first_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
    ) -> SubIntervalInPartition<Self>;

    /// Constructs a middle sub-interval from the partition coordinates and the interval index.
    ///
    /// Reads `coords[interval_id]` and `coords[interval_id + 1]` as the lower and upper bounds.
    /// The concrete boundary semantics (open/closed) depend on the domain type;
    /// see the [`Grid1DIntervalBuilder`] table for details.
    ///
    /// # Arguments
    ///
    /// - `coords`: The full coordinate slice of the partition.
    /// - `interval_id`: The index of the interval (must satisfy `1 ≤ id < n - 1`).
    fn build_middle_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
        interval_id: &IntervalId,
    ) -> SubIntervalInPartition<Self>;

    /// Constructs the last sub-interval from the partition coordinates.
    ///
    /// Reads `coords[n-2]` and `coords[n-1]` as the lower and upper bounds,
    /// where `n = coords.len()`.
    /// The concrete boundary semantics (open/closed) depend on the domain type;
    /// see the [`Grid1DIntervalBuilder`] table for details.
    ///
    /// # Arguments
    ///
    /// - `coords`: The full coordinate slice of the partition (at least 2 points).
    fn build_last_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
    ) -> SubIntervalInPartition<Self>;
}

// Implementations of Grid1DIntervalBuilder for concrete interval types
#[duplicate_item(
    I                              type_first_interval            type_middle_interval           type_last_interval;
    [IntervalClosed]               [IntervalLowerClosedUpperOpen] [IntervalLowerClosedUpperOpen] [IntervalClosed];
    [IntervalOpen]                 [IntervalLowerOpenUpperClosed] [IntervalLowerOpenUpperClosed] [IntervalOpen];
    [IntervalLowerClosedUpperOpen] [IntervalLowerClosedUpperOpen] [IntervalLowerClosedUpperOpen] [IntervalLowerClosedUpperOpen];
    [IntervalLowerOpenUpperClosed] [IntervalLowerOpenUpperClosed] [IntervalLowerOpenUpperClosed] [IntervalLowerOpenUpperClosed];
)]
impl<RealType: RealScalar> Grid1DIntervalBuilder for I<RealType> {
    type FirstIntervalInPartition = type_first_interval<RealType>;
    type MiddleIntervalInPartition = type_middle_interval<RealType>;
    type LastIntervalInPartition = type_last_interval<RealType>;

    fn build_first_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
    ) -> SubIntervalInPartition<Self> {
        let coords = coords.as_ref();

        SubIntervalInPartition::First(
            <Self as Grid1DIntervalBuilder>::FirstIntervalInPartition::new(
                coords[0].clone(),
                coords[1].clone(),
            ),
        )
    }

    fn build_middle_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
        interval_id: &IntervalId,
    ) -> SubIntervalInPartition<Self> {
        let coords = coords.as_ref();
        let i = *interval_id.as_ref();

        SubIntervalInPartition::Middle(
            <Self as Grid1DIntervalBuilder>::MiddleIntervalInPartition::new(
                coords[i].clone(),
                coords[i + 1].clone(),
            ),
        )
    }

    fn build_last_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
    ) -> SubIntervalInPartition<Self> {
        let coords = coords.as_ref();
        let last_id = coords.len() - 1;

        SubIntervalInPartition::Last(
            <Self as Grid1DIntervalBuilder>::LastIntervalInPartition::new(
                coords[last_id - 1].clone(),
                coords[last_id].clone(),
            ),
        )
    }
}

impl<RealType: RealScalar> Grid1DIntervalBuilder for IntervalFinitePositiveLength<RealType> {
    type FirstIntervalInPartition = IntervalFinitePositiveLength<RealType>;
    type MiddleIntervalInPartition = IntervalFinitePositiveLength<RealType>;
    type LastIntervalInPartition = IntervalFinitePositiveLength<RealType>;

    fn build_first_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
    ) -> SubIntervalInPartition<Self> {
        let coords = coords.as_ref();
        let lower_bound = coords[0].clone();
        let upper_bound = coords[1].clone();

        let first_interval = match self {
            IntervalFinitePositiveLength::Open(_)
            | IntervalFinitePositiveLength::LowerOpenUpperClosed(_) => {
                IntervalLowerOpenUpperClosed::new(lower_bound, upper_bound).into()
            }
            IntervalFinitePositiveLength::Closed(_)
            | IntervalFinitePositiveLength::LowerClosedUpperOpen(_) => {
                IntervalLowerClosedUpperOpen::new(lower_bound, upper_bound).into()
            }
        };
        SubIntervalInPartition::First(first_interval)
    }

    fn build_middle_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
        interval_id: &IntervalId,
    ) -> SubIntervalInPartition<Self> {
        let coords = coords.as_ref();
        let i = *interval_id.as_ref();
        let lower_bound = coords[i].clone();
        let upper_bound = coords[i + 1].clone();

        let middle_interval = match self {
            IntervalFinitePositiveLength::Open(_)
            | IntervalFinitePositiveLength::LowerOpenUpperClosed(_) => {
                IntervalLowerOpenUpperClosed::new(lower_bound, upper_bound).into()
            }
            IntervalFinitePositiveLength::Closed(_)
            | IntervalFinitePositiveLength::LowerClosedUpperOpen(_) => {
                IntervalLowerClosedUpperOpen::new(lower_bound, upper_bound).into()
            }
        };
        SubIntervalInPartition::Middle(middle_interval)
    }

    fn build_last_interval(
        &self,
        coords: &Coords1D<Self::RealType>,
    ) -> SubIntervalInPartition<Self> {
        let coords = coords.as_ref();
        let last_id = coords.len() - 1;
        let lower_bound = coords[last_id - 1].clone();
        let upper_bound = coords[last_id].clone();

        let last_interval = match self {
            IntervalFinitePositiveLength::Open(_) => {
                IntervalOpen::new(lower_bound, upper_bound).into()
            }
            IntervalFinitePositiveLength::LowerOpenUpperClosed(_) => {
                IntervalLowerOpenUpperClosed::new(lower_bound, upper_bound).into()
            }
            IntervalFinitePositiveLength::Closed(_) => {
                IntervalClosed::new(lower_bound, upper_bound).into()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(_) => {
                IntervalLowerClosedUpperOpen::new(lower_bound, upper_bound).into()
            }
        };
        SubIntervalInPartition::Last(last_interval)
    }
}

// ============================================================================
// Grid1DTrait
// ============================================================================

/// A trait for types that can locate a point within their interval partition.
///
/// Implemented by all grid types. The containment rule at interior boundary points
/// depends on the domain type (see [`Grid1DIntervalBuilder`]).
pub trait FindIntervalIdOfPoint {
    /// The type of the 1D points.
    type Point1DType: RealScalar;

    /// Finds the ID of the interval containing `x`, returning `None` if outside the domain.
    ///
    /// The containment rule for boundary points depends on the domain type:
    /// - For left-closed sub-intervals (`[a,b]` and `[a,b)` domains), interior boundary point
    ///   `pₖ` belongs to the **right** interval (the one that starts at `pₖ`).
    /// - For right-closed sub-intervals (`(a,b)` and `(a,b]` domains), `pₖ` belongs to the
    ///   **left** interval (the one that ends at `pₖ`).
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1DUniform::new(
    ///     IntervalClosed::new(0.0, 2.0),
    ///     NumIntervals::try_new(2).unwrap()
    /// );
    ///
    /// assert_eq!(grid.find_interval_id_of_point(&1.0), Some(IntervalId::new(1))); // Boundary → right for [a,b]
    /// assert_eq!(grid.find_interval_id_of_point(&3.0), None); // Outside domain
    /// ```
    fn find_interval_id_of_point(&self, x: &Self::Point1DType) -> Option<IntervalId>;

    /// Finds all intervals that contain any of the given points.
    ///
    /// Returns `None` for points outside the domain.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 10.0),
    ///     NumIntervals::try_new(10).unwrap()
    /// );
    /// let points = vec![1.5, 3.2, 11.0];
    /// let intervals = grid.find_intervals_for_points(&points);
    ///
    /// assert_eq!(intervals[0], Some(IntervalId::new(1)));
    /// assert_eq!(intervals[1], Some(IntervalId::new(3)));
    /// assert_eq!(intervals[2], None); // Outside domain
    /// ```
    fn find_intervals_for_points(&self, points: &[Self::Point1DType]) -> Vec<Option<IntervalId>> {
        points
            .iter()
            .map(|point| self.find_interval_id_of_point(point))
            .collect()
    }

    /// Parallel version of [`find_intervals_for_points`](Self::find_intervals_for_points).
    ///
    /// Uses Rayon to parallelize the point location across multiple threads.
    /// This is beneficial when locating many points, typically more than ~1000.
    ///
    /// # Performance
    ///
    /// - **Uniform grids**: O(n/p) where n = number of points, p = number of threads
    /// - **Non-uniform grids**: O((n/p) log m) where m = number of intervals
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 10.0),
    ///     NumIntervals::try_new(100).unwrap()
    /// );
    ///
    /// // Generate many test points
    /// let points: Vec<f64> = (0..10_000).map(|i| i as f64 / 1000.0).collect();
    ///
    /// // Parallel point location
    /// let intervals = grid.find_intervals_for_points_parallel(&points);
    /// assert_eq!(intervals.len(), 10_000);
    /// ```
    fn find_intervals_for_points_parallel(
        &self,
        points: &[Self::Point1DType],
    ) -> Vec<Option<IntervalId>>
    where
        Self: Sync,
        Self::Point1DType: Send + Sync,
    {
        points
            .par_iter()
            .map(|point| self.find_interval_id_of_point(point))
            .collect()
    }
}

/// A trait for types that expose a contiguous range of interval indices.
///
/// `HasIntervalIdRange` provides a minimal, composable interface for querying which interval
/// indices are in scope — whether that is an entire grid or just a window into it.
/// It is a supertrait of [`Grid1DTrait`], so every full partition automatically
/// satisfies this interface, but lighter-weight types such as
/// [`Grid1DWindow`](crate::Grid1DWindow) and
/// [`IntervalIndexSpace1D`](crate::topology::IntervalIndexSpace1D) implement it too.
///
/// # Provided Method
///
/// [`num_intervals`](HasIntervalIdRange::num_intervals) is derived from `first_interval` and
/// `last_interval` and does not need to be overridden unless a more efficient
/// implementation is available.
///
/// # Implementors
///
/// | Type | `first_interval` | `last_interval` |
/// |------|-----------------|----------------|
/// | [`Grid1DUniform`](crate::Grid1DUniform) | always `0` | `n - 1` |
/// | [`Grid1DNonUniform`](crate::Grid1DNonUniform) | always `0` | `n - 1` |
/// | [`Grid1D`](crate::Grid1D) | always `0` | `n - 1` |
/// | [`Grid1DUnion`](crate::Grid1DUnion) | always `0` | `n - 1` |
/// | [`Grid1DWindow`](crate::Grid1DWindow) | window start | window end |
/// | [`IntervalIndexSpace1D`](crate::topology::IntervalIndexSpace1D) | always `0` | `n - 1` |
///
/// # Example
///
/// ```rust
/// use grid1d::{
///     Grid1D, HasIntervalIdRange,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
///
/// assert_eq!(grid.first_interval_id().as_ref(), &0);
/// assert_eq!(grid.last_interval_id().as_ref(), &3);
/// assert_eq!(grid.num_intervals().as_ref(), &4);
/// ```
///
/// # See Also
///
/// - [`HasIntervals`]: Extends `HasIntervalIdRange` with interval lookup (`interval()`) and iteration
/// - [`Grid1DTrait`]: Full partition interface, a superset of both `HasIntervalIdRange` and `HasIntervals`
/// - [`crate::Grid1DWindow`]: A windowed view that narrows the range
pub trait HasIntervalIdRange {
    /// Returns the [`IntervalId`] of the first interval in scope.
    fn first_interval_id(&self) -> IntervalId;

    /// Returns the [`IntervalId`] of the last interval in scope.
    fn last_interval_id(&self) -> IntervalId;

    /// Returns the total number of intervals in scope.
    ///
    /// The default implementation derives this from `first_interval_id` and `last_interval_id`:
    ///
    /// ```text
    /// num_intervals = last_interval_id - first_interval_id + 1
    /// ```
    #[inline]
    fn num_intervals(&self) -> NumIntervals {
        NumIntervals::try_new(
            self.last_interval_id().into_inner() - self.first_interval_id().into_inner() + 1,
        )
        .unwrap()
    }

    /// Returns `true` if `interval_id` falls within this range.
    #[inline]
    fn range_contains_interval_id(&self, interval_id: &IntervalId) -> bool {
        let first_id = self.first_interval_id();
        let last_id = self.last_interval_id();

        interval_id >= &first_id && interval_id <= &last_id
    }
}

/// A trait for types that expose a contiguous range of coordinate indices.
///
/// `HasCoordIdRange` provides a minimal, composable interface for querying which coordinate
/// indices are in scope — whether that is an entire grid or just a window into it.
/// All full-partition types implement this trait, as does
/// [`CoordIndexSpace1D`](crate::topology::CoordIndexSpace1D).
///
/// # Provided Methods
///
/// [`num_coord_ids`](HasCoordIdRange::num_coord_ids) and
/// [`range_contains_coord_id`](HasCoordIdRange::range_contains_coord_id) are derived
/// from `first_coord_id` and `last_coord_id` and do not need to be overridden unless
/// a more efficient implementation is available.
///
/// # Implementors
///
/// | Type | `first_coord_id` | `last_coord_id` |
/// |------|-----------------|----------------|
/// | [`Grid1DUniform`](crate::Grid1DUniform) | `0` (non-periodic) or `1` (lower-open periodic) | `n` (non-periodic) or `n-1` (periodic) |
/// | [`Grid1DNonUniform`](crate::Grid1DNonUniform) | always `0` | `n` |
/// | [`Grid1D`](crate::Grid1D) | depends on variant | depends on variant |
/// | [`CoordIndexSpace1D`](crate::topology::CoordIndexSpace1D) | `0` (non-periodic) or `1` (lower-open periodic) | `n` (non-periodic) or `n-1` (periodic) |
///
/// # Example
///
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// let num_intervals = NumIntervals::try_new(4).unwrap();
///
/// let grid = Grid1DUniform::new(IntervalClosed::new(0.0, 1.0), num_intervals);
/// assert_eq!(grid.first_coord_id().as_ref(), &0);
/// assert_eq!(grid.last_coord_id().as_ref(), &4);
/// assert_eq!(grid.num_coord_ids().as_ref(), &5);
///
/// let grid_periodic_upper_open = Grid1DUniform::new_periodic(
///     IntervalLowerClosedUpperOpen::new(0.0, 1.0),
///     num_intervals,
/// );
/// assert_eq!(grid_periodic_upper_open.first_coord_id().as_ref(), &0);
/// assert_eq!(grid_periodic_upper_open.last_coord_id().as_ref(), &3);
/// assert_eq!(grid_periodic_upper_open.num_coord_ids().as_ref(), &4);
///
/// let grid_periodic_lower_open = Grid1DUniform::new_periodic(
///     IntervalLowerOpenUpperClosed::new(0.0, 1.0),
///     num_intervals,
/// );
/// assert_eq!(grid_periodic_lower_open.first_coord_id().as_ref(), &1);
/// assert_eq!(grid_periodic_lower_open.last_coord_id().as_ref(), &4);
/// assert_eq!(grid_periodic_lower_open.num_coord_ids().as_ref(), &4);
/// ```
///
/// # See Also
///
/// - [`HasIntervalIdRange`]: Analogous trait for the interval index range
/// - [`HasCoords1D`]: Provides direct access to the coordinate values
/// - [`Grid1DTrait`]: Full partition interface
pub trait HasCoordIdRange {
    /// Returns the [`CoordId`] of the first coordinate in scope.
    fn first_coord_id(&self) -> CoordId;

    /// Returns the [`CoordId`] of the last coordinate in scope.
    fn last_coord_id(&self) -> CoordId;

    /// Returns the total number of coordinates in scope.
    ///
    /// The default implementation derives this from `first_coord_id` and `last_coord_id`:
    ///
    /// ```text
    /// num_coord_ids = last_coord_id - first_coord_id + 1
    /// ```
    #[inline]
    fn num_coord_ids(&self) -> PositiveNumPoints1D {
        PositiveNumPoints1D::try_new(
            self.last_coord_id().into_inner() - self.first_coord_id().into_inner() + 1,
        )
        .unwrap()
    }

    /// Returns `true` if `coord_id` falls within this range.
    #[inline]
    fn range_contains_coord_id(&self, coord_id: &CoordId) -> bool {
        let first_id = self.first_coord_id();
        let last_id = self.last_coord_id();

        coord_id >= &first_id && coord_id <= &last_id
    }
}

/// A trait for types that expose individual intervals by ID and support iteration over them.
///
/// `HasIntervals` extends [`HasIntervalIdRange`] with the ability to retrieve a concrete interval
/// value given its [`IntervalId`], and to iterate over all intervals in scope. It sits between
/// [`HasIntervalIdRange`] (which only exposes the index range) and [`Grid1DTrait`] (which
/// additionally exposes coordinates, a domain, and point-location).
///
/// This intermediate position makes `HasIntervals` the right bound for generic code that needs
/// to inspect interval geometry (e.g. containment checks, length queries) without requiring a
/// full partition — most notably for windowed views via
/// [`Grid1DWindow`](crate::Grid1DWindow).
///
/// # Associated Type
///
/// [`IntervalType`](HasIntervals::IntervalType) must implement [`IntervalTrait`]. For all
/// [`Grid1DTrait`] implementors it is [`SubIntervalInPartition`], which is a concrete enum
/// that encodes the open/closed boundary semantics of each sub-interval position (first, middle,
/// last).
///
/// # Implementors
///
/// | Type | `IntervalType` | Notes |
/// |------|---------------|-------|
/// | [`Grid1DUniform`](crate::Grid1DUniform) | [`SubIntervalInPartition`] | computed on the fly from the uniform step |
/// | [`Grid1DNonUniform`](crate::Grid1DNonUniform) | [`SubIntervalInPartition`] | constructed from stored coordinates |
/// | [`Grid1D`](crate::Grid1D) | [`SubIntervalInPartition`] | delegates to the active variant |
/// | [`Grid1DUnion`](crate::Grid1DUnion) | [`SubIntervalInPartition`] | delegates to the unified grid |
/// | [`Grid1DWindow`](crate::Grid1DWindow) | same as the underlying partition | delegates via the absolute ID |
///
/// # Example: interval geometry
///
/// ```rust
/// use grid1d::{Grid1D, HasIntervals, FindIntervalIdOfPoint, intervals::*, scalars::NumIntervals};
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(
///     IntervalClosed::new(0.0, 3.0),
///     NumIntervals::try_new(3).unwrap(),
/// );
///
/// let id = grid.find_interval_id_of_point(&1.5).unwrap();
/// let interval = grid.interval(&id);
/// assert!(interval.contains_point(&1.5));
/// ```
///
/// # Example: iterating over all intervals
///
/// ```rust
/// use grid1d::{Grid1D, HasIntervals, intervals::*, scalars::NumIntervals};
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(
///     IntervalClosed::new(0.0, 2.0),
///     NumIntervals::try_new(2).unwrap(),
/// );
///
/// for (id, interval) in grid.iter_intervals() {
///     println!("Interval {}: {:?}", id.as_ref(), interval);
/// }
/// ```
///
/// # See Also
///
/// - [`HasIntervalIdRange`]: Lighter-weight supertrait that only exposes the index range
/// - [`Grid1DTrait`]: Full partition interface, a superset of `HasIntervals`
/// - [`SubIntervalInPartition`]: The concrete interval type returned by partition implementors
pub trait HasIntervals: HasIntervalIdRange {
    /// The concrete interval type returned by [`interval`](HasIntervals::interval).
    ///
    /// Must implement [`IntervalTrait`]. For all [`Grid1DTrait`] implementors this
    /// is [`SubIntervalInPartition`], which encodes the open/closed boundary semantics for
    /// first, middle, and last sub-intervals.
    type IntervalType: IntervalTrait;

    /// Returns the interval corresponding to the given `interval_id`.
    ///
    /// # Panics
    ///
    /// In debug builds, panics if `interval_id` is outside the range
    /// `[first_interval_id(), last_interval_id()]`.
    fn interval(&self, interval_id: &IntervalId) -> Self::IntervalType;

    /// Returns an iterator over all `(IntervalId, IntervalType)` pairs in scope,
    /// in ascending ID order.
    ///
    /// The iterator visits every interval from [`first_interval_id`](HasIntervalIdRange::first_interval_id)
    /// to [`last_interval_id`](HasIntervalIdRange::last_interval_id) inclusive.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::NumIntervals};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 2.0),
    ///     NumIntervals::try_new(2).unwrap()
    /// );
    ///
    /// for (id, interval) in grid.iter_intervals() {
    ///     println!("Interval {}: {:?}", id.as_ref(), interval);
    /// }
    /// ```
    fn iter_intervals(&self) -> impl Iterator<Item = (IntervalId, Self::IntervalType)> + '_ {
        let first_id = self.first_interval_id().into_inner();
        let last_id = self.last_interval_id().into_inner();
        (first_id..=last_id).map(move |i| {
            let id = IntervalId::new(i);
            let interval = self.interval(&id);
            (id, interval)
        })
    }
}

impl<G: Grid1DTrait> HasIntervals for G {
    type IntervalType = SubIntervalInPartition<G::Domain1D>;

    /// Returns the `i`-th interval as [`SubIntervalInPartition`].
    ///
    /// # Panics
    ///
    /// In debug mode, panics if the interval ID is not in the interval range of the partition.
    fn interval(&self, i: &IntervalId) -> SubIntervalInPartition<G::Domain1D> {
        debug_assert!(
            self.range_contains_interval_id(i),
            "Interval ID {:?} is out of range [{:?}, {:?}]",
            i,
            self.first_interval_id(),
            self.last_interval_id()
        );

        let domain = self.domain();

        // Case with only one interval: it's the domain itself
        let num_intervals = *self.num_intervals().as_ref();
        if num_intervals == 1 {
            return domain.clone().build_single_interval();
        }

        let coords = self.coords();

        let i_usize = *i.as_ref();
        if i_usize == 0 {
            // First interval
            domain.build_first_interval(coords)
        } else if i_usize == num_intervals - 1 {
            // Last interval
            domain.build_last_interval(coords)
        } else {
            // Middle interval
            domain.build_middle_interval(coords, i)
        }
    }
}

/// A trait for modeling a partition of an interval as a set of adjacent non-overlapping intervals.
///
/// This trait provides a unified interface for partitioning any type of interval with finite,
/// positive length into sub-intervals. It supports closed, open, and semi-open intervals.
///
/// # Partition Rules
///
/// Given points `{p_0, p_1, ..., p_n}` that define partition boundaries:
///
/// | Domain Type | First Interval `I_0` | Middle Intervals `I_k` | Last Interval `I_{n-1}` |
/// |-------------|----------------------|------------------------|-------------------------|
/// | `[a, b]`    | `[p_0, p_1)`         | `[p_k, p_{k+1})`       | `[p_{n-1}, p_n]`        |
/// | `(a, b)`    | `(p_0, p_1]`         | `(p_k, p_{k+1}]`       | `(p_{n-1}, p_n)`        |
/// | `[a, b)`    | `[p_0, p_1)`         | `[p_k, p_{k+1})`       | `[p_{n-1}, p_n)`        |
/// | `(a, b]`    | `(p_0, p_1]`         | `(p_k, p_{k+1}]`       | `(p_{n-1}, p_n]`        |
///
/// # Mathematical Properties
///
/// - **Completeness**: The union of all sub-intervals equals the original domain
/// - **Non-overlapping**: Sub-intervals are pairwise disjoint
/// - **Point uniqueness**: Every point in the domain belongs to exactly one sub-interval
///
/// # Example
///
/// ```rust
/// use grid1d::{*, intervals::*};
/// use try_create::TryNew;
///
/// let grid = Grid1D::uniform(
///     IntervalClosed::new(0.0, 10.0),
///     NumIntervals::try_new(5).unwrap()
/// );
///
/// // Basic partition properties
/// assert_eq!(grid.num_intervals().as_ref(), &5);
/// assert_eq!(grid.num_points().as_ref(), &6);
///
/// // Point location
/// let interval_id = grid.find_interval_id_of_point(&3.5).unwrap();
/// assert_eq!(interval_id.as_ref(), &1);
///
/// // Interval access
/// let interval = grid.interval(&interval_id);
/// ```
///
/// # See Also
///
/// - [`HasDomain1D`]: Provides domain access
/// - [`HasCoords1D`]: Provides coordinate access
/// - [`HasIntervals`]: Provides interval lookup and iteration (supertrait of `Grid1DTrait`)
/// - [`Grid1DIntervalBuilder`]: Defines sub-interval construction rules
pub trait Grid1DTrait:
    HasCoords1D
    + HasIntervals<IntervalType = SubIntervalInPartition<Self::Domain1D>>
    + FindIntervalIdOfPoint<Point1DType = <Self as HasCoords1D>::CoordType>
    + HasDomain1D<Domain1D: Grid1DIntervalBuilder<RealType = <Self as HasCoords1D>::CoordType>>
    + Serialize
    + for<'a> Deserialize<'a>
where
    Self::Domain1D: Grid1DIntervalBuilder<RealType = Self::CoordType>,
{
    /// Returns the index spaces associated with this partition.
    ///
    /// The returned [`Grid1DIndexSpaces`] bundles the [`IntervalIndexSpace1D`](crate::topology::IntervalIndexSpace1D)
    /// (one entry per sub-interval) and the [`CoordIndexSpace1D`](crate::topology::CoordIndexSpace1D) (one entry per
    /// partition vertex), both sharing the same [`Topology1D`](crate::topology::Topology1D).
    ///
    /// Use this to query topology information or perform adjacency lookups:
    ///
    /// ```rust
    /// use grid1d::{Grid1D, Grid1DTrait, IndexSpace1D, Topology1D, intervals::*, scalars::NumIntervals};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
    /// let spaces = grid.index_spaces();
    /// assert_eq!(spaces.interval_index_space().topology(), &Topology1D::RealLine);
    /// assert_eq!(*spaces.interval_index_space().count().as_ref(), 4);
    /// assert_eq!(*spaces.coord_index_space().count().as_ref(), 5);
    /// ```
    fn index_spaces(&self) -> &Grid1DIndexSpaces;

    /// Returns the left neighbor interval ID, if it exists.
    ///
    /// For periodic topology, the left neighbor of interval `0` is the last
    /// interval. For non-periodic topology, interval `0` has no left neighbor.
    fn get_interval_id_left_neighbor(&self, interval_id: &IntervalId) -> Option<IntervalId> {
        self.index_spaces()
            .interval_index_space()
            .left_neighbor(interval_id)
    }

    /// Returns the right neighbor interval ID, if it exists.
    ///
    /// For periodic topology, the right neighbor of the last interval is
    /// interval `0`. For non-periodic topology, the last interval has no
    /// right neighbor.
    fn get_interval_id_right_neighbor(&self, interval_id: &IntervalId) -> Option<IntervalId> {
        self.index_spaces()
            .interval_index_space()
            .right_neighbor(interval_id)
    }

    /// Returns the neighbor interval ID on the requested [`Side`], if it exists.
    fn get_interval_id_neighbor(
        &self,
        interval_id: &IntervalId,
        side: &Side,
    ) -> Option<IntervalId> {
        self.index_spaces()
            .interval_index_space()
            .neighbor(interval_id, side)
    }

    /// Returns the interval length at the specified index.
    ///
    /// # Panics
    ///
    /// In debug mode, panics if the interval ID is out of bounds.
    fn interval_length(&self, interval_id: &IntervalId) -> PositiveRealScalar<Self::Point1DType> {
        let coords = self.coords().as_ref();
        let i = *interval_id.as_ref();
        debug_assert_lt!(i, coords.len() - 1, "Interval ID out of bounds");

        PositiveRealScalar::try_new(coords[i + 1].clone() - &coords[i])
            .expect("Non-positive interval length!")
    }

    /// Returns the maximum interval length in the partition.
    fn max_interval_length(&self) -> PositiveRealScalar<Self::Point1DType> {
        (0..*self.num_intervals().as_ref())
            .map(|i| self.interval_length(&IntervalId::new(i)))
            .max_by(|a, b| a.as_ref().partial_cmp(b.as_ref()).unwrap())
            .unwrap()
    }

    /// Returns the minimum interval length in the partition.
    fn min_interval_length(&self) -> PositiveRealScalar<Self::Point1DType> {
        (0..*self.num_intervals().as_ref())
            .map(|i| self.interval_length(&IntervalId::new(i)))
            .min_by(|a, b| a.as_ref().partial_cmp(b.as_ref()).unwrap())
            .unwrap()
    }

    /// Returns the ratio between maximum and minimum interval lengths.
    ///
    /// A value close to 1.0 indicates a nearly uniform grid,
    /// while larger values indicate significant non-uniformity.
    fn uniformity_ratio(&self) -> PositiveRealScalar<Self::Point1DType> {
        let max_len = self.max_interval_length();
        let min_len = self.min_interval_length();

        PositiveRealScalar::try_new(max_len.as_ref().clone() / min_len.as_ref())
            .expect("Uniformity ratio must be positive")
    }

    /// Returns all intervals that intersect with the given domain.
    ///
    /// Only intervals with positive-length intersections are returned.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 10.0),
    ///     NumIntervals::try_new(10).unwrap()
    /// );
    ///
    /// let subdomain = IntervalClosed::new(2.5, 4.5);
    /// let intersections = grid.intervals_in_intersection(&subdomain);
    ///
    /// // Returns intervals 2, 3, 4 with their intersection regions
    /// assert_eq!(intersections.len(), 3);
    /// ```
    fn intervals_in_intersection<
        IntervalType: IntervalFinitePositiveLengthTrait<RealType = Self::Point1DType>,
    >(
        &self,
        domain_in: &IntervalType,
    ) -> Vec<(IntervalId, IntervalFinitePositiveLength<Self::Point1DType>)> {
        // Find the positive-length overlap between the grid's domain and the input domain
        let overlap = match self.domain().intersection(domain_in) {
            Some(Interval::FiniteLength(IntervalFiniteLength::PositiveLength(p))) => p,
            _ => return Vec::new(),
        };

        // Find the start and end interval IDs that cover the overlap
        let id_min = self
            .find_interval_id_of_point(overlap.lower_bound_value())
            .expect("Lower bound of overlap must be in domain");
        let id_max = self
            .find_interval_id_of_point(overlap.upper_bound_value())
            .expect("Upper bound of overlap must be in domain");

        let id_min_as_usize = *id_min.as_ref();
        let id_max_as_usize = *id_max.as_ref();

        let n_max_intersecting_intervals = id_max_as_usize - id_min_as_usize + 1;
        let mut intersecting_intervals = Vec::with_capacity(n_max_intersecting_intervals);

        for i in id_min_as_usize..=id_max_as_usize {
            let current_id = IntervalId::new(i);
            let grid_interval = self.interval(&current_id);

            if let Some(Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                intersection,
            ))) = grid_interval.intersection(&overlap)
            {
                intersecting_intervals.push((current_id, intersection));
            }
        }

        debug_assert!(
            !intersecting_intervals.is_empty(),
            "Expected at least one interval with positive length in the intersection"
        );

        intersecting_intervals
    }

    /// The type of grid resulting from uniform refinement.
    type UniformlyRefinedGrid1DType: Grid1DTrait<
            CoordType = Self::CoordType,
            Point1DType = Self::Point1DType,
            Domain1D = Self::Domain1D,
        >;

    /// Performs uniform refinement where all intervals are subdivided equally.
    ///
    /// # Parameters
    ///
    /// - `num_extra_points_each_interval`: Number of additional points to insert in each interval
    ///
    /// # Returns
    ///
    /// A [`crate::Grid1DUniformRefinement`] containing the refined grid
    /// and bidirectional mappings.
    fn refine_uniform(
        self,
        num_extra_points_each_interval: &PositiveNumPoints1D,
    ) -> Grid1DUniformRefinement<Self>;

    /// Performs selective refinement of specified intervals.
    ///
    /// Only the intervals in the map are refined; others remain unchanged.
    /// The [`BTreeMap`] ensures unique interval IDs and ordered iteration.
    ///
    /// # Parameters
    ///
    /// - `intervals_to_refine`: Map from interval IDs to number of extra points
    ///
    /// # Panics
    ///
    /// In debug mode, panics if any interval ID is out of bounds.
    fn refine(
        self,
        intervals_to_refine: &BTreeMap<IntervalId, PositiveNumPoints1D>,
    ) -> Grid1DNonUniformRefinement<Self>;
}

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

#[cfg(test)]
mod tests {
    use super::*;
    use crate::intervals::{Contains, bounded::IntervalFromBounds};
    use crate::{Grid1D, Grid1DNonUniform, Grid1DUniform, IndexSpace1D, Topology1D};
    use try_create::TryNew;

    #[test]
    fn test_has_domain_1d() {
        let domain = IntervalClosed::new(0.0, 10.0);
        let grid = Grid1DUniform::new(domain.clone(), NumIntervals::try_new(5).unwrap());

        assert_eq!(grid.domain(), &domain);
    }

    #[test]
    fn test_has_coords_1d() {
        let grid = Grid1DUniform::new(
            IntervalClosed::new(0.0, 1.0),
            NumIntervals::try_new(4).unwrap(),
        );

        assert_eq!(grid.num_points().as_ref(), &5);
        assert_eq!(grid.coords().len(), 5);
    }

    #[test]
    fn test_has_coords_1d_coord_uniform() {
        let grid = Grid1DUniform::new(
            IntervalClosed::new(0.0, 1.0),
            NumIntervals::try_new(4).unwrap(),
        );

        assert_eq!(grid.coord(CoordId::new(0)), &0.0);
        assert_eq!(grid.coord(CoordId::new(2)), &0.5);
        assert_eq!(grid.coord(CoordId::new(4)), &1.0);
    }

    #[test]
    fn test_has_coords_1d_coord_non_uniform() {
        let grid = Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(
            Coords1D::try_from(sorted_vec::partial::SortedSet::from_unsorted(vec![
                0.0, 0.25, 1.0,
            ]))
            .unwrap(),
        )
        .unwrap();

        assert_eq!(grid.coord(CoordId::new(0)), &0.0);
        assert_eq!(grid.coord(CoordId::new(1)), &0.25);
        assert_eq!(grid.coord(CoordId::new(2)), &1.0);
    }

    #[test]
    fn test_build_interval_in_partition_closed() {
        let partition_domain = IntervalClosed::new(0.0, 3.0);

        let coords = Coords1D::new_uniform(&partition_domain, &NumIntervals::try_new(3).unwrap());

        let first = partition_domain.build_first_interval(&coords);
        let middle = partition_domain.build_middle_interval(&coords, &IntervalId::new(1));
        let last = partition_domain.build_last_interval(&coords);

        // First sub-interval for closed domain is [a, b)
        assert!(first.contains_point(&0.0));
        assert!(!first.contains_point(&1.0));

        // Middle intervals are always [a, b)
        assert!(middle.contains_point(&1.0));
        assert!(!middle.contains_point(&2.0));

        // Last sub-interval for closed domain is [a, b]
        assert!(last.contains_point(&2.0));
        assert!(last.contains_point(&3.0));
    }

    #[test]
    fn test_build_interval_in_partition_open() {
        let partition_domain = IntervalOpen::new(0.0, 3.0);

        let coords = Coords1D::new_uniform(&partition_domain, &NumIntervals::try_new(3).unwrap());

        let first = partition_domain.build_first_interval(&coords);
        let middle = partition_domain.build_middle_interval(&coords, &IntervalId::new(1));
        let last = partition_domain.build_last_interval(&coords);

        // First interval for open domain is (a, b]
        assert!(!first.contains_point(&0.0));
        assert!(first.contains_point(&1.0));

        // Middle intervals are always (a, b]
        assert!(!middle.contains_point(&1.0));
        assert!(middle.contains_point(&2.0));

        // Last interval for open domain is (a, b)
        assert!(!last.contains_point(&2.0));
        assert!(!last.contains_point(&3.0));
    }

    #[test]
    fn test_build_interval_in_partition_lower_closed_upper_open() {
        let partition_domain = IntervalLowerClosedUpperOpen::new(0.0, 3.0);

        let coords = Coords1D::new_uniform(&partition_domain, &NumIntervals::try_new(3).unwrap());

        let first = partition_domain.build_first_interval(&coords);
        let middle = partition_domain.build_middle_interval(&coords, &IntervalId::new(1));
        let last = partition_domain.build_last_interval(&coords);

        // First interval for lower closed upper open domain is [a, b)
        assert!(first.contains_point(&0.0));
        assert!(!first.contains_point(&1.0));

        // Middle intervals are always [a, b)
        assert!(middle.contains_point(&1.0));
        assert!(!middle.contains_point(&2.0));

        // Last interval for lower closed upper open domain is [a, b)
        assert!(last.contains_point(&2.0));
        assert!(!last.contains_point(&3.0));
    }

    #[test]
    fn test_build_interval_in_partition_lower_open_upper_closed() {
        let partition_domain = IntervalLowerOpenUpperClosed::new(0.0, 3.0);

        let coords = Coords1D::new_uniform(&partition_domain, &NumIntervals::try_new(3).unwrap());

        let first = partition_domain.build_first_interval(&coords);
        let middle = partition_domain.build_middle_interval(&coords, &IntervalId::new(1));
        let last = partition_domain.build_last_interval(&coords);

        // First interval for lower open upper closed domain is (a, b]
        assert!(!first.contains_point(&0.0));
        assert!(first.contains_point(&1.0));

        // Middle intervals are always (a, b]
        assert!(!middle.contains_point(&1.0));
        assert!(middle.contains_point(&2.0));

        // Last interval for lower open upper closed domain is (a, b]
        assert!(!last.contains_point(&2.0));
        assert!(last.contains_point(&3.0));
    }

    #[test]
    fn test_build_interval_in_partition_enum_first_middle_open_branch() {
        let coords = Coords1D::new_uniform(
            &IntervalClosed::new(0.0, 3.0),
            &NumIntervals::try_new(3).unwrap(),
        );

        let domain_enum: IntervalFinitePositiveLength<f64> = IntervalOpen::new(0.0, 3.0).into();

        let first = domain_enum.build_first_interval(&coords);
        let middle = domain_enum.build_middle_interval(&coords, &IntervalId::new(1));

        match first {
            SubIntervalInPartition::First(interval) => {
                assert!(matches!(
                    interval,
                    IntervalFinitePositiveLength::LowerOpenUpperClosed(_)
                ));
            }
            _ => panic!("Expected first interval"),
        }

        match middle {
            SubIntervalInPartition::Middle(interval) => {
                assert!(matches!(
                    interval,
                    IntervalFinitePositiveLength::LowerOpenUpperClosed(_)
                ));
            }
            _ => panic!("Expected middle interval"),
        }
    }

    #[test]
    fn test_build_interval_in_partition_enum_first_middle_closed_branch() {
        let coords = Coords1D::new_uniform(
            &IntervalClosed::new(0.0, 3.0),
            &NumIntervals::try_new(3).unwrap(),
        );

        let domain_enum: IntervalFinitePositiveLength<f64> = IntervalClosed::new(0.0, 3.0).into();

        let first = domain_enum.build_first_interval(&coords);
        let middle = domain_enum.build_middle_interval(&coords, &IntervalId::new(1));

        match first {
            SubIntervalInPartition::First(interval) => {
                assert!(matches!(
                    interval,
                    IntervalFinitePositiveLength::LowerClosedUpperOpen(_)
                ));
            }
            _ => panic!("Expected first interval"),
        }

        match middle {
            SubIntervalInPartition::Middle(interval) => {
                assert!(matches!(
                    interval,
                    IntervalFinitePositiveLength::LowerClosedUpperOpen(_)
                ));
            }
            _ => panic!("Expected middle interval"),
        }
    }

    #[test]
    fn test_build_interval_in_partition_enum_last_open_branch() {
        let coords = Coords1D::new_uniform(
            &IntervalClosed::new(0.0, 3.0),
            &NumIntervals::try_new(3).unwrap(),
        );

        let domain_enum: IntervalFinitePositiveLength<f64> = IntervalOpen::new(0.0, 3.0).into();
        let last = domain_enum.build_last_interval(&coords);

        match last {
            SubIntervalInPartition::Last(interval) => {
                assert!(matches!(interval, IntervalFinitePositiveLength::Open(_)));
            }
            _ => panic!("Expected last interval"),
        }
    }

    #[test]
    fn test_build_interval_in_partition_enum_last_lower_open_upper_closed_branch() {
        let coords = Coords1D::new_uniform(
            &IntervalClosed::new(0.0, 3.0),
            &NumIntervals::try_new(3).unwrap(),
        );

        let domain_enum: IntervalFinitePositiveLength<f64> =
            IntervalLowerOpenUpperClosed::new(0.0, 3.0).into();
        let last = domain_enum.build_last_interval(&coords);

        match last {
            SubIntervalInPartition::Last(interval) => {
                assert!(matches!(
                    interval,
                    IntervalFinitePositiveLength::LowerOpenUpperClosed(_)
                ));
            }
            _ => panic!("Expected last interval"),
        }
    }

    #[test]
    fn test_build_interval_in_partition_enum_last_closed_branch() {
        let coords = Coords1D::new_uniform(
            &IntervalClosed::new(0.0, 3.0),
            &NumIntervals::try_new(3).unwrap(),
        );

        let domain_enum: IntervalFinitePositiveLength<f64> = IntervalClosed::new(0.0, 3.0).into();
        let last = domain_enum.build_last_interval(&coords);

        match last {
            SubIntervalInPartition::Last(interval) => {
                assert!(matches!(interval, IntervalFinitePositiveLength::Closed(_)));
            }
            _ => panic!("Expected last interval"),
        }
    }

    #[test]
    fn test_build_interval_in_partition_enum_last_lower_closed_upper_open_branch() {
        let coords = Coords1D::new_uniform(
            &IntervalClosed::new(0.0, 3.0),
            &NumIntervals::try_new(3).unwrap(),
        );

        let domain_enum: IntervalFinitePositiveLength<f64> =
            IntervalLowerClosedUpperOpen::new(0.0, 3.0).into();
        let last = domain_enum.build_last_interval(&coords);

        match last {
            SubIntervalInPartition::Last(interval) => {
                assert!(matches!(
                    interval,
                    IntervalFinitePositiveLength::LowerClosedUpperOpen(_)
                ));
            }
            _ => panic!("Expected last interval"),
        }
    }

    #[test]
    fn test_interval_partition_num_intervals() {
        let grid = Grid1D::uniform(
            IntervalClosed::new(0.0, 10.0),
            NumIntervals::try_new(5).unwrap(),
        );

        assert_eq!(grid.num_intervals().as_ref(), &5);
        assert_eq!(grid.num_points().as_ref(), &6);
    }

    #[test]
    fn test_interval_partition_find_interval() {
        let grid = Grid1DUniform::new(
            IntervalClosed::new(0.0, 4.0),
            NumIntervals::try_new(4).unwrap(),
        );

        // Interior points
        assert_eq!(
            grid.find_interval_id_of_point(&0.5).unwrap(),
            IntervalId::new(0)
        );
        assert_eq!(
            grid.find_interval_id_of_point(&1.5).unwrap(),
            IntervalId::new(1)
        );
        assert_eq!(
            grid.find_interval_id_of_point(&3.5).unwrap(),
            IntervalId::new(3)
        );

        // Boundary points belong to RIGHT interval for [a,b] domain (sub-intervals are [pₖ, pₖ₊₁))
        assert_eq!(
            grid.find_interval_id_of_point(&1.0).unwrap(),
            IntervalId::new(1)
        );
        assert_eq!(
            grid.find_interval_id_of_point(&2.0).unwrap(),
            IntervalId::new(2)
        );

        // Domain boundaries
        assert_eq!(
            grid.find_interval_id_of_point(&0.0).unwrap(),
            IntervalId::new(0)
        );
        assert_eq!(
            grid.find_interval_id_of_point(&4.0).unwrap(),
            IntervalId::new(3)
        );
    }

    #[test]
    fn test_interval_partition_try_find() {
        let grid = Grid1DUniform::new(
            IntervalClosed::new(0.0, 2.0),
            NumIntervals::try_new(2).unwrap(),
        );

        assert_eq!(
            grid.find_interval_id_of_point(&1.0),
            Some(IntervalId::new(1)) // Boundary → right for [a,b] domain
        );
        assert_eq!(grid.find_interval_id_of_point(&-1.0), None);
        assert_eq!(grid.find_interval_id_of_point(&3.0), None);
    }

    #[test]
    fn test_interval_partition_uniformity() {
        let uniform_grid = Grid1DUniform::new(
            IntervalClosed::new(0.0, 1.0),
            NumIntervals::try_new(10).unwrap(),
        );

        assert_eq!(uniform_grid.uniformity_ratio().as_ref(), &1.0);
        assert_eq!(
            uniform_grid.min_interval_length(),
            uniform_grid.max_interval_length()
        );
    }

    #[test]
    fn test_side_helpers() {
        assert!(Side::Left.is_left());
        assert!(!Side::Left.is_right());
        assert!(Side::Right.is_right());
        assert!(!Side::Right.is_left());
    }

    #[test]
    fn test_interval_neighbors_real_line_no_wraparound() {
        let grid = Grid1DUniform::new(
            IntervalClosed::new(0.0, 1.0),
            NumIntervals::try_new(4).unwrap(),
        );

        let grid_index_spaces = grid.index_spaces();
        assert_eq!(
            grid_index_spaces.interval_index_space().topology(),
            &Topology1D::RealLine
        );
        assert_eq!(
            grid_index_spaces.coord_index_space().topology(),
            &Topology1D::RealLine
        );

        assert_eq!(
            grid.get_interval_id_left_neighbor(&IntervalId::new(0)),
            None
        );
        assert_eq!(
            grid.get_interval_id_right_neighbor(&IntervalId::new(0)),
            Some(IntervalId::new(1))
        );

        assert_eq!(
            grid.get_interval_id_left_neighbor(&IntervalId::new(3)),
            Some(IntervalId::new(2))
        );
        assert_eq!(
            grid.get_interval_id_right_neighbor(&IntervalId::new(3)),
            None
        );

        assert_eq!(
            grid.get_interval_id_neighbor(&IntervalId::new(2), &Side::Left),
            Some(IntervalId::new(1))
        );
        assert_eq!(
            grid.get_interval_id_neighbor(&IntervalId::new(2), &Side::Right),
            Some(IntervalId::new(3))
        );
    }

    #[test]
    fn test_interval_neighbors_circle_wraparound_uniform() {
        let grid = Grid1DUniform::new_periodic(
            IntervalLowerClosedUpperOpen::new(0.0, 1.0),
            NumIntervals::try_new(4).unwrap(),
        );

        let grid_index_spaces = grid.index_spaces();
        assert_eq!(
            grid_index_spaces.interval_index_space().topology(),
            &Topology1D::Circle
        );
        assert_eq!(
            grid_index_spaces.coord_index_space().topology(),
            &Topology1D::Circle
        );

        assert_eq!(
            grid.get_interval_id_left_neighbor(&IntervalId::new(0)),
            Some(IntervalId::new(3))
        );
        assert_eq!(
            grid.get_interval_id_right_neighbor(&IntervalId::new(3)),
            Some(IntervalId::new(0))
        );

        assert_eq!(
            grid.get_interval_id_neighbor(&IntervalId::new(0), &Side::Left),
            Some(IntervalId::new(3))
        );
        assert_eq!(
            grid.get_interval_id_neighbor(&IntervalId::new(3), &Side::Right),
            Some(IntervalId::new(0))
        );
    }

    #[test]
    fn test_interval_neighbors_circle_wraparound_non_uniform() {
        let coords = Coords1D::try_from(sorted_vec::partial::SortedSet::from_unsorted(vec![
            0.0, 0.2, 0.7, 1.0,
        ]))
        .unwrap();

        let grid =
            Grid1DNonUniform::<IntervalLowerClosedUpperOpen<f64>>::try_new_periodic_from_coords(
                coords,
            )
            .unwrap();

        let grid_index_spaces = grid.index_spaces();
        assert_eq!(
            grid_index_spaces.interval_index_space().topology(),
            &Topology1D::Circle
        );
        assert_eq!(
            grid_index_spaces.coord_index_space().topology(),
            &Topology1D::Circle
        );

        assert_eq!(
            grid.get_interval_id_left_neighbor(&IntervalId::new(0)),
            Some(IntervalId::new(2))
        );
        assert_eq!(
            grid.get_interval_id_right_neighbor(&IntervalId::new(2)),
            Some(IntervalId::new(0))
        );
    }
}