grid1d 0.5.1

#![deny(rustdoc::broken_intra_doc_links)]

//! Concrete 1D grid types with performance specialization.
//!
//! This module defines the three concrete grid types and their shared error enum:
//!
//! | Type | Point location | Memory | Best for |
//! |------|---------------|--------|----------|
//! | [`Grid1DUniform`] | **O(1)** analytical | O(n) | Equal spacing, FDM, spectral methods |
//! | [`Grid1DNonUniform`] | O(log n) binary search | O(n) | Adaptive spacing, boundary layers |
//! | [`Grid1D`] | O(1) or O(log n) | O(n) | Enum dispatch over the two variants above |
//!
//! The submodule [`traits`] defines the core grid traits ([`Grid1DTrait`],
//! [`HasDomain1D`], [`HasCoords1D`], etc.) that are re-exported at the crate root.
//! The submodule [`window`] provides [`Grid1DWindow`], a zero-copy borrowed view
//! over a contiguous slice of intervals within any [`Grid1DTrait`].
//!
//! # Examples
//!
//! ## Uniform grid
//!
//! ```rust
//! use grid1d::{*, intervals::*};
//! use try_create::TryNew;
//!
//! let grid = Grid1DUniform::new(
//!     IntervalClosed::new(0.0_f64, 1.0),
//!     NumIntervals::try_new(4).unwrap(),
//! );
//! assert_eq!(grid.num_intervals().as_ref(), &4);
//! ```
//!
//! ## Non-uniform grid from sorted coordinates
//!
//! ```rust
//! use grid1d::{*, intervals::*};
//! use sorted_vec::partial::SortedSet;
//!
//! let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
//!     SortedSet::from_unsorted(vec![0.0, 0.1, 0.5, 1.0]),
//! ).unwrap();
//! assert_eq!(grid.num_intervals().as_ref(), &3);
//! ```
//!
//! ## Uniform refinement
//!
//! ```rust
//! use grid1d::{*, intervals::*};
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(IntervalClosed::new(0.0_f64, 1.0), NumIntervals::try_new(4).unwrap());
//! let refined = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
//! assert_eq!(refined.refined_grid().num_intervals().as_ref(), &8);
//! ```

pub mod traits;
pub use traits::{
    FindIntervalIdOfPoint, Grid1DIntervalBuilder, Grid1DTrait, HasCoordIdRange, HasCoords1D,
    HasDomain1D, HasIntervalIdRange, HasIntervals,
};

pub mod window;
pub use window::{ErrorsGrid1DWindow, Grid1DWindow};

use crate::{
    CoordId, Coords1D, ErrorsCoords1D, ErrorsIntervalConstruction, ErrorsPositiveInt,
    Grid1DIndexSpaces, IntervalFinitePositiveLengthTrait, IntervalFromBounds, IntervalId,
    NumIntervals, PositiveNumPoints1D, SupportsCircularTopology, Topology1D,
    operations::refinement::{
        Grid1DNonUniformRefinement, Grid1DUniformRefinement, build_non_uniform_grid_refinement,
        refine_all_intervals_in_coords1d_same_subdivisions,
    },
};
use getset::Getters;
use num_valid::{
    RealScalar,
    core::errors::capture_backtrace,
    functions::{Reciprocal, Rounding},
    scalars::{PositiveRealScalar, RelativeTolerance},
};
use serde::{Deserialize, Serialize};
use sorted_vec::partial::SortedSet;
use std::backtrace::Backtrace;
use thiserror::Error;
use try_create::TryNew;

//------------------------------------------------------------------------------------------------------------------------------
///  Comprehensive error types for 1D grid construction and manipulation.
///
/// [`ErrorsGrid1D`] provides specific error information for all failure modes
/// in grid construction, enabling robust error handling and debugging.
///
/// ## Error Variants
///
/// ### Construction Errors
/// - [`ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints`]: Grid needs minimum 2 distinct coordinates
/// - [`ErrorsGrid1D::FromCoords1DConstructor`]: Error from underlying coordinate construction
/// - [`ErrorsGrid1D::FromIntervalConstructor`]: Error from underlying interval construction
///
/// # Examples
///
/// ## Handling Insufficient Points
///
/// ```rust
/// use grid1d::{
///     Coords1D, Grid1D, ErrorsGrid1D,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// let coords = SortedSet::from_unsorted(vec![0.0]); // Only one point provided
/// let err = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap_err();
/// assert!(matches!(err, ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { .. }), "{:?}", err);
/// ```
#[derive(Debug, Error)]
pub enum ErrorsGrid1D<Domain1D: IntervalFinitePositiveLengthTrait> {
    /// Indicates that fewer than two distinct points were provided, which is insufficient to define a valid grid.
    ///
    /// This error occurs when attempting to construct a [`Grid1D`] with fewer than two points, as at least two points
    /// are required to define a partition of a 1D interval.
    #[error("The number of distinct points required to define an interval partition must be >= 2!")]
    RequiredAtLeastTwoDistinctPoints {
        /// The actual number of points that were provided.
        num_points_provided: PositiveNumPoints1D,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },

    /// Error propagated from the underlying [`Coords1D`] construction.
    #[error("Error from Coords1D construction")]
    FromCoords1DConstructor {
        /// The underlying error from [`Coords1D`] construction.
        #[from]
        source: ErrorsCoords1D,
    },

    /// Error propagated from the construction of the underlying interval.
    #[error("Error from the construction of the interval")]
    FromIntervalConstructor {
        /// The underlying error from interval construction.
        #[from]
        source: ErrorsIntervalConstruction<Domain1D::RealType>,
    },

    #[error(
        "The first point of the grid ({first_point}) is not equal to the lower bound of the interval representing the domain {domain:?}!"
    )]
    /// Indicates that the first coordinate does not match the domain's lower bound.
    FirstPointIsNotDomainLowerBound {
        /// The actual first coordinate found in the grid.
        first_point: Domain1D::RealType,
        /// The domain interval whose lower bound was expected.
        domain: Domain1D,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },

    #[error(
        "The last point of the grid ({last_point}) is not equal to the upper bound of the interval representing the domain {domain:?}!"
    )]
    /// Indicates that the last coordinate does not match the domain's upper bound.
    LastPointIsNotDomainUpperBound {
        /// The actual last coordinate found in the grid.
        last_point: Domain1D::RealType,
        /// The domain interval whose upper bound was expected.
        domain: Domain1D,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },

    #[error(transparent)]
    /// Wraps a [`ErrorsPositiveInt`] from an invalid interval count.
    NumIntervals(#[from] ErrorsPositiveInt),
}
//------------------------------------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------
pub(crate) fn validate_coords_in_domain_closure<Domain1D: IntervalFinitePositiveLengthTrait>(
    coords: &Coords1D<Domain1D::RealType>,
    domain: &Domain1D,
) -> Result<NumIntervals, ErrorsGrid1D<Domain1D>> {
    let lower_bound = domain.lower_bound_value();
    let upper_bound = domain.upper_bound_value();

    let first_point = coords.first();
    if first_point != lower_bound {
        return Err(ErrorsGrid1D::FirstPointIsNotDomainLowerBound {
            first_point: first_point.clone(),
            domain: domain.clone(),
            backtrace: capture_backtrace(),
        });
    }

    let last_point = coords.last();
    if last_point != upper_bound {
        return Err(ErrorsGrid1D::LastPointIsNotDomainUpperBound {
            last_point: last_point.clone(),
            domain: domain.clone(),
            backtrace: capture_backtrace(),
        });
    }

    let num_intervals = NumIntervals::try_new(coords.num_points().as_ref() - 1)?;

    Ok(num_intervals)
}
//------------------------------------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------
///  A general-purpose, non-uniform partition of a finite interval with positive length.
///
/// The [`Grid1DNonUniform`] struct represents a partition of a bounded interval into
/// adjacent, non-overlapping sub-intervals where **the spacing between points can be arbitrary**.
/// This provides maximum flexibility for creating adaptive grids that can resolve complex
/// features, such as boundary layers, shock fronts, or regions with high solution gradients.
///
/// Unlike [`Grid1DUniform`] which enforces equal spacing, [`Grid1DNonUniform`] is defined by an
/// explicit set of coordinates, allowing for fine-grained control over the mesh resolution.
///
/// ## Mathematical Definition
///
/// Given a domain `D` and a set of `n+1` unique, sorted coordinates `{p₀, p₁, ..., pₙ}`,
/// a non-uniform grid is formed if and only if:
/// 1.  `p₀` equals the lower bound of `D`.
/// 2.  `pₙ` equals the upper bound of `D`.
///
/// This creates `n` sub-intervals `{I₀, I₁, ..., I_{n-1}}` where `Iᵢ` is the interval
/// between `pᵢ` and `pᵢ₊₁`. The length of each interval, `δᵢ = pᵢ₊₁ - pᵢ`, can vary.
///
/// ## Core Mathematical Properties
///
/// - **Arbitrary Spacing**: Interval lengths `δᵢ` can be different, allowing for adaptive resolution.
/// - **Boundary Matching**: The first and last grid points are guaranteed to match the domain's lower and upper bounds exactly. This is a critical validation step.
/// - **Complete Coverage**: `D = I₀ ∪ I₁ ∪ ... ∪ I_{n-1}` (the union of sub-intervals perfectly reconstructs the domain).
/// - **Non-Overlapping**: `Iᵢ ∩ Iⱼ = ∅` for `i ≠ j` (sub-intervals are disjoint).
/// - **Inherited Properties**: Inherits uniqueness, sorting, and non-emptiness from [`Coords1D`].
///
/// ## Performance Characteristics vs. Uniform Grids
///
/// | Operation | [`Grid1DNonUniform`] | [`Grid1DUniform`] | Advantage |
/// |-----------|----------------------|-------------------|-----------|
/// | **Construction** | O(n) + validation | O(n) | Equal |
/// | **Memory Storage** | O(n) | O(n) | Equal (both store coordinates) |
/// | **Point Location** | O(log n) binary search | O(1) analytical | ✅ **Uniform much faster** |
/// | **Interval Access** | O(1) | O(1) | Equal |
/// | **Coordinate Access** | O(1) | O(1) | Equal (pre-computed) |
/// | **Length Queries** | O(1) computed | O(1) cached | ✅ **Uniform simpler** |
///
/// ## Design Philosophy
///
/// ### Correctness Through Strict Validation
/// The primary design principle is **mathematical correctness**. The `try_new` constructor is strict and performs critical runtime checks to ensure the grid is a valid partition *of the specified domain*. This prevents common numerical errors, such as:
/// - Grids that do not fully cover the intended domain.
/// - Off-by-one errors at boundaries.
/// - Inconsistencies between the grid's span and the domain's length.
///
/// ### Flexibility and Generality
/// By accepting any valid `Coords1D`, this struct can represent any one-dimensional partition, making it the general-purpose backbone of the `Grid1D` system. It is suitable for:
/// - Manually defined meshes.
/// - Programmatically generated adaptive meshes.
/// - Results of grid union operations.
///
/// ### Integration with numerical Ecosystem
/// - **Unified Interface**: Implements the `Grid1DTrait` trait, providing the same API as `Grid1DUniform`.
/// - **Grid Enum Integration**: Serves as the general-purpose variant in the `Grid1D` enum.
/// - **Domain-Awareness**: Correctly handles open, closed, and semi-open domains through the `Grid1DIntervalBuilder` trait.
///
/// ## Construction and Usage Examples
///
/// ### Basic Construction
/// ```rust
/// use grid1d::{
///     Grid1DNonUniform, HasCoords1D, Grid1DTrait, HasIntervalIdRange,
///     intervals::*, HasDomain1D,
///     scalars::NumIntervals,
/// };
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// // Define a non-uniform set of coordinates
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 0.1, 0.5, 0.9, 1.0])
/// ).unwrap();
///
/// // Create the non-uniform grid
/// let grid = Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(coords).unwrap();
///
/// // Verify properties
/// assert_eq!(grid.num_intervals().as_ref(), &4);
/// assert_eq!(grid.coords().deref(), &[0.0, 0.1, 0.5, 0.9, 1.0]);
/// ```
///
/// ### Adaptive Mesh for a Boundary Layer
/// ```rust
/// use grid1d::{
///     Coords1D, Grid1DNonUniform, Grid1DTrait, HasCoords1D,
///     intervals::*,
///     scalars::IntervalId,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // Simulate a boundary layer by clustering points near x=0
/// let mut points = vec![0.0, 1.0];
/// for i in 1..=5 {
///     points.push(1.0 / (2.0_f64.powi(i)));
/// }
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(points)).unwrap();
/// let grid = Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(coords).unwrap();
///
/// println!("Adaptive grid points: {:?}", grid.coords());
///
/// // Interval lengths decrease towards the boundary
/// let id_0 = IntervalId::new(0);
/// let id_1 = IntervalId::new(1);
/// more_asserts::assert_lt!((grid.interval_length(&id_0).into_inner() - grid.interval_length(&id_1).as_ref()).abs(), 1e-15);
/// ```
///
/// ### Error Handling During Construction
/// ```rust
/// use grid1d::{
///     Grid1DNonUniform, ErrorsGrid1D, Coords1D,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // Not enough points
/// let coords_too_few = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0])).unwrap();
/// let result = Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(coords_too_few);
/// assert!(matches!(result, Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { .. })));
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`Grid1DNonUniform`]
/// - **Adaptive Mesh Refinement (AMR)**: When you need higher resolution in specific areas (e.g., near boundaries, shocks, or high-gradient regions).
/// - **Complex Geometries**: When mapping from a simple computational domain to a complex physical domain results in non-uniform spacing.
/// - **Multi-Physics Coupling**: When combining grids from different simulations, the resulting union is typically non-uniform.
/// - **Error-Driven Adaptation**: When grid points are placed dynamically to minimize numerical error.
///
/// ### When to Prefer [`Grid1DUniform`]
/// - For problems where the solution is smooth and does not have sharp features.
/// - When using standard finite difference or Fourier spectral methods that rely on equal spacing.
/// - When point location performance is the absolute top priority (O(1) vs. O(log n)).
///
/// ## Mathematical Properties and Guarantees
///
/// ### Boundary Matching (Runtime Validated)
/// ```text
/// coords.first() == domain.lower_bound_value()
/// coords.last()  == domain.upper_bound_value()
/// ```
/// This is the core invariant enforced by `try_new`.
///
/// Note: these are equality checks on boundary *values*. Containment of those
/// boundary values still follows the domain boundary semantics (open/closed).
/// For example, with an open domain `(a, b)`, `a` and `b` may appear as endpoint
/// coordinate values while remaining outside the domain as points.
///
/// ### Partition Completeness
/// The union of all sub-intervals perfectly reconstructs the domain.
/// ```text
/// D = I₀ ∪ I₁ ∪ ... ∪ Iₙ₋₁
/// ```
///
/// ### Point Location Uniqueness
/// Every point `x` in the domain belongs to exactly one sub-interval.
///
/// ## See Also
///
/// - [`Grid1DUniform`]: The specialized, high-performance alternative for equally spaced grids.
/// - [`Grid1D`]: The enum that unifies both uniform and non-uniform grids under a single API.
/// - [`Grid1DTrait`]: The core trait providing partition operations.
/// - [`Coords1D`]: The underlying coordinate container.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
#[serde(bound(deserialize = "Domain1D: for<'d> serde::Deserialize<'d>, \
                  Domain1D::RealType: for<'d> serde::Deserialize<'d>"))]
pub struct Grid1DNonUniform<Domain1D: IntervalFinitePositiveLengthTrait> {
    /// Domain covered by the grid.
    pub(crate) domain: Domain1D,

    /// The index spaces associated with the grid, which define the topology and interval indexing.
    pub(crate) index_spaces: Grid1DIndexSpaces,

    /// Strictly sorted, duplicate-free coordinates of partition points.
    pub(crate) coords: Coords1D<Domain1D::RealType>,
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> Grid1DNonUniform<Domain1D> {
    /// Creates a new non-uniform grid from an explicit domain and pre-validated coordinates.
    ///
    /// Validates that the coordinates are consistent with the domain (first and last
    /// coordinate must match the domain bounds). Prefer [`Grid1D::try_from_sorted`]
    /// or [`Grid1D::try_from_coords`] when the domain can be inferred from the coordinates.
    pub fn try_new(
        domain: Domain1D,
        coords: Coords1D<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>> {
        let num_intervals = validate_coords_in_domain_closure(&coords, &domain)?;

        let index_spaces = Grid1DIndexSpaces::new_non_periodic(num_intervals);

        Ok(Self {
            domain,
            coords,
            index_spaces,
        })
    }
}

impl<Domain1D: IntervalFromBounds> Grid1DNonUniform<Domain1D> {
    fn create_domain_from_coords(
        coords: &Coords1D<Domain1D::RealType>,
    ) -> Result<Domain1D, ErrorsGrid1D<Domain1D>>
//    where
//        Domain1D: IntervalFromBounds,
    {
        let n_pts = coords.num_points();
        if n_pts.as_ref() < &2 {
            Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints {
                num_points_provided: n_pts,
                backtrace: capture_backtrace(),
            })
        } else {
            Ok(Domain1D::try_new(
                coords.first().clone(),
                coords.last().clone(),
            )?)
        }
    }

    /// Attempts to create an [`Grid1DNonUniform`] from a set of coordinates modeled by a [`Coords1D`] object.
    ///
    /// # Parameters
    ///
    /// - `coords`: A [`Coords1D`] object containing the sorted coordinates of the points.
    ///
    /// # Errors
    ///
    /// - Returns [`ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints`] if fewer than 2 points are provided.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Coords1D, Grid1DNonUniform, HasDomain1D,
    ///     intervals::*,
    /// };
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 1.0])).unwrap();
    /// let partition = Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(coords).unwrap();
    ///
    /// let domain = partition.domain();
    /// assert_eq!(domain.lower_bound_value(), &0.0);
    /// assert_eq!(domain.upper_bound_value(), &1.0);
    /// ```
    pub fn try_new_from_coords(
        coords: Coords1D<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>> {
        let domain = Self::create_domain_from_coords(&coords)?;

        let num_intervals = NumIntervals::try_new(coords.num_points().as_ref() - 1)
            .expect("At least 2 points are required to form a valid grid!");
        let index_spaces = Grid1DIndexSpaces::new_non_periodic(num_intervals);

        Ok(Self {
            domain,
            coords,
            index_spaces,
        })
    }
}

impl<Domain1D> Grid1DNonUniform<Domain1D>
where
    Domain1D: IntervalFinitePositiveLengthTrait + SupportsCircularTopology,
{
    /// Attempts to create a periodic non-uniform grid from sorted coordinates.
    ///
    /// The resulting grid uses [`Topology1D::Circle`], so the first and last
    /// intervals are considered neighbors.
    ///
    /// # Errors
    ///
    /// Returns [`ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints`] if fewer than
    /// two distinct points are provided.
    pub fn try_new_periodic(
        domain: Domain1D,
        coords: Coords1D<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>> {
        let num_intervals = validate_coords_in_domain_closure(&coords, &domain)?;

        let index_spaces = if domain.is_lower_bound_open() {
            Grid1DIndexSpaces::new_periodic_lower_open_interval(num_intervals)
        } else {
            Grid1DIndexSpaces::new_periodic_upper_open_interval(num_intervals)
        };

        Ok(Self {
            domain,
            coords,
            index_spaces,
        })
    }
}

impl<Domain1D> Grid1DNonUniform<Domain1D>
where
    Domain1D: IntervalFromBounds + SupportsCircularTopology,
{
    /// Attempts to create a periodic non-uniform grid from sorted coordinates.
    ///
    /// The resulting grid uses [`Topology1D::Circle`], so the first and last
    /// intervals are considered neighbors.
    ///
    /// # Errors
    ///
    /// Returns [`ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints`] if fewer than
    /// two distinct points are provided.
    pub fn try_new_periodic_from_coords(
        coords: Coords1D<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>> {
        let domain = Self::create_domain_from_coords(&coords)?;

        let num_intervals = NumIntervals::try_new(coords.num_points().as_ref() - 1)
            .expect("At least 2 points are required to form a valid grid!");

        let index_spaces = if domain.is_lower_bound_open() {
            Grid1DIndexSpaces::new_periodic_lower_open_interval(num_intervals)
        } else {
            Grid1DIndexSpaces::new_periodic_upper_open_interval(num_intervals)
        };

        Ok(Self {
            domain,
            coords,
            index_spaces,
        })
    }
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> HasDomain1D for Grid1DNonUniform<Domain1D> {
    type Domain1D = Domain1D;

    #[inline(always)]
    fn domain(&self) -> &Domain1D {
        &self.domain
    }
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> HasCoords1D for Grid1DNonUniform<Domain1D> {
    type CoordType = Domain1D::RealType;

    #[inline(always)]
    fn coords(&self) -> &Coords1D<Self::CoordType> {
        &self.coords
    }
}

impl<Domain1D: Grid1DIntervalBuilder> FindIntervalIdOfPoint for Grid1DNonUniform<Domain1D> {
    type Point1DType = Domain1D::RealType;

    fn find_interval_id_of_point(&self, x: &Self::Point1DType) -> Option<IntervalId> {
        if !self.domain().contains_point(x) {
            return None;
        }

        // Here we use the domain's lower bound closure property to determine the assignment of interior partition points.
        // This is because the domain type dictates the sub-interval construction rules, which in turn determine how boundary points are assigned to intervals.
        // For example:
        // - For a closed domain [a, b], all sub-intervals are left-closed, so interior points belong to the right interval.
        // - For a half-open domain (a, b], all sub-intervals are right-closed, so interior points belong to the left interval.
        // - For a half-open domain [a, b), the first interval is left-closed and the last interval is open, so interior points belong to the right interval.
        // - For a open domain (a, b), all sub-intervals are open, so interior points do not belong to any interval, but since we are checking for exact matches with grid points, we can still assign them based on the closure of the lower bound.
        // This logic ensures that the point location is consistent with the interval construction rules defined by the domain type, while still leveraging the uniform spacing for O(1) access.
        // The key point is that the domain's lower bound closure property directly informs how we assign interior partition points to intervals, which is crucial for maintaining the correct semantics of the partitioning.
        if self.domain().is_lower_bound_closed() {
            // Sub-intervals are left-closed [pₖ, pₖ₊₁): boundary point pₖ belongs to interval k.
            // find_floor_index returns the largest index j where coords[j] ≤ x.
            // Clamped to num_intervals-1 to handle the closed upper boundary of [a,b] domains.
            let floor = self
                .coords()
                .find_floor_index(x)
                .expect("Point is in domain but has no floor coordinate");
            let final_id = floor.min(*self.num_intervals().as_ref() - 1);
            Some(IntervalId::new(final_id))
        } else {
            // Sub-intervals are right-closed (pₖ, pₖ₊₁]: boundary point pₖ belongs to interval k-1.
            // partition_point(c < x) returns the first index where coord ≥ x.
            let idx = self.coords().find_insertion_index(x);
            if idx == 0 {
                Some(IntervalId::new(0))
            } else {
                Some(IntervalId::new(idx - 1))
            }
        }
    }
}

impl<Domain1D> HasIntervalIdRange for Grid1DNonUniform<Domain1D>
where
    Domain1D: Grid1DIntervalBuilder,
{
    #[inline(always)]
    fn first_interval_id(&self) -> IntervalId {
        self.index_spaces.interval_index_space().first_interval_id()
    }

    #[inline(always)]
    fn last_interval_id(&self) -> IntervalId {
        self.index_spaces.interval_index_space().last_interval_id()
    }
}

impl<Domain1D> HasCoordIdRange for Grid1DNonUniform<Domain1D>
where
    Domain1D: Grid1DIntervalBuilder,
{
    #[inline(always)]
    fn first_coord_id(&self) -> CoordId {
        self.index_spaces.coord_index_space().first_coord_id()
    }

    #[inline(always)]
    fn last_coord_id(&self) -> CoordId {
        self.index_spaces.coord_index_space().last_coord_id()
    }
}

impl<Domain1D> Grid1DTrait for Grid1DNonUniform<Domain1D>
where
    Domain1D: Grid1DIntervalBuilder,
{
    #[inline(always)]
    fn index_spaces(&self) -> &Grid1DIndexSpaces {
        &self.index_spaces
    }

    type UniformlyRefinedGrid1DType = Self;

    fn refine_uniform(
        self,
        num_extra_points_each_interval: &PositiveNumPoints1D,
    ) -> Grid1DUniformRefinement<Self> {
        let (
            refined_coords,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        ) = refine_all_intervals_in_coords1d_same_subdivisions(
            self.coords(),
            num_extra_points_each_interval,
        );

        let domain = self.domain().clone();
        let num_refined_intervals = validate_coords_in_domain_closure(&refined_coords, &domain)
            .expect("Refined coordinates do not form a valid grid!");

        let topology = self.index_spaces.topology();
        let refined_index_spaces = match topology {
            Topology1D::RealLine => Grid1DIndexSpaces::new_non_periodic(num_refined_intervals),
            Topology1D::Circle => {
                if self.domain().is_lower_bound_open() {
                    Grid1DIndexSpaces::new_periodic_lower_open_interval(num_refined_intervals)
                } else {
                    Grid1DIndexSpaces::new_periodic_upper_open_interval(num_refined_intervals)
                }
            }
        };

        let refined_grid = Self {
            domain,
            coords: refined_coords,
            index_spaces: refined_index_spaces,
        };

        Grid1DUniformRefinement::new(
            refined_grid,
            self,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        )
    }

    fn refine(
        self,
        intervals_to_refine: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
    ) -> Grid1DNonUniformRefinement<Self> {
        build_non_uniform_grid_refinement(self, intervals_to_refine)
    }
}
//-------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------
///  A uniform partition of a finite interval with positive length into equal-sized sub-intervals.
///
/// The [`Grid1DUniform`] struct represents a partition of a bounded interval into
/// adjacent, non-overlapping sub-intervals where **all sub-intervals have exactly the same length**.
/// This provides optimal performance and memory efficiency for regular grids commonly used in
/// finite difference methods, spectral computations, and uniform discretizations.
///
/// Unlike [`Grid1DNonUniform`] which allows arbitrary point distributions, [`Grid1DUniform`]
/// enforces equal spacing while providing specialized optimizations for uniform access patterns.
///
/// ## Mathematical Definition
///
/// Given a domain `D` with finite positive length and a number of intervals `n`, the uniform
/// partition creates `n+1` equally spaced points:
/// ```text
/// p_k = lower_bound(D) + k * δ  for k = 0, 1, ..., n
/// ```
/// where `δ = length(D) / n` is the uniform spacing.
///
/// This creates `n` sub-intervals `{I₀, I₁, ..., I_{n-1}}` where each interval has length `δ`.
///
/// ## Core Mathematical Properties
///
/// - **Uniform Spacing**: All sub-intervals have identical length `δ = domain_length / num_intervals`
/// - **Complete Coverage**: `D = I₀ ∪ I₁ ∪ ... ∪ I_{n-1}` (complete domain coverage)
/// - **Non-Overlapping**: `Iᵢ ∩ Iⱼ = ∅` for `i ≠ j` (disjoint sub-intervals)
/// - **Boundary Matching**: First point equals domain lower bound, last point equals upper bound
///   (value matching only). For open or semi-open domains, those endpoint values
///   are not necessarily contained points of the domain.
/// - **Regularity**: Perfect geometric regularity enables analytical optimizations
/// - **Predictable Access**: Any coordinate can be computed as `p_k = lower + k * δ`
///
/// ## Performance Characteristics vs. Non-Uniform Grids
///
/// | Operation | [`Grid1DUniform`] | [`Grid1DNonUniform`] | Advantage |
/// |-----------|-------------------|----------------------|-----------|
/// | **Construction** | O(n) | O(n) | Equal |
/// | **Memory Storage** | O(n) | O(n) | Equal (both cache coordinates) |
/// | **Point Location** | O(1) analytical | O(log n) binary search | ✅ **Uniform much faster** |
/// | **Interval Access** | O(1) | O(1) | Equal |
/// | **Coordinate Access** | O(1) | O(1) | Equal (both pre-computed) |
/// | **Length Queries** | O(1) cached | O(1) computed | ✅ **Uniform simpler** |
/// | **Cache Efficiency** | Excellent | Good | ✅ **Uniform better** |
///
/// ## Design Philosophy
///
/// ### Performance Through Uniformity
/// - **Analytical Formulas**: Many operations can use closed-form expressions
/// - **Predictable Memory Patterns**: Enables excellent CPU cache utilization
/// - **SIMD Optimization**: Regular patterns are ideal for vectorized operations
/// - **Minimal Branching**: Uniform structure reduces conditional logic
///
/// ### Memory vs. Computation Trade-offs
/// **Current Design (O(n) memory)**:
/// - ✅ **Compatibility**: Same interface as [`Grid1DNonUniform`]
/// - ✅ **Zero computation overhead**: Pre-computed coordinates
/// - ✅ **Simple access patterns**: Direct array indexing
/// - ❌ **Memory overhead**: Stores redundant coordinate data
///
/// **Alternative O(1) Design** (not implemented):
/// - ✅ **Memory efficient**: Only stores `(domain, num_intervals, δ)`
/// - ✅ **Analytical point location**: No binary search needed
/// - ❌ **Computation overhead**: Must compute coordinates on-demand
/// - ❌ **Interface complexity**: Different from non-uniform grids
///
/// ### Integration with other numerical Ecosystem
/// - **Unified Interface**: Implements [`Grid1DTrait`] trait like non-uniform grids
/// - **Trait Compatibility**: Same [`HasDomain1D`] and [`HasCoords1D`] interfaces
/// - **Grid Enum Integration**: Seamlessly works within [`Grid1D`] enum
/// - **Optimization Opportunities**: Specialized methods for uniform-specific operations
///
/// ## Construction and Usage Examples
///
/// ### Basic Construction
/// ```rust
/// use grid1d::{
///     Grid1DUniform, HasCoords1D, HasDomain1D, Grid1DTrait, HasIntervalIdRange,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// use std::ops::Deref;
/// use try_create::TryNew;
///
/// // Create uniform grid with 4 intervals
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let num_intervals = NumIntervals::try_new(4).unwrap();
/// let grid = Grid1DUniform::new(domain.clone(), num_intervals);
///
/// // Verify properties
/// assert_eq!(grid.domain(), &domain);
/// assert_eq!(grid.num_intervals().as_ref(), &4);
/// assert_eq!(grid.coords().deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
/// assert_eq!(grid.delta_points().as_ref(), &0.25);
/// ```
///
/// ### High-Resolution Grids
/// ```rust
/// use grid1d::{*, intervals::*};
/// use try_create::TryNew;
///
/// // Create high-resolution grid for numerical simulations
/// let domain = IntervalClosed::new(0.0, 10.0);
/// let num_intervals = NumIntervals::try_new(10000).unwrap();
/// let grid = Grid1DUniform::new(domain, num_intervals);
///
/// // Properties
/// assert_eq!(grid.coords().len(), 10001); // n+1 points
/// more_asserts::assert_lt!((*grid.delta_points().as_ref() - 0.001).abs(), 1.0e-15); // Δx = 10.0/10000
///
/// // Efficient operations on large grids
/// let interval_id = grid.find_interval_id_of_point(&0.).unwrap();
/// assert_eq!(interval_id, IntervalId::new(0));
///
/// let interval_id = grid.find_interval_id_of_point(&0.0009).unwrap();
/// assert_eq!(interval_id, IntervalId::new(0));
///
/// let interval_id = grid.find_interval_id_of_point(&0.001).unwrap();
/// assert_eq!(interval_id, IntervalId::new(1)); // p₁ = 0.001 belongs to interval 1 (left-closed)
///
/// let interval_id = grid.find_interval_id_of_point(&0.0015).unwrap();
/// assert_eq!(interval_id, IntervalId::new(1));
///
/// let interval_id = grid.find_interval_id_of_point(&0.002).unwrap();
/// assert_eq!(interval_id, IntervalId::new(2)); // p₂ = 0.002 belongs to interval 2 (left-closed)
///
/// let interval_id = grid.find_interval_id_of_point(&2.537).unwrap();
/// assert_eq!(interval_id, IntervalId::new(2537)); // p₂₅₃₇ = 2.537 belongs to interval 2537 (left-closed)
///
/// let interval_id = grid.find_interval_id_of_point(&9.999).unwrap();
/// assert_eq!(interval_id, IntervalId::new(9999)); // p₉₉₉₉ = 9.999 belongs to interval 9999 (left-closed)
///
/// let interval_id = grid.find_interval_id_of_point(&9.9991).unwrap();
/// assert_eq!(interval_id, IntervalId::new(9999));
///
/// let interval_id = grid.find_interval_id_of_point(&10.000).unwrap();
/// assert_eq!(interval_id, IntervalId::new(9999));
///
/// ```
///
/// ### Different Domain Types
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// // Closed interval [0, 2π] for Fourier analysis
/// let fourier_domain = IntervalClosed::new(0.0, 2.0 * std::f64::consts::PI);
/// let fourier_grid = Grid1DUniform::new(fourier_domain, NumIntervals::try_new(256).unwrap());
///
/// // Semi-open interval [0, 1) for periodic boundaries
/// let periodic_domain = IntervalLowerClosedUpperOpen::new(0.0, 1.0);
/// let periodic_grid = Grid1DUniform::new(periodic_domain, NumIntervals::try_new(100).unwrap());
///
/// // Different boundary semantics in partition
/// let closed_interval_0 = fourier_grid.interval(&IntervalId::new(0));
/// let periodic_interval_0 = periodic_grid.interval(&IntervalId::new(0));
///
/// // Both have uniform spacing but different boundary conditions
/// let diff_delta_fourier_grid = (fourier_grid.delta_points().clone().into_inner() - 2.0 * std::f64::consts::PI / 256.0).abs();
/// more_asserts::assert_lt!(diff_delta_fourier_grid, 1.0e-15);
///
/// let diff_delta_periodic_grid = (periodic_grid.delta_points().clone().into_inner() - &0.01).abs();
/// more_asserts::assert_lt!(diff_delta_periodic_grid, 1.0e-15);
/// ```
///
/// ### Working with Different Scalar Types
/// ```rust
/// use grid1d::{
///     Grid1DUniform, HasCoords1D,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Create grid with validated scalar type
/// let domain = IntervalClosed::new(
///     Real::try_new(-1.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let grid = Grid1DUniform::new(domain, NumIntervals::try_new(8).unwrap());
///
/// // All operations work with the validated type
/// assert_eq!(grid.coords().num_points().as_ref(), &9);
/// assert_eq!(grid.delta_points().as_ref(), &Real::try_new(0.25).unwrap());
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Finite Difference Method Integration
/// ```rust
/// use grid1d::{
///     Grid1DUniform, HasCoords1D, Grid1DTrait,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// use try_create::TryNew;
///
/// fn finite_difference_setup(
///     domain: IntervalClosed<f64>,
///     resolution: usize
/// ) -> (Grid1DUniform<IntervalClosed<f64>>, Vec<f64>) {
///     let grid = Grid1DUniform::new(domain, NumIntervals::try_new(resolution).unwrap());
///     
///     // Initialize solution vector
///     let solution = vec![0.0; grid.coords().len()];
///     
///     (grid, solution)
/// }
///
/// // Usage for heat equation on [0, 1]
/// let (grid, mut u) = finite_difference_setup(IntervalClosed::new(0.0, 1.0), 100);
/// let dx = *grid.delta_points().as_ref();
/// let dt = 0.5 * dx * dx; // Stability condition
///
/// // Apply finite difference stencil
/// for i in 1..u.len()-1 {
///     let d2u_dx2 = (u[i-1] - 2.0*u[i] + u[i+1]) / (dx * dx);
///     // Time stepping logic here...
/// }
/// ```
///
/// ### Spectral Method Grid Setup
/// ```rust
/// use grid1d::{
///     Grid1DUniform, HasCoords1D,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use try_create::TryNew;
///
/// fn spectral_grid_setup(n_modes: usize) -> Grid1DUniform<IntervalLowerClosedUpperOpen<f64>> {
///     // Periodic domain [0, 2π) for Fourier spectral methods
///     let domain = IntervalLowerClosedUpperOpen::new(0.0, 2.0 * std::f64::consts::PI);
///     Grid1DUniform::new(domain, NumIntervals::try_new(n_modes).unwrap())
/// }
///
/// let spectral_grid = spectral_grid_setup(256);
/// let coords = spectral_grid.coords();
///
/// // Create FFT-ready data
/// let mut function_values: Vec<f64> = coords.iter()
///     .map(|&x| (x.sin() + 0.5 * (2.0 * x).cos()))
///     .collect();
///
/// // Grid spacing is perfect for FFT
/// let dx = *spectral_grid.delta_points().as_ref();
/// assert!((dx - 2.0 * std::f64::consts::PI / 256.0).abs() < 1e-15);
/// ```
///
/// ### Performance-Critical Point Location
/// ```rust
/// use grid1d::{*, intervals::*};
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalClosed::new(0.0, 100.0),
///     NumIntervals::try_new(10000).unwrap()
/// );
///
/// // Batch point location - O(1) per point for uniform grids
/// let test_points: Vec<f64> = (0..1000).map(|i| i as f64 * 0.1).collect();
/// let start = std::time::Instant::now();
///
/// let interval_ids: Vec<IntervalId> = test_points.iter()
///     .map(|&point| grid.find_interval_id_of_point(&point).unwrap())
///     .collect();
///
/// let elapsed = start.elapsed();
/// println!("Located {} points in {:.2}μs ({:.2}ns per point)",
///          test_points.len(), elapsed.as_micros(), elapsed.as_nanos() as f64 / test_points.len() as f64);
///
/// // For uniform grids, all intervals have the same length
/// for &id in &interval_ids[0..10] {
///     let length_diff = (grid.interval_length(&id).into_inner()  - grid.delta_points().as_ref()).abs();
///     more_asserts::assert_lt!(length_diff, 1.0e-15);
/// }
/// ```
///
/// ### Grid Quality and Analysis
/// ```rust
/// use grid1d::{
///     Grid1DUniform, Grid1DTrait, HasIntervalIdRange, HasDomain1D,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalClosed::new(-1.0, 1.0),
///     NumIntervals::try_new(1000).unwrap()
/// );
///
/// // Perfect uniformity properties
/// let min_length = grid.min_interval_length();
/// let max_length = grid.max_interval_length();
/// let uniformity_ratio = grid.uniformity_ratio();
///
/// // For uniform grids, these are always equal
/// assert_eq!(min_length, max_length);
/// assert_eq!(min_length, *grid.delta_points());
/// assert_eq!(uniformity_ratio.as_ref(), &1.0); // Perfect uniformity
///
/// println!("Grid Analysis:");
/// println!("  Domain: [{}, {}]",
///          grid.domain().lower_bound_value(), grid.domain().upper_bound_value());
/// println!("  Intervals: {}", grid.num_intervals().as_ref());
/// println!("  Spacing: {:.6}", grid.delta_points().as_ref());
/// println!("  Uniformity ratio: {:.1} (perfect)", uniformity_ratio.as_ref());
/// ```
///
/// ### Integration with Grid1D Enum
/// ```rust
/// use grid1d::{*, intervals::*};
/// use try_create::TryNew;
///
/// // Grid1DUniform is used internally by Grid1D for uniform cases
/// let domain = IntervalClosed::new(0.0, 10.0);
/// let uniform_grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(100).unwrap());
///
/// // Grid1D automatically uses Grid1DUniform for uniform spacing
/// match &uniform_grid {
///     Grid1D::Uniform(uniform_partition) => {
///         println!("Using Grid1DUniform with spacing: {:.3}",
///                  uniform_partition.delta_points().as_ref());
///         
///         // Access uniform-specific properties
///         assert_eq!(uniform_partition.delta_points().as_ref(), &0.1);
///         
///         // All standard Grid1DTrait methods work
///         assert_eq!(uniform_partition.uniformity_ratio().as_ref(), &1.0);
///     }
///     Grid1D::NonUniform(_) => {
///         panic!("Expected uniform grid");
///     }
/// }
///
/// // Unified interface works regardless of internal type
/// let point = 5.37;
/// let interval_id = uniform_grid.find_interval_id_of_point(&point).unwrap();
/// let length = uniform_grid.interval_length(&interval_id);
/// println!("Point {} in interval {} with length {:.3}",
///          point, interval_id.as_ref(), length.as_ref());
/// ```
///
/// ### Memory Layout and Performance Optimization
/// ```rust
/// use grid1d::{
///     Grid1DUniform, HasCoords1D,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use std::mem;
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalClosed::new(0.0, 1000.0),
///     NumIntervals::try_new(10000).unwrap()
/// );
///
/// // Memory analysis
/// let coords_size = mem::size_of_val(grid.coords());
/// let delta_size = mem::size_of_val(grid.delta_points());
/// let total_size = mem::size_of_val(&grid);
///
/// println!("Memory breakdown:");
/// println!("  Coordinates: {} bytes", coords_size);
/// println!("  Delta storage: {} bytes", delta_size);
/// println!("  Total grid: {} bytes", total_size);
/// println!("  Overhead: {} bytes", total_size - coords_size);
///
/// // Cache-friendly access patterns
/// let coords = grid.coords().as_ref();
/// let mut sum = 0.0;
///
/// // Sequential access is very cache-friendly for uniform grids
/// for &coord in coords.iter() {
///     sum += coord; // Perfect memory locality
/// }
///
/// let last_coord = coords.last().unwrap();
/// let analytical_sum = coords.len() as f64 * (coords[0] + last_coord) / 2.0;
/// assert!((sum - analytical_sum).abs() < 1e-10);
/// ```
///
/// ## Domain-Specific Behavior
///
/// The [`Grid1DUniform`] works with different domain types, each with specific
/// sub-interval construction rules:
///
/// ### Closed Domain `[a, b]`
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalClosed::new(0.0, 3.0),
///     NumIntervals::try_new(3).unwrap()
/// );
///
/// // First interval: [0, 1) (closed on lower bound, open on upper bound)
/// let interval_0 = grid.interval(&IntervalId::new(0));
/// assert!(interval_0.contains_point(&0.0)); // ✓ includes lower bound
/// assert!(!interval_0.contains_point(&1.0)); // ✗ excludes upper bound
///
/// // Middle interval: [1, 2) (closed on lower bound, open on upper bound for middle intervals)
/// let interval_1 = grid.interval(&IntervalId::new(1));
/// assert!(interval_1.contains_point(&1.0)); // ✓ includes left bound
/// assert!(!interval_1.contains_point(&2.0));  // ✗ excludes right bound
///
/// // Last interval: [2, 3] (closed on both ends)
/// let interval_2 = grid.interval(&IntervalId::new(2));
/// assert!(interval_2.contains_point(&2.0)); // ✓ includes left bound
/// assert!(interval_2.contains_point(&3.0));  // ✓ includes right bound
/// ```
///
/// ### Semi-Open Domain `[a, b)`
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalLowerClosedUpperOpen::new(0.0, 2.0),
///     NumIntervals::try_new(2).unwrap()
/// );
///
/// // First interval: [0, 1) (closed on lower bound, open on upper bound for first interval)
/// let interval_0 = grid.interval(&IntervalId::new(0));
/// assert!(interval_0.contains_point(&0.0)); // ✓ includes lower bound
/// assert!(!interval_0.contains_point(&1.0)); // ✗ excludes upper bound
///
/// // Last interval: [1, 2) (respects open upper bound of domain)
/// let interval_1 = grid.interval(&IntervalId::new(1));
/// assert!(interval_1.contains_point(&1.0)); // ✓ includes lower bound
/// assert!(!interval_1.contains_point(&2.0)); // ✗ excludes upper bound
/// ```
///
/// ## Specialized Uniform Grid Operations
///
/// ### Optimized Range Queries
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalClosed::new(0.0, 10.0),
///     NumIntervals::try_new(1000).unwrap()
/// );
///
/// // For uniform grids, range queries can be analytical
/// fn find_intervals_in_range(
///     grid: &Grid1DUniform<IntervalClosed<f64>>,
///     x_min: f64,
///     x_max: f64
/// ) -> Vec<IntervalId> {
///     let delta = *grid.delta_points().as_ref();
///     let domain_start = *grid.domain().lower_bound_value();
///     
///     // Analytical calculation - O(1) instead of O(log n)
///     let start_interval = ((x_min - domain_start) / delta).floor() as usize;
///     let end_interval = ((x_max - domain_start) / delta).ceil() as usize;
///     
///     let start_clamped = start_interval.min(*grid.num_intervals().as_ref() - 1);
///     let end_clamped = end_interval.min(*grid.num_intervals().as_ref() - 1);
///     
///     (start_clamped..=end_clamped)
///         .map(IntervalId::new)
///         .collect()
/// }
///
/// let range_intervals = find_intervals_in_range(&grid, 2.5, 7.3);
/// println!("Range [2.5, 7.3] spans {} intervals", range_intervals.len());
/// ```
///
/// ### Analytical Derivatives and Integration
/// ```rust
/// use grid1d::{
///     Grid1DUniform, HasCoords1D,
///     intervals::*,
///     scalars::{NumIntervals, CoordId},
/// };
/// use try_create::TryNew;
///
/// let grid = Grid1DUniform::new(
///     IntervalClosed::new(0.0, 1.0),
///     NumIntervals::try_new(100).unwrap()
/// );
///
/// // Function values on the grid
/// let coords = grid.coords();
/// let values: Vec<f64> = coords.iter().map(|&x| x.powi(2)).collect();
///
/// // Analytical integration using trapezoidal rule
/// let dx = *grid.delta_points().as_ref();
/// let integral = dx * (values[0]/2.0 + values.iter().skip(1).take(values.len()-2).sum::<f64>() + values.last().unwrap()/2.0);
/// let analytical = 1.0/3.0; // ∫₀¹ x² dx = 1/3
///
/// println!("Numerical integral: {:.6}", integral);
/// println!("Analytical result: {:.6}", analytical);
/// println!("Error: {:.2e}", (integral - analytical).abs());
///
/// // Finite difference derivatives
/// let mut derivatives = vec![0.0; values.len()];
/// for i in 1..values.len()-1 {
///     derivatives[i] = (values[i+1] - values[i-1]) / (2.0 * dx);
/// }
///
/// // At x=0.5, analytical derivative is 2*0.5 = 1.0
/// let mid_idx = values.len() / 2;
/// let numerical_deriv = derivatives[mid_idx];
/// let analytical_deriv = 2.0 * coords[CoordId::new(mid_idx)];
/// println!("Derivative at x=0.5: numerical={:.6}, analytical={:.6}",
///          numerical_deriv, analytical_deriv);
/// ```
///
/// ## Error Handling and Validation
///
/// ### Construction Validation
/// ```rust
/// use grid1d::{
///     Grid1DTrait, HasIntervalIdRange, Grid1DUniform, HasCoords1D,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use try_create::TryNew;
///
/// // Valid construction
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let num_intervals = NumIntervals::try_new(10).unwrap();
/// let grid = Grid1DUniform::new(domain, num_intervals);
///
/// // The constructor doesn't fail for Grid1DUniform since it computes
/// // coordinates analytically, but domain validation happens in the
/// // IntervalClosed construction
/// assert_eq!(grid.coords().len(), 11); // n+1 points
/// assert_eq!(grid.num_intervals().as_ref(), &10);
/// ```
///
/// ### Robust Integration Pattern
/// ```rust
/// use grid1d::{
///     Grid1D, Grid1DUniform, HasCoords1D, Grid1DTrait, HasIntervalIdRange,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use try_create::TryNew;
///
/// fn create_simulation_grid(
///     x_min: f64,
///     x_max: f64,
///     target_resolution: f64
/// ) -> Result<Grid1DUniform<IntervalClosed<f64>>, String> {
///     if x_min >= x_max {
///         return Err("Invalid domain: x_min must be less than x_max".to_string());
///     }
///     
///     if target_resolution <= 0.0 {
///         return Err("Invalid resolution: must be positive".to_string());
///     }
///     
///     let domain_length = x_max - x_min;
///     let num_intervals = (domain_length / target_resolution).ceil() as usize;
///     
///     if num_intervals < 1 {
///         return Err("Domain too small for target resolution".to_string());
///     }
///     
///     let domain = IntervalClosed::new(x_min, x_max);
///     let num_intervals = NumIntervals::try_new(num_intervals)
///         .map_err(|e| format!("Invalid number of intervals: {}", e))?;
///     
///     Ok(Grid1DUniform::new(domain, num_intervals))
/// }
///
/// // Usage
/// let grid = create_simulation_grid(0.0, 1.0, 0.01).unwrap();
/// println!("Created grid with {} intervals, actual spacing: {:.6}",
///          grid.num_intervals().as_ref(), grid.delta_points().as_ref());
///
/// // Error cases
/// assert!(create_simulation_grid(1.0, 0.0, 0.01).is_err()); // Invalid domain
/// assert!(create_simulation_grid(0.0, 1.0, -0.01).is_err()); // Invalid resolution
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`Grid1DUniform`]
/// - **Finite difference methods**: Equal spacing required for standard stencils
/// - **Spectral methods**: Uniform grids for Fourier transforms and Chebyshev methods
/// - **Regular sampling**: Data collection or function evaluation at uniform intervals
/// - **Performance-critical applications**: When O(1) point location is essential
/// - **Memory-constrained systems**: When uniform spacing allows analytical computation
/// - **Simple problems**: When uniform resolution is sufficient across the domain
///
/// ### When to Use [`Grid1DNonUniform`] Instead
/// - **Adaptive mesh refinement**: Variable resolution based on solution gradients
/// - **Boundary layers**: High resolution near walls or interfaces
/// - **Multi-scale problems**: Different length scales in the same domain
/// - **Complex geometries**: Non-uniform point distributions for irregular shapes
/// - **Error-based refinement**: Points placed based on error indicators
/// - **Shock capturing**: Concentrated points near discontinuities
///
/// ## Mathematical Properties and Guarantees
///
/// The [`Grid1DUniform`] struct maintains these mathematical properties:
///
/// ### Uniform Spacing Property
/// ```text
/// ∀i ∈ [0, n-1]: length(Iᵢ) = δ = domain_length / n
/// ```
/// All sub-intervals have exactly the same length.
///
/// ### Perfect Regularity
/// ```text
/// ∀i ∈ [0, n]: pᵢ = lower_bound + i * δ
/// ```
/// Coordinates follow a perfect arithmetic progression.
///
/// ### Optimal Uniformity
/// ```text
/// uniformity_ratio = max_length / min_length = 1.0
/// ```
/// The uniformity ratio is always exactly 1.0.
///
/// ### Boundary Consistency
/// ```text
/// p₀ = domain.lower_bound ∧ pₙ = domain.upper_bound
/// ```
/// First and last points exactly match domain boundaries.
///
/// ### Analytical Predictability
/// ```text
/// point_location(x) = floor((x - lower_bound) / δ)  [O(1) operation]
/// ```
/// Point location can be computed analytically without search.
///
/// ### Cache Efficiency
/// ```text
/// memory_pattern = contiguous_array[n+1]  [optimal for CPU cache]
/// ```
/// Coordinates are stored in a contiguous array for maximum cache efficiency.
///
/// These properties make [`Grid1DUniform`] ideal for numerical methods that benefit
/// from regular spacing and predictable access patterns.
///
/// ## See Also
///
/// - [`Grid1DNonUniform`]: Variable spacing alternative for adaptive applications
/// - [`Grid1D`]: Unified interface supporting both uniform and non-uniform partitions
/// - [`Grid1DTrait`]: Core trait providing partition operations
/// - [`Coords1D`]: Underlying coordinate container
/// - [`HasDomain1D`], [`HasCoords1D`]: Core domain access traits
/// - [`NumIntervals`]: Type for specifying the number of intervals
/// - [`PositiveRealScalar`]: Type for the uniform spacing value (`delta_points`)
#[derive(Debug, Clone, Getters, PartialEq, Serialize, Deserialize)]
#[serde(bound(deserialize = "Domain1D: for<'d> serde::Deserialize<'d>, \
                  Domain1D::RealType: for<'d> serde::Deserialize<'d>"))]
pub struct Grid1DUniform<Domain1D: IntervalFinitePositiveLengthTrait> {
    /// The points that defines the partition of the domain.
    coords: Coords1D<Domain1D::RealType>,

    /// Domain covered by the grid.
    domain: Domain1D,

    /// The index spaces associated with the grid, which define the topology and interval indexing.
    index_spaces: Grid1DIndexSpaces,

    /// The distance between two consecutive points (i.e. the length of the interval defined by two consecutive points).
    #[getset(get = "pub")]
    delta_points: PositiveRealScalar<Domain1D::RealType>,
}

fn compute_delta_points<Domain1D: IntervalFinitePositiveLengthTrait>(
    domain: &Domain1D,
    num_intervals: &NumIntervals,
) -> PositiveRealScalar<Domain1D::RealType> {
    let n_intervals = *num_intervals.as_ref();

    PositiveRealScalar::try_new(
        domain.length().into_inner()
            * Domain1D::RealType::try_from_f64(n_intervals as f64)
                .unwrap()
                .try_reciprocal()
                .unwrap(),
    )
    .expect("Delta points must be positive!")
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> Grid1DUniform<Domain1D> {
    /// Creates a non-periodic uniform grid on a real-line topology.
    ///
    /// # Parameters
    ///
    /// - `domain`: Bounded interval domain to partition.
    /// - `num_intervals`: Number of equal sub-intervals.
    ///
    /// # Returns
    ///
    /// A [`Grid1DUniform`] with [`Topology1D::RealLine`].
    #[must_use]
    pub fn new(domain: Domain1D, num_intervals: NumIntervals) -> Self {
        let delta_points = compute_delta_points(&domain, &num_intervals);

        let coords = Coords1D::new_uniform(&domain, &num_intervals);

        let index_spaces = Grid1DIndexSpaces::new_non_periodic(num_intervals);

        Self {
            coords,
            domain,
            delta_points,
            index_spaces,
        }
    }
}

impl<Domain1D> Grid1DUniform<Domain1D>
where
    Domain1D: IntervalFinitePositiveLengthTrait + SupportsCircularTopology,
{
    /// Creates a periodic uniform grid with circular topology.
    ///
    /// This constructor is available only for domain interval types that support
    /// circular topology (currently half-open bounded intervals).
    ///
    /// # Parameters
    ///
    /// - `domain`: Bounded interval domain to partition.
    /// - `num_intervals`: Number of equal sub-intervals.
    ///
    /// # Returns
    ///
    /// A [`Grid1DUniform`] with [`Topology1D::Circle`].
    #[must_use]
    pub fn new_periodic(domain: Domain1D, num_intervals: NumIntervals) -> Self {
        let delta_points = compute_delta_points(&domain, &num_intervals);

        let coords = Coords1D::new_uniform(&domain, &num_intervals);

        let index_spaces = if domain.is_lower_bound_open() {
            Grid1DIndexSpaces::new_periodic_lower_open_interval(num_intervals)
        } else {
            Grid1DIndexSpaces::new_periodic_upper_open_interval(num_intervals)
        };

        Self {
            coords,
            domain,
            delta_points,
            index_spaces,
        }
    }
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> HasDomain1D for Grid1DUniform<Domain1D> {
    type Domain1D = Domain1D;

    #[inline(always)]
    fn domain(&self) -> &Domain1D {
        &self.domain
    }
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> HasCoords1D for Grid1DUniform<Domain1D> {
    type CoordType = Domain1D::RealType;

    /// Get the coordinates of the 1D points.
    ///
    /// # Note
    /// The returned values does not contain duplicates and are sorted in ascending order.
    #[inline(always)]
    fn coords(&self) -> &Coords1D<Self::CoordType> {
        &self.coords
    }
}

impl<Domain1D: Grid1DIntervalBuilder> FindIntervalIdOfPoint for Grid1DUniform<Domain1D> {
    type Point1DType = Domain1D::RealType;

    /// Optimized O(1) point location for uniform grids.
    ///
    /// This implementation leverages the uniform spacing to compute the interval ID analytically,
    /// while ensuring consistency with the boundary point assignment rules defined by the
    /// [`Grid1DTrait`] trait.
    ///
    /// # Boundary Point Assignment Rules
    ///
    /// - **Interior partition points**: For left-closed sub-intervals (`[a,b]` and `[a,b)` domains),
    ///   `pₖ` belongs to the **right** interval; for right-closed sub-intervals, to the **left** interval.
    /// - **Domain lower bound**: Belongs to the first interval (ID 0)
    /// - **Domain upper bound**: Belongs to the last interval (ID n-1)
    ///
    /// # Example
    /// ```rust
    /// use grid1d::{*, intervals::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1DUniform::new(
    ///     IntervalClosed::new(0.0, 2.0),
    ///     NumIntervals::try_new(2).unwrap()
    /// );
    /// // Grid points: [0.0, 1.0, 2.0]
    /// // Intervals: [0.0, 1.0), [1.0, 2.0]  (left-closed for [a,b] domain)
    ///
    /// assert_eq!(grid.find_interval_id_of_point(&0.0).unwrap(), IntervalId::new(0)); // Domain lower bound
    /// assert_eq!(grid.find_interval_id_of_point(&0.5).unwrap(), IntervalId::new(0)); // Interior point
    /// assert_eq!(grid.find_interval_id_of_point(&1.0).unwrap(), IntervalId::new(1)); // p₁ → right interval (left-closed)
    /// assert_eq!(grid.find_interval_id_of_point(&1.5).unwrap(), IntervalId::new(1)); // Interior point
    /// assert_eq!(grid.find_interval_id_of_point(&2.0).unwrap(), IntervalId::new(1)); // Domain upper bound
    /// ```
    fn find_interval_id_of_point(&self, x: &Self::Point1DType) -> Option<IntervalId> {
        let domain = self.domain();

        if !domain.contains_point(x) {
            return None;
        }

        let num_intervals = *self.num_intervals().as_ref();
        // Check if x is exactly equal to any grid coordinate
        match self.coords().find_index(x) {
            Some(coord_id) => {
                let coord_id = coord_id.into_inner();

                // Exact match with a grid point
                if coord_id == 0 {
                    // This corresponds to the domain lower bound
                    Some(IntervalId::new(0))
                } else if coord_id == num_intervals {
                    // This corresponds to the domain upper bound
                    Some(IntervalId::new(num_intervals - 1))
                } else {
                    // Here we use the domain's lower bound closure property to determine the assignment of interior partition points.
                    // This is because the domain type dictates the sub-interval construction rules, which in turn determine how boundary points are assigned to intervals.
                    // For example:
                    // - For a closed domain [a, b], all sub-intervals are left-closed, so interior points belong to the right interval.
                    // - For a half-open domain (a, b], all sub-intervals are right-closed, so interior points belong to the left interval.
                    // - For a half-open domain [a, b), the first interval is left-closed and the last interval is open, so interior points belong to the right interval.
                    // - For a open domain (a, b), all sub-intervals are open, so interior points do not belong to any interval, but since we are checking for exact matches with grid points, we can still assign them based on the closure of the lower bound.
                    // This logic ensures that the point location is consistent with the interval construction rules defined by the domain type, while still leveraging the uniform spacing for O(1) access.
                    // The key point is that the domain's lower bound closure property directly informs how we assign interior partition points to intervals, which is crucial for maintaining the correct semantics of the partitioning.
                    if self.domain().is_lower_bound_closed() {
                        Some(IntervalId::new(coord_id))
                    } else {
                        Some(IntervalId::new(coord_id - 1))
                    }
                }
            }
            None => {
                // No exact match, so x lies strictly within an interval

                // Compute the analytical interval index
                let domain_start = domain.lower_bound_value();
                let delta = self.delta_points.as_ref();
                let offset = x.clone() - domain_start;
                let offset_div_delta = offset / delta;

                // Convert to usize and clamp to valid range
                let index = offset_div_delta
                    .kernel_floor()
                    .truncate_to_usize()
                    .unwrap()
                    .min(num_intervals - 1);

                Some(IntervalId::new(index))
            }
        }
    }
}

impl<Domain1D: Grid1DIntervalBuilder> HasIntervalIdRange for Grid1DUniform<Domain1D> {
    #[inline(always)]
    fn first_interval_id(&self) -> IntervalId {
        self.index_spaces.interval_index_space().first_interval_id()
    }

    #[inline(always)]
    fn last_interval_id(&self) -> IntervalId {
        self.index_spaces.interval_index_space().last_interval_id()
    }
}

impl<Domain1D: Grid1DIntervalBuilder> HasCoordIdRange for Grid1DUniform<Domain1D> {
    #[inline(always)]
    fn first_coord_id(&self) -> CoordId {
        self.index_spaces.coord_index_space().first_coord_id()
    }

    #[inline(always)]
    fn last_coord_id(&self) -> CoordId {
        self.index_spaces.coord_index_space().last_coord_id()
    }
}

impl<Domain1D: Grid1DIntervalBuilder> Grid1DTrait for Grid1DUniform<Domain1D> {
    #[inline(always)]
    fn index_spaces(&self) -> &Grid1DIndexSpaces {
        &self.index_spaces
    }

    /// Returns the maximum interval length in the partition.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 1.5, 3.0]);
    /// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
    /// let max_length = grid.max_interval_length();
    /// assert_eq!(max_length.as_ref(), &1.5); // Interval [1.5, 3.0] has length 1.5
    /// ```
    #[inline(always)]
    fn max_interval_length(&self) -> PositiveRealScalar<Self::Point1DType> {
        self.delta_points.clone()
    }

    /// Returns the minimum interval length in the partition.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 1.5, 3.0]);
    /// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
    /// let min_length = grid.min_interval_length();
    /// assert_eq!(min_length.as_ref(), &0.5); // Interval [1.0, 1.5] has length 0.5
    /// ```
    #[inline(always)]
    fn min_interval_length(&self) -> PositiveRealScalar<Self::Point1DType> {
        self.delta_points.clone()
    }

    type UniformlyRefinedGrid1DType = Self;

    fn refine_uniform(
        self,
        num_extra_points_each_interval: &PositiveNumPoints1D,
    ) -> Grid1DUniformRefinement<Self> {
        let num_old_intervals = *self.num_intervals().as_ref();
        let num_sub_intervals = *num_extra_points_each_interval.as_ref() + 1;
        let num_new_intervals = num_old_intervals * num_sub_intervals;

        let refined_grid = Grid1DUniform::new(
            self.domain.clone(),
            NumIntervals::try_new(num_new_intervals).unwrap(),
        );

        // Maps each interval in the refined grid to the corresponding interval in the original grid.
        let mut refined_to_original_interval_mapping: Vec<IntervalId> =
            Vec::with_capacity(num_new_intervals);

        // Maps each interval in the original grid to the corresponding intervals in the refined grid.
        let mut original_to_refined_interval_mapping: Vec<Vec<IntervalId>> =
            Vec::with_capacity(num_old_intervals);

        let mut refined_id = 0;
        for original_id in 0..num_old_intervals {
            refined_to_original_interval_mapping.resize(
                refined_to_original_interval_mapping.len() + num_sub_intervals,
                IntervalId::new(original_id),
            );

            original_to_refined_interval_mapping.push(
                (refined_id..refined_id + num_sub_intervals)
                    .map(IntervalId::new)
                    .collect(),
            );

            refined_id += num_sub_intervals;
        }

        Grid1DUniformRefinement::new(
            refined_grid,
            self,
            refined_to_original_interval_mapping,
            original_to_refined_interval_mapping,
        )
    }

    fn refine(
        self,
        intervals_to_refine: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
    ) -> Grid1DNonUniformRefinement<Self> {
        build_non_uniform_grid_refinement(self, intervals_to_refine)
    }
}
//-------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------
/// An enum representing a one-dimensional grid, which defines a partition of a bounded interval with positive length.
///
/// The [`Grid1D`] enum models a set of points that partition a bounded interval into adjacent, non-overlapping subintervals.
/// It supports both ***uniform grids***, where all subintervals have the same length, and ***non-uniform grids***, where
/// subintervals can have varying lengths.
///
/// # Type Parameters
///
/// - `Point1DType`: The type of the points in the grid. It must implement the [`RealScalar`] trait, ensuring that the
///   points are real numbers with additional mathematical properties.
///
/// # Variants
///
/// - `Uniform`: Represents a grid where all subintervals have the same length. This is modeled by the
///   [`Grid1DUniform`] struct.
/// - `NonUniform`: Represents a grid where subintervals can have varying lengths. This is modeled by the
///   [`Grid1DNonUniform`] struct.
///
/// # Features
///
/// - **Uniform Grids**: All subintervals have the same length, making them suitable for evenly spaced computations.
/// - **Non-Uniform Grids**: Subintervals can have varying lengths, allowing for more flexibility in grid design.
/// - **Integration with Grid1DTrait**: The [`Grid1D`] enum implements the [`Grid1DTrait`] trait, providing
///   methods for accessing intervals and performing operations on the grid.
///
/// # Examples
///
/// ## Creating a Uniform Grid
///
/// ```rust
/// use grid1d::{
///     Grid1D, HasCoords1D, Grid1DTrait, HasIntervalIdRange,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use std::ops::Deref;
/// use try_create::TryNew;
///
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let grid = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());
///
/// assert_eq!(grid.coords().deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
/// assert_eq!(grid.num_intervals().as_ref(), &4);
/// ```
///
/// ## Creating a Non-Uniform Grid
///
/// ```rust
/// use grid1d::{
///     Coords1D, Grid1D, HasCoords1D,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// let coords = SortedSet::from_unsorted(vec![0.0, 0.5, 1.0, 2.0]);
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
/// assert_eq!(grid.coords().deref(), &[0.0, 0.5, 1.0, 2.0]);
/// ```
///
/// ## Handling Errors
///
/// ```rust
/// use grid1d::{
///     Grid1D, ErrorsGrid1D,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// let coords = SortedSet::from_unsorted(vec![0.0]); // Only one point provided
/// let err = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap_err();
/// assert!(matches!(err, ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { .. }), "{:?}", err);
/// ```
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
#[serde(bound(deserialize = "Domain1D: for<'d> serde::Deserialize<'d>, \
                  Domain1D::RealType: for<'d> serde::Deserialize<'d>"))]
pub enum Grid1D<Domain1D: IntervalFinitePositiveLengthTrait> {
    /// All the intervals defined by the [`Grid1D`] object have the same length.
    Uniform(Grid1DUniform<Domain1D>),

    /// The intervals defined by the [`Grid1D`] object have not the same length.
    NonUniform(Grid1DNonUniform<Domain1D>),
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> Grid1D<Domain1D> {
    /// Creates a new uniform grid over the specified domain with the given number of intervals.
    ///
    /// # Parameters
    ///
    /// - `domain`: The bounded interval over which the grid is defined.
    /// - `num_intervals`: The number of subintervals in the grid.
    ///
    /// # Returns
    ///
    /// A [`Grid1D`] instance representing a uniform grid.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Grid1D, HasCoords1D, Grid1DTrait,
    ///     intervals::*,
    ///     scalars::NumIntervals,
    /// };
    /// use std::ops::Deref;
    /// use try_create::TryNew;
    ///
    /// let domain = IntervalClosed::new(0., 1.);
    /// let grid = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());
    ///
    /// assert_eq!(grid.coords().deref(), &[0., 0.25, 0.5, 0.75, 1.]);
    /// ```
    #[must_use]
    pub fn uniform(domain: Domain1D, num_intervals: NumIntervals) -> Self {
        Grid1D::Uniform(Grid1DUniform::new(domain, num_intervals))
    }

    /// Creates a new grid from a sorted set of points (at least 2).
    ///
    /// # Parameters
    ///
    /// - `values`: A sorted set of points defining the grid.
    ///
    /// # Returns
    ///
    /// - `Ok(Self)`: If at least two points are provided.
    /// - [`Err<ErrorsGrid1D>`]: If fewer than two points are provided.
    ///
    /// # Errors
    ///
    /// - Returns [`ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints`] if fewer than two points are provided.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Grid1D, HasCoords1D,
    ///     intervals::*,
    /// };
    /// use sorted_vec::partial::SortedSet;
    /// use std::ops::Deref;
    ///
    /// let coords = SortedSet::from_unsorted(vec![0., 0.5, 1., 2.]);
    /// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
    /// assert_eq!(grid.coords().deref(), &[0., 0.5, 1., 2.]);
    /// ```
    pub fn try_from_sorted(
        values: SortedSet<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>>
    where
        Domain1D: IntervalFromBounds,
    {
        Self::try_from_coords(Coords1D::try_from(values)?)
    }

    /// Creates a new grid from a domain and pre-validated coordinates.
    ///
    /// This is a more efficient constructor when you already have a [`Coords1D`]
    /// object, avoiding redundant validation.
    pub fn try_from_coords(
        coords: Coords1D<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>>
    where
        Domain1D: IntervalFromBounds,
    {
        let num_coords = coords.num_points();
        if *num_coords.as_ref() < 2 {
            return Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints {
                num_points_provided: num_coords,
                backtrace: capture_backtrace(),
            });
        }

        let domain = Domain1D::try_new(coords.first().clone(), coords.last().clone())?;

        Self::try_new(domain, coords)
    }

    /// Creates a new [`Grid1D`] from an explicit domain and pre-validated coordinates.
    ///
    /// Validates that the coordinates are consistent with the domain (first and last
    /// coordinate must match the domain bounds, and at least two distinct points must
    /// be present). Prefer [`Grid1D::try_from_sorted`] or [`Grid1D::try_from_coords`]
    /// when the domain can be inferred automatically.
    pub fn try_new(
        domain: Domain1D,
        coords: Coords1D<Domain1D::RealType>,
    ) -> Result<Self, ErrorsGrid1D<Domain1D>> {
        let num_intervals = validate_coords_in_domain_closure(&coords, &domain)?;

        let n_pts = coords.num_points();

        // In order to build a NumIntervals object, we need at least 2 distinct points
        if n_pts.as_ref() < &2 {
            return Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints {
                num_points_provided: n_pts,
                backtrace: capture_backtrace(),
            });
        }
        let n_pts = n_pts.into_inner();

        // Check if the coordinates are uniform
        let coords_slice = coords.as_ref();
        let length_first_interval = coords_slice[1].clone() - &coords_slice[0];

        // Use a small absolute tolerance for cases where interval lengths are near zero.
        let rel_tol = RelativeTolerance::epsilon();

        let is_uniform = (2..n_pts).all(|i| {
            let length_current_interval = coords_slice[i].clone() - &coords_slice[i - 1];
            approx_eq(&length_current_interval, &length_first_interval, &rel_tol)
        });

        if is_uniform {
            Ok(Grid1D::Uniform(Grid1DUniform::new(domain, num_intervals)))
        } else {
            let index_spaces = Grid1DIndexSpaces::new_non_periodic(num_intervals);

            let grid1d = Grid1DNonUniform {
                domain,
                coords,
                index_spaces,
            };

            Ok(Grid1D::NonUniform(grid1d))
        }
    }
}

/// **INTERNAL** Checks approximate equality using relative tolerance.
///
/// This function determines if two real numbers are approximately equal by comparing
/// their relative difference to a specified tolerance. It's used internally to
/// detect uniform spacing in coordinate arrays.
///
/// ## Mathematical Definition
/// Two numbers `a` and `b` are considered approximately equal if:
/// ```text
/// |a - b| ≤ max(|a|, |b|) * relative_tolerance
/// ```
///
/// ## Parameters
/// - `a`: First value to compare
/// - `b`: Second value to compare  
/// - `rel_tol`: Relative tolerance for the comparison
///
/// ## Returns
/// `true` if the values are approximately equal within the specified tolerance.
///
/// ## Use Cases
/// - **Uniform grid detection**: Checking if intervals have consistent spacing
/// - **Numerical comparison**: Handling floating-point precision issues
/// - **Grid classification**: Determining if a grid should be treated as uniform
fn approx_eq<T: RealScalar>(a: &T, b: &T, rel_tol: &RelativeTolerance<T>) -> bool {
    let diff = (a.clone() - b).abs();

    // Use the max of the absolute values for the relative comparison
    let abs_a = a.clone().abs();
    let abs_b = b.clone().abs();
    let max_abs = abs_a.max_by_ref(&abs_b).clone();
    diff <= max_abs * rel_tol.as_ref()
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> HasDomain1D for Grid1D<Domain1D> {
    type Domain1D = Domain1D;

    #[inline]
    fn domain(&self) -> &Domain1D {
        match self {
            Grid1D::Uniform(pts) => pts.domain(),
            Grid1D::NonUniform(pts) => pts.domain(),
        }
    }
}

impl<Domain1D: IntervalFinitePositiveLengthTrait> HasCoords1D for Grid1D<Domain1D> {
    type CoordType = Domain1D::RealType;

    #[inline]
    fn coords(&self) -> &Coords1D<Self::CoordType> {
        match self {
            Grid1D::Uniform(pts) => pts.coords(),
            Grid1D::NonUniform(pts) => pts.coords(),
        }
    }
}

impl<Domain1D: Grid1DIntervalBuilder> FindIntervalIdOfPoint for Grid1D<Domain1D> {
    type Point1DType = Domain1D::RealType;

    /// Dispatches to the inner grid's point-location method.
    ///
    /// For the [`Grid1D::Uniform`] variant this resolves to the O(1) analytical
    /// override defined on [`Grid1DUniform`]; for [`Grid1D::NonUniform`] it uses
    /// the default O(log n) binary-search implementation.
    fn find_interval_id_of_point(&self, x: &Self::Point1DType) -> Option<IntervalId> {
        match self {
            Grid1D::Uniform(grid) => grid.find_interval_id_of_point(x),
            Grid1D::NonUniform(grid) => grid.find_interval_id_of_point(x),
        }
    }
}

impl<Domain1D: Grid1DIntervalBuilder> HasIntervalIdRange for Grid1D<Domain1D> {
    #[inline]
    fn first_interval_id(&self) -> IntervalId {
        match self {
            Grid1D::Uniform(grid) => grid.first_interval_id(),
            Grid1D::NonUniform(grid) => grid.first_interval_id(),
        }
    }

    #[inline]
    fn last_interval_id(&self) -> IntervalId {
        match self {
            Grid1D::Uniform(grid) => grid.last_interval_id(),
            Grid1D::NonUniform(grid) => grid.last_interval_id(),
        }
    }
}

impl<Domain1D: Grid1DIntervalBuilder> HasCoordIdRange for Grid1D<Domain1D> {
    #[inline]
    fn first_coord_id(&self) -> CoordId {
        match self {
            Grid1D::Uniform(grid) => grid.first_coord_id(),
            Grid1D::NonUniform(grid) => grid.first_coord_id(),
        }
    }

    #[inline]
    fn last_coord_id(&self) -> CoordId {
        match self {
            Grid1D::Uniform(grid) => grid.last_coord_id(),
            Grid1D::NonUniform(grid) => grid.last_coord_id(),
        }
    }
}

impl<Domain1D: Grid1DIntervalBuilder> Grid1DTrait for Grid1D<Domain1D> {
    #[inline]
    fn index_spaces(&self) -> &Grid1DIndexSpaces {
        match self {
            Grid1D::Uniform(grid) => grid.index_spaces(),
            Grid1D::NonUniform(grid) => grid.index_spaces(),
        }
    }

    type UniformlyRefinedGrid1DType = Self;

    fn refine_uniform(
        self,
        num_extra_points_each_interval: &PositiveNumPoints1D,
    ) -> Grid1DUniformRefinement<Self> {
        match self {
            Grid1D::Uniform(grid_uniform) => grid_uniform
                .refine_uniform(num_extra_points_each_interval)
                .into(),
            Grid1D::NonUniform(grid_non_uniform) => grid_non_uniform
                .refine_uniform(num_extra_points_each_interval)
                .into(),
        }
    }

    fn refine(
        self,
        intervals_to_refine: &std::collections::BTreeMap<IntervalId, PositiveNumPoints1D>,
    ) -> Grid1DNonUniformRefinement<Self> {
        match self {
            Grid1D::Uniform(grid_uniform) => grid_uniform.refine(intervals_to_refine).into(),
            Grid1D::NonUniform(grid_non_uniform) => {
                grid_non_uniform.refine(intervals_to_refine).into()
            }
        }
    }
}

// ============================================================================
// From trait implementations for Grid1D enum
// ============================================================================

/// Converts a [`Grid1DUniform`] into a [`Grid1D`] enum variant.
///
/// This allows seamless use of uniform grids in contexts expecting
/// the unified [`Grid1D`] interface.
///
/// # Example
///
/// ```rust
/// use grid1d::{Grid1D, Grid1DUniform, intervals::*, scalars::NumIntervals};
/// use try_create::TryNew;
///
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let uniform_grid = Grid1DUniform::new(domain, NumIntervals::try_new(10).unwrap());
///
/// // Using From::from
/// let grid: Grid1D<IntervalClosed<f64>> = Grid1D::from(uniform_grid.clone());
///
/// // Using .into()
/// let grid: Grid1D<IntervalClosed<f64>> = uniform_grid.into();
/// ```
impl<Domain1D: IntervalFinitePositiveLengthTrait> From<Grid1DUniform<Domain1D>>
    for Grid1D<Domain1D>
{
    fn from(grid: Grid1DUniform<Domain1D>) -> Self {
        Grid1D::Uniform(grid)
    }
}

/// Converts a [`Grid1DNonUniform`] into a [`Grid1D`] enum variant.
///
/// This allows seamless use of non-uniform grids in contexts expecting
/// the unified [`Grid1D`] interface.
///
/// # Example
///
/// ```rust
/// use grid1d::{Grid1D, Grid1DNonUniform, coords::Coords1D, intervals::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.3, 0.7, 1.0])).unwrap();
/// let non_uniform_grid = Grid1DNonUniform::<IntervalClosed<f64>>::try_new_from_coords(coords).unwrap();
///
/// // Using From::from
/// let grid: Grid1D<IntervalClosed<f64>> = Grid1D::from(non_uniform_grid.clone());
///
/// // Using .into()
/// let grid: Grid1D<IntervalClosed<f64>> = non_uniform_grid.into();
/// ```
impl<Domain1D: IntervalFinitePositiveLengthTrait> From<Grid1DNonUniform<Domain1D>>
    for Grid1D<Domain1D>
{
    fn from(grid: Grid1DNonUniform<Domain1D>) -> Self {
        Grid1D::NonUniform(grid)
    }
}
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
#[cfg(test)]
mod tests {
    use crate::{coords::*, grids::*, intervals::*, *};
    use std::ops::Deref;

    #[cfg(feature = "rug")]
    use num_valid::vec_f64_into_vec_real;

    mod native64_strict_finite {
        use super::*;
        use num_valid::RealNative64StrictFinite;

        type Real = RealNative64StrictFinite;

        mod grid1d {

            use super::*;

            #[test]
            fn test_grid1d_uniform_creation() {
                let domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());

                assert_eq!(grid.domain(), &domain);
                assert_eq!(
                    grid.coords().deref(),
                    &[
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(0.25).unwrap(),
                        Real::try_new(0.5).unwrap(),
                        Real::try_new(0.75).unwrap(),
                        Real::try_new(1.0).unwrap()
                    ]
                );
            }

            #[test]
            fn test_grid1d_non_uniform_creation() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(0.5).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords.clone()).unwrap();

                assert_eq!(grid.coords().deref(), &coords.to_vec());
            }

            #[test]
            fn test_grid1d_insufficient_points() {
                let coords = SortedSet::from_unsorted(vec![Real::try_new(0.0).unwrap()]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords);
                matches!(
                    grid,
                    Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { .. })
                );
            }

            #[test]
            fn grid1d_01() {
                let v = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]);
                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();
                assert_eq!(grid1d.coords().deref(), &v.to_vec());
            }

            #[test]
            fn interval_01() {
                let v = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]);

                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();

                let interval = grid1d.interval(&IntervalId::new(0));

                match interval {
                    SubIntervalInPartition::Single(interval) => {
                        assert_eq!(
                            TryInto::<IntervalClosed<Real>>::try_into(interval).unwrap(),
                            IntervalClosed::new(
                                Real::try_new(0.0).unwrap(),
                                Real::try_new(1.0).unwrap()
                            )
                        );
                    }
                    _ => {
                        panic!("There should be only one interval")
                    }
                }
            }

            #[test]
            #[cfg(debug_assertions)]
            #[should_panic]
            fn interval_02() {
                let v = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]);
                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();
                let _interval = grid1d.interval(&IntervalId::new(1));
            }

            #[test]
            fn find_interval_id_01() {
                let v = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();

                let i0 = grid1d
                    .find_interval_id_of_point(&Real::try_new(0.0).unwrap())
                    .unwrap();
                assert_eq!(i0, IntervalId::new(0));

                let i1 = grid1d
                    .find_interval_id_of_point(&Real::try_new(0.5).unwrap())
                    .unwrap();
                assert_eq!(i1, IntervalId::new(0));

                let i2 = grid1d
                    .find_interval_id_of_point(&Real::try_new(1.0).unwrap())
                    .unwrap();
                assert_eq!(i2, IntervalId::new(1)); // p₁ → right interval

                let i3 = grid1d
                    .find_interval_id_of_point(&Real::try_new(1.5).unwrap())
                    .unwrap();
                assert_eq!(i3, IntervalId::new(1));

                let i4 = grid1d
                    .find_interval_id_of_point(&Real::try_new(2.0).unwrap())
                    .unwrap();
                assert_eq!(i4, IntervalId::new(2)); // p₂ → right interval

                let i5 = grid1d
                    .find_interval_id_of_point(&Real::try_new(2.5).unwrap())
                    .unwrap();
                assert_eq!(i5, IntervalId::new(2));

                let i6 = grid1d
                    .find_interval_id_of_point(&Real::try_new(3.0).unwrap())
                    .unwrap();
                assert_eq!(i6, IntervalId::new(2));
            }

            #[test]
            fn uniform_01() {
                let interval =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let grid1d = Grid1D::uniform(interval, NumIntervals::try_new(4).unwrap());

                let expected_coords = vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(0.25).unwrap(),
                    Real::try_new(0.5).unwrap(),
                    Real::try_new(0.75).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ];

                assert_eq!(grid1d.coords().deref(), &expected_coords);
                //        assert_eq!(grid1d.coords(), &Grid1DCoords::Uniform);
            }

            #[test]
            fn grid_intersection_01() {
                let grid1d_a = Grid1D::<IntervalClosed<Real>>::try_from_sorted(
                    SortedSet::from_unsorted(vec![
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                    ]),
                )
                .unwrap();
                let grid1d_b = Grid1D::<IntervalClosed<Real>>::try_from_sorted(
                    SortedSet::from_unsorted(vec![
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(0.5).unwrap(),
                        Real::try_new(2.0).unwrap(),
                    ]),
                )
                .unwrap();

                let union = Grid1DUnion::try_new(&grid1d_a, &grid1d_b).unwrap();
                let (grid1d_fine, map_a, map_b) = union.into_parts();
                assert_eq!(
                    grid1d_fine.coords().deref(),
                    &[
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(0.5).unwrap(),
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap()
                    ]
                );
                assert_eq!(map_a[0], IntervalId::new(0));
                assert_eq!(map_a[1], IntervalId::new(0));
                assert_eq!(map_a[2], IntervalId::new(1));

                assert_eq!(map_b[0], IntervalId::new(0));
                assert_eq!(map_b[1], IntervalId::new(1));
                assert_eq!(map_b[2], IntervalId::new(1));
            }

            #[test]
            fn intervals_in_intersection_single_interval() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                let domain =
                    IntervalClosed::new(Real::try_new(1.0).unwrap(), Real::try_new(1.5).unwrap());
                let intervals = grid.intervals_in_intersection(&domain);

                assert_eq!(intervals.len(), 1);
                assert_eq!(intervals[0].0, IntervalId::new(1));
                assert_eq!(
                    intervals[0].1,
                    IntervalFinitePositiveLength::from(IntervalClosed::new(
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(1.5).unwrap()
                    ))
                );
            }

            #[test]
            fn intervals_in_intersection_multiple_intervals() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                let domain =
                    IntervalClosed::new(Real::try_new(0.5).unwrap(), Real::try_new(2.5).unwrap());
                let intervals = grid.intervals_in_intersection(&domain);

                assert_eq!(intervals.len(), 3);

                // First interval
                assert_eq!(intervals[0].0, IntervalId::new(0));
                assert_eq!(
                    intervals[0].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(0.5).unwrap(),
                        Real::try_new(1.0).unwrap()
                    ))
                );

                // Middle interval
                assert_eq!(intervals[1].0, IntervalId::new(1));
                assert_eq!(
                    intervals[1].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap()
                    ))
                );

                // Last interval
                assert_eq!(intervals[2].0, IntervalId::new(2));
                assert_eq!(
                    intervals[2].1,
                    IntervalFinitePositiveLength::from(IntervalClosed::new(
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(2.5).unwrap()
                    ))
                );
            }

            #[test]
            fn intervals_in_intersection_full_overlap() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                let intervals = grid.intervals_in_intersection(grid.domain());

                assert_eq!(intervals.len(), 3);

                // First interval
                assert_eq!(intervals[0].0, IntervalId::new(0));
                assert_eq!(
                    intervals[0].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(1.0).unwrap()
                    ))
                );

                // Middle interval
                assert_eq!(intervals[1].0, IntervalId::new(1));
                assert_eq!(
                    intervals[1].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap()
                    ))
                );

                // Last interval
                assert_eq!(intervals[2].0, IntervalId::new(2));
                assert_eq!(
                    intervals[2].1,
                    IntervalFinitePositiveLength::from(IntervalClosed::new(
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(3.0).unwrap()
                    ))
                );
            }

            #[test]
            fn intervals_in_intersection_no_overlap() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                let domain =
                    IntervalClosed::new(Real::try_new(3.5).unwrap(), Real::try_new(4.0).unwrap());
                let intervals = grid.intervals_in_intersection(&domain);

                assert!(intervals.is_empty());
            }

            #[test]
            fn intervals_in_intersection_out_of_bounds_01() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                let domain =
                    IntervalClosed::new(Real::try_new(2.5).unwrap(), Real::try_new(4.0).unwrap());
                let intervals = grid.intervals_in_intersection(&domain);
                assert_eq!(
                    intervals[0],
                    (
                        IntervalId::new(2),
                        IntervalFinitePositiveLength::from(IntervalClosed::new(
                            Real::try_new(2.5).unwrap(),
                            Real::try_new(3.0).unwrap()
                        ))
                    )
                );
            }

            #[test]
            fn intervals_in_intersection_out_of_bounds_02() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                let domain = IntervalLowerOpenUpperClosed::new(
                    Real::try_new(2.5).unwrap(),
                    Real::try_new(4.0).unwrap(),
                );
                let intervals = grid.intervals_in_intersection(&domain);
                assert_eq!(
                    intervals[0],
                    (
                        IntervalId::new(2),
                        IntervalFinitePositiveLength::from(IntervalLowerOpenUpperClosed::new(
                            Real::try_new(2.5).unwrap(),
                            Real::try_new(3.0).unwrap()
                        ))
                    )
                );
            }

            #[test]
            fn test_from_grid1d_uniform_to_grid1d() {
                let domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let uniform_grid =
                    Grid1DUniform::new(domain.clone(), NumIntervals::try_new(4).unwrap());

                // Convert using From trait
                let grid: Grid1D<IntervalClosed<Real>> = Grid1D::from(uniform_grid.clone());

                // Verify it's a Uniform variant
                assert!(matches!(grid, Grid1D::Uniform(_)));

                // Verify properties are preserved
                assert_eq!(grid.domain(), &domain);
                assert_eq!(grid.num_intervals().as_ref(), &4);
                assert_eq!(grid.coords().deref(), uniform_grid.coords().deref());
            }

            #[test]
            fn test_from_grid1d_non_uniform_to_grid1d() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(0.3).unwrap(),
                    Real::try_new(0.7).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]))
                .unwrap();
                let non_uniform_grid =
                    Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords).unwrap();

                // Convert using From trait
                let grid: Grid1D<IntervalClosed<Real>> = Grid1D::from(non_uniform_grid.clone());

                // Verify it's a NonUniform variant
                assert!(matches!(grid, Grid1D::NonUniform(_)));

                // Verify properties are preserved
                assert_eq!(grid.domain(), non_uniform_grid.domain());
                assert_eq!(grid.num_intervals().as_ref(), &3);
                assert_eq!(grid.coords().deref(), non_uniform_grid.coords().deref());
            }

            #[test]
            fn test_from_grid1d_uniform_into_syntax() {
                let domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let uniform_grid =
                    Grid1DUniform::new(domain.clone(), NumIntervals::try_new(4).unwrap());

                // Test .into() syntax
                let grid: Grid1D<IntervalClosed<Real>> = uniform_grid.into();
                assert!(matches!(grid, Grid1D::Uniform(_)));
            }

            #[test]
            fn test_from_grid1d_non_uniform_into_syntax() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(0.5).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]))
                .unwrap();
                let non_uniform_grid =
                    Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords).unwrap();

                // Test .into() syntax
                let grid: Grid1D<IntervalClosed<Real>> = non_uniform_grid.into();
                assert!(matches!(grid, Grid1D::NonUniform(_)));
            }

            #[test]
            fn test_grid1d_non_uniform_refine_uniform_dispatch() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(0.2).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]))
                .unwrap();
                let grid = Grid1D::NonUniform(
                    Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords).unwrap(),
                );

                let refinement = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());

                assert!(matches!(refinement.original_grid(), Grid1D::NonUniform(_)));
                assert!(matches!(refinement.refined_grid(), Grid1D::NonUniform(_)));
                assert_eq!(refinement.original_grid().num_intervals().as_ref(), &2);
                assert_eq!(refinement.refined_grid().num_intervals().as_ref(), &4);

                assert_eq!(
                    refinement.find_original_interval(&IntervalId::new(0)),
                    IntervalId::new(0)
                );
                assert_eq!(
                    refinement.find_original_interval(&IntervalId::new(1)),
                    IntervalId::new(0)
                );
                assert_eq!(
                    refinement.find_original_interval(&IntervalId::new(2)),
                    IntervalId::new(1)
                );
                assert_eq!(
                    refinement.find_original_interval(&IntervalId::new(3)),
                    IntervalId::new(1)
                );
            }
        }

        mod grid1d_open_intervals {
            use super::*;

            #[test]
            fn test_open_interval_partition() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalOpen<Real>>::try_from_sorted(coords).unwrap();

                // Test first interval: should be (0, 1]
                let interval_0 = grid.interval(&IntervalId::new(0));
                match interval_0 {
                    SubIntervalInPartition::First(interval) => {
                        assert!(!interval.contains_point(&Real::try_new(0.0).unwrap())); // Open lower bound
                        assert!(interval.contains_point(&Real::try_new(0.5).unwrap())); // Inside
                        assert!(interval.contains_point(&Real::try_new(1.0).unwrap())); // Closed upper bound
                    }
                    _ => panic!("Expected first interval"),
                }

                // Test last interval: should be (1, 2)
                let interval_1 = grid.interval(&IntervalId::new(1));
                match interval_1 {
                    SubIntervalInPartition::Last(interval) => {
                        assert!(!interval.contains_point(&Real::try_new(1.0).unwrap())); // Open lower bound
                        assert!(interval.contains_point(&Real::try_new(1.5).unwrap())); // Inside
                        assert!(!interval.contains_point(&Real::try_new(2.0).unwrap())); // Open upper bound
                    }
                    _ => panic!("Expected last interval"),
                }
            }

            #[test]
            fn test_find_interval_on_open_domain() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalOpen<Real>>::try_from_sorted(coords).unwrap();

                // Points inside the domain should be found correctly
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.5).unwrap())
                        .unwrap(),
                    IntervalId::new(0)
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.5).unwrap())
                        .unwrap(),
                    IntervalId::new(1)
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(2.5).unwrap())
                        .unwrap(),
                    IntervalId::new(2)
                );

                // Point exactly on boundary
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.0).unwrap())
                        .unwrap(),
                    IntervalId::new(0)
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(2.0).unwrap())
                        .unwrap(),
                    IntervalId::new(1)
                );
            }

            #[test]
            fn test_find_interval_boundary_point_not_in_open_domain() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalOpen<Real>>::try_from_sorted(coords).unwrap();

                let res = grid.find_interval_id_of_point(&Real::try_new(0.0).unwrap());
                assert!(res.is_none(), "Expected None for point on open boundary");
            }
        }

        mod grid1d_semi_open_intervals {
            use super::*;

            #[test]
            fn test_lower_closed_upper_open_partition() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]);
                let grid =
                    Grid1D::<IntervalLowerClosedUpperOpen<Real>>::try_from_sorted(coords).unwrap();

                // Test first interval: should be [0, 1)
                let interval_0 = grid.interval(&IntervalId::new(0));
                match interval_0 {
                    SubIntervalInPartition::First(interval) => {
                        assert!(interval.contains_point(&Real::try_new(0.0).unwrap())); // Closed lower bound
                        assert!(interval.contains_point(&Real::try_new(0.5).unwrap())); // Inside
                        assert!(!interval.contains_point(&Real::try_new(1.0).unwrap())); // Open upper bound
                    }
                    _ => panic!("Expected first interval"),
                }

                // Test last interval: should be [1, 2)
                let interval_1 = grid.interval(&IntervalId::new(1));
                match interval_1 {
                    SubIntervalInPartition::Last(interval) => {
                        assert!(interval.contains_point(&Real::try_new(1.0).unwrap())); // Closed lower bound
                        assert!(interval.contains_point(&Real::try_new(1.5).unwrap())); // Inside
                        assert!(!interval.contains_point(&Real::try_new(2.0).unwrap())); // Open upper bound
                    }
                    _ => panic!("Expected last interval"),
                }
            }

            #[test]
            fn test_lower_open_upper_closed_partition() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]);
                let grid =
                    Grid1D::<IntervalLowerOpenUpperClosed<Real>>::try_from_sorted(coords).unwrap();

                // Test first interval: should be (0, 1]
                let interval_0 = grid.interval(&IntervalId::new(0));
                match interval_0 {
                    SubIntervalInPartition::First(interval) => {
                        assert!(!interval.contains_point(&Real::try_new(0.0).unwrap())); // Open lower bound
                        assert!(interval.contains_point(&Real::try_new(0.5).unwrap())); // Inside
                        assert!(interval.contains_point(&Real::try_new(1.0).unwrap())); // Closed upper bound
                    }
                    _ => panic!("Expected first interval"),
                }

                // Test last interval: should be (1, 2]
                let interval_1 = grid.interval(&IntervalId::new(1));
                match interval_1 {
                    SubIntervalInPartition::Last(interval) => {
                        assert!(!interval.contains_point(&Real::try_new(1.0).unwrap())); // Open lower bound
                        assert!(interval.contains_point(&Real::try_new(1.5).unwrap())); // Inside
                        assert!(interval.contains_point(&Real::try_new(2.0).unwrap())); // Closed upper bound
                    }
                    _ => panic!("Expected last interval"),
                }
            }
        }

        mod interval_partition_general {
            use super::*;

            #[test]
            fn test_interval_partition_general_creation() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(0.5).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]))
                .unwrap();
                let partition =
                    Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords.clone())
                        .unwrap();

                assert_eq!(partition.coords(), &coords);
                assert_eq!(
                    partition.domain(),
                    &IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap())
                );
            }

            #[test]
            fn test_interval_partition_general_insufficient_points() {
                let coords =
                    Coords1D::try_from(SortedSet::from_unsorted(vec![Real::try_new(0.0).unwrap()]))
                        .unwrap();
                let partition =
                    Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords);
                matches!(
                    partition,
                    Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { .. })
                );
            }
        }

        mod interval_partition_advanced {
            use super::*;

            #[test]
            fn test_interval_partition_iter_intervals() {
                let grid = Grid1D::uniform(
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(3.0).unwrap()),
                    NumIntervals::try_new(3).unwrap(),
                );

                let intervals: Vec<(IntervalId, SubIntervalInPartition<IntervalClosed<Real>>)> =
                    grid.iter_intervals().collect();
                assert_eq!(intervals.len(), 3);

                // Check that interval IDs are sequential
                for (i, (id, _)) in intervals.iter().enumerate() {
                    assert_eq!(*id.as_ref(), i);
                }
            }

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

                let points = vec![
                    Real::try_new(0.5).unwrap(),
                    Real::try_new(2.3).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(7.8).unwrap(),
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(11.0).unwrap(),
                ];
                let intervals = grid.find_intervals_for_points(&points);

                assert_eq!(intervals.len(), 6);
                assert_eq!(intervals[0], Some(IntervalId::new(0)));
                assert_eq!(intervals[1], Some(IntervalId::new(2)));
                assert_eq!(intervals[2], Some(IntervalId::new(5))); // p₅ = 5.0 → right interval
                assert_eq!(intervals[3], Some(IntervalId::new(7)));
                assert_eq!(intervals[4], Some(IntervalId::new(9)));
                assert_eq!(intervals[5], None); // Outside domain
            }

            #[test]
            fn test_interval_partition_length_statistics() {
                let coords = SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(1.5).unwrap(),
                    Real::try_new(4.0).unwrap(),
                ]);
                let grid = Grid1D::<IntervalClosed<Real>>::try_from_sorted(coords).unwrap();

                // Intervals: [0,1] (length 1), (1,1.5] (length 0.5), (1.5,4] (length 2.5)
                assert_eq!(
                    grid.min_interval_length().as_ref(),
                    &Real::try_new(0.5).unwrap()
                );
                assert_eq!(
                    grid.max_interval_length().as_ref(),
                    &Real::try_new(2.5).unwrap()
                );
                assert_eq!(
                    grid.uniformity_ratio().as_ref(),
                    &Real::try_new(5.0).unwrap()
                ); // 2.5/0.5
            }

            #[test]
            fn test_interval_partition_intersection_edge_cases() {
                let grid = Grid1D::uniform(
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(2.0).unwrap()),
                    NumIntervals::try_new(2).unwrap(),
                );

                // Test intersection with smaller domain
                let domain_small =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(0.5).unwrap());
                let intersections = grid.intervals_in_intersection(&domain_small);
                assert_eq!(
                    &intersections,
                    &[(
                        IntervalId::new(0),
                        IntervalClosed::new(
                            Real::try_new(0.0).unwrap(),
                            Real::try_new(0.5).unwrap()
                        )
                        .into()
                    )]
                );

                // Test intersection with entire domain
                let domain_full =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(2.0).unwrap());
                let intersections_full = grid.intervals_in_intersection(&domain_full);
                assert_eq!(
                    &intersections_full,
                    &[
                        (
                            IntervalId::new(0),
                            IntervalLowerClosedUpperOpen::new(
                                Real::try_new(0.0).unwrap(),
                                Real::try_new(1.0).unwrap()
                            )
                            .into()
                        ),
                        (
                            IntervalId::new(1),
                            IntervalClosed::new(
                                Real::try_new(1.0).unwrap(),
                                Real::try_new(2.0).unwrap()
                            )
                            .into()
                        )
                    ]
                );

                // Test intersection with domain slightly larger
                let domain_larger =
                    IntervalClosed::new(Real::try_new(-0.1).unwrap(), Real::try_new(2.1).unwrap());
                let intersections_larger = grid.intervals_in_intersection(&domain_larger);
                assert_eq!(intersections_larger.len(), 2); // Should clip to actual domain
                assert_eq!(
                    &intersections_larger,
                    &[
                        (
                            IntervalId::new(0),
                            IntervalLowerClosedUpperOpen::new(
                                Real::try_new(0.0).unwrap(),
                                Real::try_new(1.0).unwrap()
                            )
                            .into()
                        ),
                        (
                            IntervalId::new(1),
                            IntervalClosed::new(
                                Real::try_new(1.0).unwrap(),
                                Real::try_new(2.0).unwrap()
                            )
                            .into()
                        )
                    ]
                );
            }
        }

        mod grid1d_uniform {
            use super::*;

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

                assert_eq!(grid.domain(), &domain);
                assert_eq!(grid.num_intervals(), num_intervals);
                assert_eq!(
                    grid.coords().deref(),
                    &[
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(4.0).unwrap(),
                        Real::try_new(6.0).unwrap(),
                        Real::try_new(8.0).unwrap(),
                        Real::try_new(10.0).unwrap()
                    ]
                );
                assert_eq!(*grid.delta_points().as_ref(), Real::try_new(2.0).unwrap());
            }

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

                // Test interval lengths and uniformity
                assert_eq!(
                    *grid.interval_length(&IntervalId::new(0)).as_ref(),
                    Real::try_new(1.0).unwrap()
                );
                assert_eq!(
                    *grid.min_interval_length().as_ref(),
                    Real::try_new(1.0).unwrap()
                );
                assert_eq!(
                    *grid.max_interval_length().as_ref(),
                    Real::try_new(1.0).unwrap()
                );
                assert_eq!(
                    *grid.uniformity_ratio().as_ref(),
                    Real::try_new(1.0).unwrap()
                );

                // Test interval types
                let interval_0 = grid.interval(&IntervalId::new(0));
                let interval_1 = grid.interval(&IntervalId::new(1));

                // First interval: [0, 1)
                match interval_0 {
                    SubIntervalInPartition::First(i) => {
                        assert_eq!(
                            i,
                            IntervalLowerClosedUpperOpen::new(
                                Real::try_new(0.0).unwrap(),
                                Real::try_new(1.0).unwrap()
                            )
                        );
                        assert!(i.contains_point(&Real::try_new(0.0).unwrap()));
                        assert!(!i.contains_point(&Real::try_new(1.0).unwrap()));
                    }
                    _ => panic!("Expected first interval"),
                }

                // Last interval: [1, 2]
                match interval_1 {
                    SubIntervalInPartition::Last(i) => {
                        assert_eq!(
                            i,
                            IntervalClosed::new(
                                Real::try_new(1.0).unwrap(),
                                Real::try_new(2.0).unwrap()
                            )
                        );
                        assert!(i.contains_point(&Real::try_new(1.0).unwrap()));
                        assert!(i.contains_point(&Real::try_new(2.0).unwrap()));
                    }
                    _ => panic!("Expected last interval"),
                }
            }

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

                // Test boundaries
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.0).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(10.0).unwrap()),
                    Some(IntervalId::new(9))
                );

                // Test interior points
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.5).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.5).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(9.9).unwrap()),
                    Some(IntervalId::new(9))
                );

                // Test points exactly on grid nodes (belong to RIGHT interval for [a,b] domain)
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(1))
                ); // p₁ → right interval
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(2.0).unwrap()),
                    Some(IntervalId::new(2))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(9.0).unwrap()),
                    Some(IntervalId::new(9))
                );

                // Test points outside domain
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(-0.1).unwrap()),
                    None
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(10.1).unwrap()),
                    None
                );
            }

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

                // Point near start
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.00005).unwrap()),
                    Some(IntervalId::new(0))
                );
                // Point on first node (p₁ = 0.0001 belongs to right interval for left-closed)
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.0001).unwrap()),
                    Some(IntervalId::new(1))
                );
                // Point just after first node
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.00011).unwrap()),
                    Some(IntervalId::new(1))
                );
                // Point in the middle
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.54321).unwrap()),
                    Some(IntervalId::new(5432))
                );
                // Point near end
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.99999).unwrap()),
                    Some(IntervalId::new(9999))
                );
                // Point on end
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(9999))
                );
            }

            #[test]
            fn test_uniform_grid_open_domain() {
                let domain =
                    IntervalOpen::new(Real::try_new(0.0).unwrap(), Real::try_new(2.0).unwrap());
                let grid = Grid1DUniform::new(domain, NumIntervals::try_new(2).unwrap());

                let expected_coords = [
                    Real::try_new(0.).unwrap(),
                    Real::try_new(1.).unwrap(),
                    Real::try_new(2.).unwrap(),
                ];
                assert_eq!(grid.coords().deref(), expected_coords.as_slice());

                // First interval: (0, 1]
                let interval_0 = grid.interval(&IntervalId::new(0));
                match interval_0 {
                    SubIntervalInPartition::First(i) => {
                        assert_eq!(
                            i,
                            IntervalLowerOpenUpperClosed::new(
                                Real::try_new(0.0).unwrap(),
                                Real::try_new(1.0).unwrap()
                            )
                        );
                        assert!(!i.contains_point(&Real::try_new(0.0).unwrap()));
                        assert!(i.contains_point(&Real::try_new(1.0).unwrap()));
                    }
                    _ => panic!("Expected first interval"),
                }

                // Last interval: (1, 2)
                let interval_1 = grid.interval(&IntervalId::new(1));
                match interval_1 {
                    SubIntervalInPartition::Last(i) => {
                        assert_eq!(
                            i,
                            IntervalOpen::new(
                                Real::try_new(1.0).unwrap(),
                                Real::try_new(2.0).unwrap()
                            )
                        );
                        assert!(!i.contains_point(&Real::try_new(1.0).unwrap()));
                        assert!(!i.contains_point(&Real::try_new(2.0).unwrap()));
                    }
                    _ => panic!("Expected last interval"),
                }

                // Point location
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.1).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.9).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.0).unwrap()),
                    None
                ); // Outside open domain
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(2.0).unwrap()),
                    None
                ); // Outside open domain
            }

            #[test]
            fn test_uniform_grid_intersection() {
                let grid = Grid1DUniform::new(
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(10.0).unwrap()),
                    NumIntervals::try_new(10).unwrap(),
                );
                let domain_in =
                    IntervalClosed::new(Real::try_new(2.5).unwrap(), Real::try_new(7.5).unwrap());
                let intersections = grid.intervals_in_intersection(&domain_in);

                assert_eq!(intersections.len(), 6);

                // Check first intersection: interval 2, [2.0, 3.0) intersects with [2.5, 7.5] -> [2.5, 3.0)
                assert_eq!(intersections[0].0, IntervalId::new(2));
                assert_eq!(
                    intersections[0].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(2.5).unwrap(),
                        Real::try_new(3.0).unwrap()
                    ))
                );

                // Check a middle intersection: interval 4, [4.0, 5.0) intersects with [2.5, 7.5] -> [4.0, 5.0)
                assert_eq!(intersections[2].0, IntervalId::new(4));
                assert_eq!(
                    intersections[2].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(4.0).unwrap(),
                        Real::try_new(5.0).unwrap()
                    ))
                );

                // Check last intersection: interval 7, [7.0, 8.0] intersects with [2.5, 7.5] -> [7.0, 7.5]
                assert_eq!(intersections[5].0, IntervalId::new(7));
                assert_eq!(
                    intersections[5].1,
                    IntervalFinitePositiveLength::from(IntervalClosed::new(
                        Real::try_new(7.0).unwrap(),
                        Real::try_new(7.5).unwrap()
                    ))
                );
            }
        }

        mod grid1d_non_uniform {
            use super::*;
            use crate::intervals::IntervalClosed;

            fn create_test_grid() -> Grid1DNonUniform<IntervalClosed<Real>> {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(7.0).unwrap(),
                    Real::try_new(10.0).unwrap(),
                ]))
                .unwrap();
                Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords).unwrap()
            }

            #[test]
            fn test_non_uniform_grid_creation() {
                let grid = create_test_grid();
                assert_eq!(
                    grid.domain(),
                    &IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(10.0).unwrap())
                );
                assert_eq!(
                    grid.coords().deref(),
                    &[
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(3.0).unwrap(),
                        Real::try_new(7.0).unwrap(),
                        Real::try_new(10.0).unwrap()
                    ]
                );
                assert_eq!(grid.num_intervals().as_ref(), &4);
            }

            #[test]
            fn test_creation_fails_with_insufficient_points() {
                let coords =
                    Coords1D::try_from(SortedSet::from_unsorted(vec![Real::try_new(0.0).unwrap()]))
                        .unwrap();
                let result = Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords);
                assert!(matches!(
                    result,
                    Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { .. })
                ));
            }

            #[test]
            fn test_non_uniform_interval_partition() {
                let grid = create_test_grid();

                // First interval: [0, 1)
                let interval_0 = grid.interval(&IntervalId::new(0));
                match interval_0 {
                    SubIntervalInPartition::First(i) => {
                        assert_eq!(
                            i,
                            IntervalLowerClosedUpperOpen::new(
                                Real::try_new(0.0).unwrap(),
                                Real::try_new(1.0).unwrap()
                            )
                        );
                    }
                    _ => panic!("Expected first interval"),
                }

                // Middle interval: [1, 3)
                let interval_1 = grid.interval(&IntervalId::new(1));
                match interval_1 {
                    SubIntervalInPartition::Middle(i) => {
                        assert_eq!(
                            i,
                            IntervalLowerClosedUpperOpen::new(
                                Real::try_new(1.0).unwrap(),
                                Real::try_new(3.0).unwrap()
                            )
                        );
                    }
                    _ => panic!("Expected middle interval"),
                }

                // Last interval: [7, 10]
                let interval_3 = grid.interval(&IntervalId::new(3));
                match interval_3 {
                    SubIntervalInPartition::Last(i) => {
                        assert_eq!(
                            i,
                            IntervalClosed::new(
                                Real::try_new(7.0).unwrap(),
                                Real::try_new(10.0).unwrap()
                            )
                        );
                    }
                    _ => panic!("Expected last interval"),
                }
            }

            #[test]
            fn test_non_uniform_statistics() {
                let grid = create_test_grid();
                assert_eq!(
                    *grid.interval_length(&IntervalId::new(0)).as_ref(),
                    Real::try_new(1.0).unwrap()
                ); // [0,1]
                assert_eq!(
                    *grid.interval_length(&IntervalId::new(1)).as_ref(),
                    Real::try_new(2.0).unwrap()
                ); // (1,3]
                assert_eq!(
                    *grid.interval_length(&IntervalId::new(2)).as_ref(),
                    Real::try_new(4.0).unwrap()
                ); // (3,7]
                assert_eq!(
                    *grid.interval_length(&IntervalId::new(3)).as_ref(),
                    Real::try_new(3.0).unwrap()
                ); // (7,10]

                assert_eq!(
                    *grid.min_interval_length().as_ref(),
                    Real::try_new(1.0).unwrap()
                );
                assert_eq!(
                    *grid.max_interval_length().as_ref(),
                    Real::try_new(4.0).unwrap()
                );
                assert_eq!(
                    *grid.uniformity_ratio().as_ref(),
                    Real::try_new(4.0).unwrap()
                );
            }

            #[test]
            fn test_non_uniform_point_location() {
                let grid = create_test_grid();

                // Boundaries
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.0).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(10.0).unwrap()),
                    Some(IntervalId::new(3))
                );

                // Interior points
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(0.5).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(2.5).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(5.0).unwrap()),
                    Some(IntervalId::new(2))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(8.5).unwrap()),
                    Some(IntervalId::new(3))
                );

                // Points exactly on grid nodes (belong to RIGHT interval for [a,b] domain)
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(3.0).unwrap()),
                    Some(IntervalId::new(2))
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(7.0).unwrap()),
                    Some(IntervalId::new(3))
                );

                // Outside domain
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(-0.1).unwrap()),
                    None
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(10.1).unwrap()),
                    None
                );
            }

            #[test]
            fn test_non_uniform_intersection() {
                let grid = create_test_grid();
                let domain_in =
                    IntervalClosed::new(Real::try_new(0.5).unwrap(), Real::try_new(8.5).unwrap());
                let intersections = grid.intervals_in_intersection(&domain_in);

                assert_eq!(intersections.len(), 4);

                // First intersection: interval 0, [0,1) intersects [0.5, 8.5] -> [0.5, 1)
                assert_eq!(intersections[0].0, IntervalId::new(0));
                assert_eq!(
                    intersections[0].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(0.5).unwrap(),
                        Real::try_new(1.0).unwrap()
                    ))
                );

                // Second intersection: interval 1, [1,3) intersects [0.5, 8.5] -> [1,3)
                assert_eq!(intersections[1].0, IntervalId::new(1));
                assert_eq!(
                    intersections[1].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(3.0).unwrap()
                    ))
                );

                // Third intersection: interval 2, [3,7) intersects [0.5, 8.5] -> [3,7)
                assert_eq!(intersections[2].0, IntervalId::new(2));
                assert_eq!(
                    intersections[2].1,
                    IntervalFinitePositiveLength::from(IntervalLowerClosedUpperOpen::new(
                        Real::try_new(3.0).unwrap(),
                        Real::try_new(7.0).unwrap()
                    ))
                );

                // Fourth intersection: interval 3, [7,10] intersects [0.5, 8.5] -> [7,8.5]
                assert_eq!(intersections[3].0, IntervalId::new(3));
                assert_eq!(
                    intersections[3].1,
                    IntervalFinitePositiveLength::from(IntervalClosed::new(
                        Real::try_new(7.0).unwrap(),
                        Real::try_new(8.5).unwrap()
                    ))
                );
            }
        }

        // NOTE: Comprehensive tests for Grid1DRefinement are in src/operations/refinement.rs
        // NOTE: Comprehensive tests for Grid1DUnion are in src/operations/union.rs

        mod find_interval_id_of_point {
            use super::*;

            fn r(v: f64) -> Real {
                Real::try_new(v).unwrap()
            }

            /// Checks that `find_interval_id_of_point` returns `expected` AND,
            /// when the result is `Some(id)`, that `grid.interval(&id)` actually
            /// contains the queried point.
            fn assert_finds<G>(grid: &G, x: G::CoordType, expected: Option<IntervalId>)
            where
                G: Grid1DTrait,
            {
                let found = grid.find_interval_id_of_point(&x);
                assert_eq!(found, expected);
                if let Some(ref id) = found {
                    assert!(
                        grid.interval(id).contains_point(&x),
                        "interval {:?} returned for point but does not contain it",
                        id
                    );
                }
            }

            // ----------------------------------------------------------------
            // Grid1DNonUniform tests
            //
            // Coords: [0.0, 1.0, 3.0, 6.0, 10.0]  → 4 intervals
            // Internal boundaries: p₁=1.0, p₂=3.0, p₃=6.0
            // ----------------------------------------------------------------
            mod grid1d_non_uniform {
                use super::*;

                mod interval_closed {
                    use super::*;

                    // Domain [0, 10]: sub-intervals [0,1), [1,3), [3,6), [6,10]
                    // Lower-bound closed → internal boundary pₖ belongs to RIGHT interval k

                    fn make_grid() -> Grid1DNonUniform<IntervalClosed<Real>> {
                        let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                            r(0.0),
                            r(1.0),
                            r(3.0),
                            r(6.0),
                            r(10.0),
                        ]))
                        .unwrap();
                        Grid1DNonUniform::<IntervalClosed<Real>>::try_new_from_coords(coords)
                            .unwrap()
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(4.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(8.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary() {
                        // p₀ = 0.0, closed lower bound → belongs to interval 0
                        assert_finds(&make_grid(), r(0.0), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn internal_boundaries_belong_to_right_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(2)));
                        assert_finds(&grid, r(6.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn upper_domain_boundary() {
                        // p₄ = 10.0, closed upper bound → belongs to last interval
                        assert_finds(&make_grid(), r(10.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(10.1), None);
                    }
                }

                mod interval_open {
                    use super::*;

                    // Domain (0, 10): sub-intervals (0,1], (1,3], (3,6], (6,10)
                    // Lower-bound open → internal boundary pₖ belongs to LEFT interval k-1

                    fn make_grid() -> Grid1DNonUniform<IntervalOpen<Real>> {
                        let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                            r(0.0),
                            r(1.0),
                            r(3.0),
                            r(6.0),
                            r(10.0),
                        ]))
                        .unwrap();
                        Grid1DNonUniform::<IntervalOpen<Real>>::try_new_from_coords(coords).unwrap()
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(4.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(8.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary_returns_none() {
                        // p₀ = 0.0 is NOT in open domain (0, 10)
                        assert_finds(&make_grid(), r(0.0), None);
                    }

                    #[test]
                    fn internal_boundaries_belong_to_left_interval() {
                        // (a,b): sub-intervals are right-closed (pₖ, pₖ₊₁]
                        // pₖ belongs to the left interval k-1
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(0)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(6.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn upper_domain_boundary_returns_none() {
                        // p₄ = 10.0 is NOT in open domain (0, 10)
                        assert_finds(&make_grid(), r(10.0), None);
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(10.1), None);
                    }
                }

                mod interval_lower_closed_upper_open {
                    use super::*;

                    // Domain [0, 10): sub-intervals [0,1), [1,3), [3,6), [6,10)
                    // Lower-bound closed → internal pₖ belongs to RIGHT interval k

                    fn make_grid() -> Grid1DNonUniform<IntervalLowerClosedUpperOpen<Real>> {
                        let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                            r(0.0),
                            r(1.0),
                            r(3.0),
                            r(6.0),
                            r(10.0),
                        ]))
                        .unwrap();
                        Grid1DNonUniform::<IntervalLowerClosedUpperOpen<Real>>::try_new_from_coords(
                            coords,
                        )
                        .unwrap()
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(4.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(8.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary() {
                        // p₀ = 0.0 is in [a, b) domain
                        assert_finds(&make_grid(), r(0.0), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn internal_boundaries_belong_to_right_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(2)));
                        assert_finds(&grid, r(6.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn upper_domain_boundary_returns_none() {
                        // p₄ = 10.0 is NOT in [a, b) domain
                        assert_finds(&make_grid(), r(10.0), None);
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(10.1), None);
                    }
                }

                mod interval_lower_open_upper_closed {
                    use super::*;

                    // Domain (0, 10]: sub-intervals (0,1], (1,3], (3,6], (6,10]
                    // Lower-bound open → internal pₖ belongs to LEFT interval k-1

                    fn make_grid() -> Grid1DNonUniform<IntervalLowerOpenUpperClosed<Real>> {
                        let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                            r(0.0),
                            r(1.0),
                            r(3.0),
                            r(6.0),
                            r(10.0),
                        ]))
                        .unwrap();
                        Grid1DNonUniform::<IntervalLowerOpenUpperClosed<Real>>::try_new_from_coords(
                            coords,
                        )
                        .unwrap()
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(4.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(8.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary_returns_none() {
                        // p₀ = 0.0 is NOT in (a, b] domain
                        assert_finds(&make_grid(), r(0.0), None);
                    }

                    #[test]
                    fn internal_boundaries_belong_to_left_interval() {
                        // (a,b]: sub-intervals are right-closed (pₖ, pₖ₊₁]
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(0)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(6.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn upper_domain_boundary() {
                        // p₄ = 10.0 is in (a, b] domain → belongs to last interval
                        assert_finds(&make_grid(), r(10.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(10.1), None);
                    }
                }

                mod interval_finite_positive_length {
                    use super::*;

                    // IntervalFinitePositiveLength::Closed creates a Closed variant,
                    // so boundary semantics are identical to the interval_closed tests.
                    // Domain [0, 10]: sub-intervals [0,1), [1,3), [3,6), [6,10]

                    fn make_grid() -> Grid1DNonUniform<IntervalFinitePositiveLength<Real>> {
                        let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                            r(0.0),
                            r(1.0),
                            r(3.0),
                            r(6.0),
                            r(10.0),
                        ]))
                        .unwrap();
                        let domain: IntervalFinitePositiveLength<Real> =
                            IntervalClosed::new(r(0.0), r(10.0)).into();

                        Grid1DNonUniform::<IntervalFinitePositiveLength<Real>>::try_new(
                            domain, coords,
                        )
                        .unwrap()
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(4.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(8.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary() {
                        assert_finds(&make_grid(), r(0.0), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn internal_boundaries_belong_to_right_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(2)));
                        assert_finds(&grid, r(6.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn upper_domain_boundary() {
                        assert_finds(&make_grid(), r(10.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(10.1), None);
                    }
                }
            }

            // ----------------------------------------------------------------
            // Grid1DUniform tests
            //
            // Domain [0.0, 4.0] with 4 intervals → coords [0.0, 1.0, 2.0, 3.0, 4.0]
            // Internal boundaries: p₁=1.0, p₂=2.0, p₃=3.0
            // ----------------------------------------------------------------
            mod grid1d_uniform {
                use super::*;

                mod interval_closed {
                    use super::*;

                    // Domain [0, 4]: sub-intervals [0,1), [1,2), [2,3), [3,4]
                    // Lower-bound closed → internal pₖ belongs to RIGHT interval k

                    fn make_grid() -> Grid1DUniform<IntervalClosed<Real>> {
                        Grid1DUniform::new(
                            IntervalClosed::new(r(0.0), r(4.0)),
                            NumIntervals::try_new(4).unwrap(),
                        )
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.5), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.5), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(3.5), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary() {
                        assert_finds(&make_grid(), r(0.0), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn internal_boundaries_belong_to_right_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(2)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn upper_domain_boundary() {
                        assert_finds(&make_grid(), r(4.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(4.1), None);
                    }
                }

                mod interval_open {
                    use super::*;

                    // Domain (0, 4): sub-intervals (0,1], (1,2], (2,3], (3,4)
                    // Lower-bound open → internal pₖ belongs to LEFT interval k-1

                    fn make_grid() -> Grid1DUniform<IntervalOpen<Real>> {
                        Grid1DUniform::new(
                            IntervalOpen::new(r(0.0), r(4.0)),
                            NumIntervals::try_new(4).unwrap(),
                        )
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.5), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.5), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(3.5), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary_returns_none() {
                        // p₀ = 0.0 is NOT in (a, b) domain
                        assert_finds(&make_grid(), r(0.0), None);
                    }

                    #[test]
                    fn internal_boundaries_belong_to_left_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(0)));
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn upper_domain_boundary_returns_none() {
                        // p₄ = 4.0 is NOT in (a, b) domain
                        assert_finds(&make_grid(), r(4.0), None);
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(4.1), None);
                    }
                }

                mod interval_lower_closed_upper_open {
                    use super::*;

                    // Domain [0, 4): sub-intervals [0,1), [1,2), [2,3), [3,4)
                    // Lower-bound closed → internal pₖ belongs to RIGHT interval k

                    fn make_grid() -> Grid1DUniform<IntervalLowerClosedUpperOpen<Real>> {
                        Grid1DUniform::new(
                            IntervalLowerClosedUpperOpen::new(r(0.0), r(4.0)),
                            NumIntervals::try_new(4).unwrap(),
                        )
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.5), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.5), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(3.5), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary() {
                        // p₀ = 0.0 is in [a, b) domain
                        assert_finds(&make_grid(), r(0.0), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn internal_boundaries_belong_to_right_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(2)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn upper_domain_boundary_returns_none() {
                        // p₄ = 4.0 is NOT in [a, b) domain
                        assert_finds(&make_grid(), r(4.0), None);
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(4.1), None);
                    }
                }

                mod interval_lower_open_upper_closed {
                    use super::*;

                    // Domain (0, 4]: sub-intervals (0,1], (1,2], (2,3], (3,4]
                    // Lower-bound open → internal pₖ belongs to LEFT interval k-1

                    fn make_grid() -> Grid1DUniform<IntervalLowerOpenUpperClosed<Real>> {
                        Grid1DUniform::new(
                            IntervalLowerOpenUpperClosed::new(r(0.0), r(4.0)),
                            NumIntervals::try_new(4).unwrap(),
                        )
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.5), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.5), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(3.5), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary_returns_none() {
                        // p₀ = 0.0 is NOT in (a, b] domain
                        assert_finds(&make_grid(), r(0.0), None);
                    }

                    #[test]
                    fn internal_boundaries_belong_to_left_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(0)));
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn upper_domain_boundary() {
                        // p₄ = 4.0 is in (a, b] domain → belongs to last interval
                        assert_finds(&make_grid(), r(4.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(4.1), None);
                    }
                }

                mod interval_finite_positive_length {
                    use super::*;

                    // IntervalFinitePositiveLength::Closed creates a Closed variant,
                    // so boundary semantics are identical to interval_closed tests.
                    // Domain [0, 4]: sub-intervals [0,1), [1,2), [2,3), [3,4]

                    fn make_grid() -> Grid1DUniform<IntervalFinitePositiveLength<Real>> {
                        Grid1DUniform::new(
                            IntervalFinitePositiveLength::Closed(IntervalClosed::new(
                                r(0.0),
                                r(4.0),
                            )),
                            NumIntervals::try_new(4).unwrap(),
                        )
                    }

                    #[test]
                    fn first_interval_interior() {
                        assert_finds(&make_grid(), r(0.5), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn middle_intervals_interior() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.5), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.5), Some(IntervalId::new(2)));
                    }

                    #[test]
                    fn last_interval_interior() {
                        assert_finds(&make_grid(), r(3.5), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn lower_domain_boundary() {
                        assert_finds(&make_grid(), r(0.0), Some(IntervalId::new(0)));
                    }

                    #[test]
                    fn internal_boundaries_belong_to_right_interval() {
                        let grid = make_grid();
                        assert_finds(&grid, r(1.0), Some(IntervalId::new(1)));
                        assert_finds(&grid, r(2.0), Some(IntervalId::new(2)));
                        assert_finds(&grid, r(3.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn upper_domain_boundary() {
                        assert_finds(&make_grid(), r(4.0), Some(IntervalId::new(3)));
                    }

                    #[test]
                    fn outside_domain_returns_none() {
                        let grid = make_grid();
                        assert_finds(&grid, r(-0.1), None);
                        assert_finds(&grid, r(4.1), None);
                    }
                }
            }
        }
    }

    #[cfg(feature = "rug")]
    mod rug100 {
        use super::*;
        use num::{One, Zero};
        use num_valid::RealRugStrictFinite;

        const PRECISION: u32 = 100;
        type Real = RealRugStrictFinite<PRECISION>;

        mod grid1d {
            use super::*;
            use num_valid::Constants;

            #[test]
            fn grid1d_01() {
                let v = SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0., 1.]));
                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();
                assert_eq!(grid1d.coords().deref(), &v.to_vec());
            }

            #[test]
            #[cfg(debug_assertions)]
            #[should_panic]
            fn grid1d_02() {
                let v = SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0.]));
                let _grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();
            }

            #[test]
            fn interval_01() {
                let v = SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0., 1.]));

                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();

                let interval = grid1d.interval(&IntervalId::new(0));

                match interval {
                    SubIntervalInPartition::Single(interval) => {
                        assert_eq!(
                            TryInto::<IntervalClosed<Real>>::try_into(interval).unwrap(),
                            IntervalClosed::new(Real::zero(), Real::one())
                        );
                    }
                    _ => {
                        panic!("There should be only one interval")
                    }
                }
            }

            #[test]
            #[cfg(debug_assertions)]
            #[should_panic]
            fn interval_02() {
                let v = SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0., 1.]));
                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();
                let _interval = grid1d.interval(&IntervalId::new(1));
            }

            #[test]
            fn find_interval_id_01() {
                let v =
                    SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0., 1., 2., 3.]));
                let grid1d = Grid1D::<IntervalClosed<Real>>::try_from_sorted(v.clone()).unwrap();

                let i0 = grid1d.find_interval_id_of_point(&Real::zero()).unwrap();
                assert_eq!(i0, IntervalId::new(0));

                let i1 = grid1d
                    .find_interval_id_of_point(
                        &Real::try_new(rug::Float::with_val(PRECISION, 0.5)).unwrap(),
                    )
                    .unwrap();
                assert_eq!(i1, IntervalId::new(0));

                let i2 = grid1d.find_interval_id_of_point(&Real::one()).unwrap();
                assert_eq!(i2, IntervalId::new(1));

                let i3 = grid1d
                    .find_interval_id_of_point(
                        &Real::try_new(rug::Float::with_val(PRECISION, 1.5)).unwrap(),
                    )
                    .unwrap();
                assert_eq!(i3, IntervalId::new(1));

                let i4 = grid1d.find_interval_id_of_point(&Real::two()).unwrap();
                assert_eq!(i4, IntervalId::new(2));

                let i5 = grid1d
                    .find_interval_id_of_point(
                        &Real::try_new(rug::Float::with_val(PRECISION, 2.5)).unwrap(),
                    )
                    .unwrap();
                assert_eq!(i5, IntervalId::new(2));

                let i6 = grid1d
                    .find_interval_id_of_point(
                        &Real::try_new(rug::Float::with_val(PRECISION, 3.)).unwrap(),
                    )
                    .unwrap();
                assert_eq!(i6, IntervalId::new(2));
            }

            #[test]
            fn uniform_01() {
                let interval = IntervalClosed::new(Real::zero(), Real::one());
                let grid1d = Grid1D::uniform(interval, NumIntervals::try_new(4).unwrap());

                let expected_coords = vec_f64_into_vec_real::<Real>(vec![0., 0.25, 0.5, 0.75, 1.0]);

                assert_eq!(grid1d.coords().deref(), &expected_coords);
            }

            #[test]
            fn intersection_01() {
                let grid1d_a = Grid1D::<IntervalClosed<Real>>::try_from_sorted(
                    SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0.0, 1.0, 2.0])),
                )
                .unwrap();
                let grid1d_b = Grid1D::<IntervalClosed<Real>>::try_from_sorted(
                    SortedSet::from_unsorted(vec_f64_into_vec_real::<Real>(vec![0.0, 0.5, 2.0])),
                )
                .unwrap();

                let union = Grid1DUnion::try_new(&grid1d_a, &grid1d_b).unwrap();
                let (grid1d_fine, map_a, map_b) = union.into_parts();
                assert_eq!(
                    grid1d_fine.coords().deref(),
                    &vec_f64_into_vec_real::<f64>(vec![0.0, 0.5, 1.0, 2.0])
                );
                assert_eq!(map_a[0], IntervalId::new(0));
                assert_eq!(map_a[1], IntervalId::new(0));
                assert_eq!(map_a[2], IntervalId::new(1));

                assert_eq!(map_b[0], IntervalId::new(0));
                assert_eq!(map_b[1], IntervalId::new(1));
                assert_eq!(map_b[2], IntervalId::new(1));
            }
        }
    }
}
//------------------------------------------------------------------------------------------------