grid1d 0.4.0

#![deny(rustdoc::broken_intra_doc_links)]
#![feature(error_generic_member_access)]

//! One-dimensional grids and interval partitions with mathematical rigor and performance optimization.
//!
//! This library provides a comprehensive framework for representing and manipulating
//! ***one-dimensional grids*** and ***interval partitions*** with strong mathematical
//! foundations, type safety, and performance optimizations. It serves as the backbone
//! for numerical computations, finite element methods, adaptive mesh refinement, and
//! scientific computing applications requiring robust 1D spatial discretizations.
//!
//! ## Core Design Philosophy
//!
//! ### Mathematical Correctness First
//! The library prioritizes **mathematical rigor** over convenience:
//! - **Type-safe intervals**: Different interval types ([`intervals::IntervalClosed`], [`intervals::IntervalOpen`], etc.) encode boundary semantics at compile time
//! - **Domain validation**: Runtime checks ensure grids completely partition their specified domains
//! - **Proper boundary handling**: Open/closed boundary semantics are correctly preserved in all operations
//! - **Generic scalar support**: Works with any scalar type implementing [`num_valid::RealScalar`]
//!
//! ### Multiple Floating-Point Backends via num-valid
//!
//! **grid1d** leverages the [`num-valid`](https://crates.io/crates/num-valid) library to provide **multiple numerical backends**
//! with different precision and performance characteristics:
//!
//! #### **Native Backend (Default)**
//! - **Base types**: [`f64`] for maximum performance with IEEE 754 compliance
//! - **Validated types**: [`num_valid::RealNative64StrictFinite`] for validation in Debug and Release mode,
//!   [`num_valid::RealNative64StrictFiniteInDebug`] for validation only in Debug mode
//! - **Performance**: ⚡⚡⚡ Maximum speed for [`num_valid::RealNative64StrictFiniteInDebug`] (same as `f64`) with O(1) operations
//! - **Precision**: Standard IEEE 754 double precision (53-bit mantissa)
//! - **Best for**: High-performance simulations, real-time applications
//!
//! #### **Arbitrary Precision Backend (Optional, activated with `--features=rug`)**
//! - **Base types**: [`rug::Float`](https://docs.rs/rug/latest/rug/float/struct.Float.html), [`rug::Complex`](https://docs.rs/rug/latest/rug/complex/struct.Complex.html) from the GNU Multiple Precision Arithmetic Library
//! - **Validated types**: `num_valid::RealRugStrictFinite<N>`
//! - **Performance**: ⚡ Configurable, depends on precision setting
//! - **Precision**: User-defined N-bit precision (e.g., 100, 256, 1024 bits)
//! - **Best for**: Scientific computing requiring exact arithmetic, symbolic computation
//!
//! **All grid operations work seamlessly across backends** - you can switch from standard
//! to arbitrary precision by simply changing the scalar type, with no changes to your
//! grid manipulation code.
//!
//! ```rust
//! use grid1d::{Grid1D, IntervalPartition, intervals::*, scalars::NumIntervals};
//! use num_valid::{RealNative64StrictFinite, RealScalar};
//! use try_create::TryNew;
//!
//! // Standard precision grid
//! let domain = IntervalClosed::new(
//!     RealNative64StrictFinite::try_from_f64(0.0).unwrap(),
//!     RealNative64StrictFinite::try_from_f64(1.0).unwrap()
//! );
//! let standard_grid = Grid1D::uniform(domain, NumIntervals::try_new(100).unwrap());
//!
//! # #[cfg(feature = "rug")]
//! # {
//! use num_valid::RealRugStrictFinite;
//!
//! // High precision grid (256-bit precision) - same API!
//! let hp_domain = IntervalClosed::new(
//!     RealRugStrictFinite::<256>::try_from_f64(0.0).unwrap(),
//!     RealRugStrictFinite::<256>::try_from_f64(1.0).unwrap()
//! );
//! let high_precision_grid = Grid1D::uniform(hp_domain, NumIntervals::try_new(100).unwrap());
//!
//! // Both grids support identical operations
//! let point = RealRugStrictFinite::<256>::try_from_f64(0.5).unwrap();
//! let interval_id = high_precision_grid.find_interval_id_of_point(&point);
//! # }
//! ```
//!
//! ### Backend Selection Guide
//!
//! | Backend | Scalar Type | Precision | Performance | Use Case |
//! |---------|-------------|-----------|-------------|----------|
//! | **Native** | [`f64`] | 53 bits | ⚡⚡⚡ Maximum | Production simulations, real-time |
//! | **Native Validated (in Debug only)** | [`num_valid::RealNative64StrictFiniteInDebug`] | 53 bits | ⚡⚡⚡ Maximum (same as `f64`) | Most applications |
//! | **Native Validated** | [`num_valid::RealNative64StrictFinite`] | 53 bits | ⚡⚡ Small penalty | Safety-critical applications |
//! | **Arbitrary Precision** | `num_valid::RealRugStrictFinite<N>` | N bits | ⚡ Configurable | Scientific research, symbolic math |
//!
//! **Enabling arbitrary precision**: Add `features = ["rug"]` to your `Cargo.toml`:
//! ```toml
//! [dependencies]
//! grid1d = { version = "0.3.3", features = ["rug"] }
//! ```
//!
//! ### Performance Through Specialization
//! - **Uniform grids**: O(1) point location with analytical formulas ([`Grid1DUniform`])
//! - **Non-uniform grids**: O(log n) binary search with optimized data structures ([`Grid1DNonUniform`])
//! - **Memory efficiency**: Zero-cost abstractions and cache-friendly memory layouts
//! - **SIMD-ready**: Regular patterns enable vectorization optimizations
//!
//! #### Performance Comparison: Grid Operations
//!
//! | Operation | [`Grid1DUniform`] | [`Grid1DNonUniform`] | [`Grid1DUnion`] | Notes |
//! |-----------|------------------|---------------------|----------------|--------|
//! | **Construction** | O(n) | O(n) | O(n+m) | Union merges coordinates from two grids |
//! | **Point location** | **O(1)** ⚡ | O(log n) | O(log(n+m)) | Uniform uses direct calculation |
//! | **Interval access** | O(1) | O(1) | O(1) | All grid types support direct access by ID |
//! | **Batch point location** | **O(k)** ⚡ | O(k log n) | O(k log(n+m)) | k = number of points |
//! | **Coordinate access** | O(1) | O(1) | O(1) | Direct array indexing |
//! | **Domain check** | O(1) | O(1) | O(1) | Simple boundary comparison |
//! | **Memory footprint** | 8n + 40 bytes | 8n + 40 bytes | 8(n+m) + 16k + 80 bytes | k = mapping table size |
//! | **Cache efficiency** | ⚡⚡⚡ Excellent | ⚡⚡ Good | ⚡⚡ Good | Sequential memory layout |
//!
//! *Where n = number of grid points, m = points in second grid (for union), k = number of intervals*
//!
//! #### Scalar Type Performance Impact
//!
//! | Scalar Type | Overhead | Point Location (uniform) | Memory per Point | When to Use |
//! |-------------|----------|--------------------------|------------------|-------------|
//! | `f64` | **0%** baseline | ~1.2 ns | 8 bytes | Maximum performance, trusted inputs |
//! | [`num_valid::RealNative64StrictFiniteInDebug`] | **0%** in release | ~1.2 ns | 8 bytes | **Recommended default** |
//! | [`num_valid::RealNative64StrictFinite`] | ~5-10% | ~1.3 ns | 8 bytes | Safety-critical applications |
//! | `RealRugStrictFinite<256>` | ~100-500% | ~5-50 ns | 32+ bytes | Arbitrary precision needs |
//!
//! *Benchmarked on x86_64 with typical workloads. Actual performance varies by CPU and workload.*
//!
//! ### Composability and Integration
//! - **Trait-based design**: Core traits ([`IntervalPartition`], [`HasDomain1D`], [`HasCoords1D`]) enable generic programming
//! - **Unified interface**: [`Grid1D`] enum provides seamless switching between uniform and non-uniform grids
//! - **Grid operations**: Union, refinement, and intersection operations with bidirectional mappings
//!
//! ## Quick Start Guide
//!
//! ### Creating Grids
//! ```rust
//! use grid1d::{
//!     Grid1D, HasCoords1D, IntervalPartition,
//!     intervals::*,
//!     scalars::{NumIntervals, IntervalId},
//! };
//! use sorted_vec::partial::SortedSet;
//! use std::ops::Deref;
//! use try_create::TryNew;
//!
//! // Create uniform grid with equal spacing
//! let domain = IntervalClosed::new(0.0, 1.0);
//! let uniform_grid = Grid1D::uniform(domain.clone(), NumIntervals::try_new(4).unwrap());
//! assert_eq!(uniform_grid.coords().deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
//!
//! // Create non-uniform grid with custom spacing
//! let coords = SortedSet::from_unsorted(vec![0.0, 0.3, 0.7, 0.9, 1.0]);
//! let non_uniform_grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
//!
//! // Both grids implement the same IntervalPartition trait
//! assert_eq!(uniform_grid.num_intervals().as_ref(), &4);
//! assert_eq!(non_uniform_grid.num_intervals().as_ref(), &4);
//! ```
//!
//! ### Point Location and Interval Access
//! ```rust
//! use grid1d::{Grid1D, IntervalPartition, intervals::*, scalars::NumIntervals};
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(IntervalClosed::new(0.0, 10.0), NumIntervals::try_new(10).unwrap());
//!
//! // Efficient point location
//! let interval_id = grid.find_interval_id_of_point(&3.7);
//! println!("Point 3.7 is in interval {}", interval_id.as_ref());
//!
//! // Access interval properties
//! let interval = grid.interval(&interval_id);
//! let length = grid.interval_length(&interval_id);
//! println!("Interval length: {}", length.as_ref());
//!
//! // Batch point location
//! let points = vec![1.5, 4.2, 7.8, 9.1];
//! let intervals = grid.find_intervals_for_points(&points);
//! ```
//!
//! ### Working with Different Scalar Types
//! ```rust
//! use grid1d::{Grid1D, intervals::*, scalars::NumIntervals};
//! use num_valid::RealNative64StrictFiniteInDebug;
//! use try_create::TryNew;
//!
//! // Use validated scalar types for debug safety
//! type Real = RealNative64StrictFiniteInDebug;
//! let domain = IntervalClosed::new(
//!     Real::try_new(0.0).unwrap(),
//!     Real::try_new(1.0).unwrap()
//! );
//! let grid = Grid1D::uniform(domain, NumIntervals::try_new(100).unwrap());
//! ```
//!
//! ## Core Concepts
//!
//! ### Grid Types and Performance Characteristics
//!
//! | Grid Type | Point Location | Memory | Best For |
//! |-----------|----------------|---------|----------|
//! | **[`Grid1DUniform`]** | O(1) analytical | O(n) | Equal spacing, maximum performance |
//! | **[`Grid1DNonUniform`]** | O(log n) binary search | O(n) | Adaptive spacing, complex features |
//! | **[`Grid1DUnion`]** | O(log n) on unified grid | O(n+m) | Multi-physics, grid combination |
//!
//! ### Interval Partition Types
//! The library correctly handles different interval boundary semantics:
//! - **Closed intervals `[a,b]`**: Include both endpoints
//! - **Open intervals `(a,b)`**: Exclude both endpoints  
//! - **Semi-open intervals**: `[a,b)` or `(a,b]` - include one endpoint
//! - **Automatic sub-interval construction**: Preserves domain semantics in partitions
//!
//! ### Mathematical Properties and Guarantees
//! Every grid maintains these invariants:
//! - **Complete domain coverage**: ⋃ intervals = domain (exact reconstruction)
//! - **Non-overlapping intervals**: Intervals are disjoint except at boundaries
//! - **Unique point location**: Every domain point belongs to exactly one interval
//! - **Boundary consistency**: First/last coordinates match domain bounds exactly
//! - **Sorted coordinates**: Points are in strict ascending order
//!
//! ## Advanced Features
//!
//! ### Adaptive Mesh Refinement
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//! use std::collections::BTreeMap;
//!
//! let base_grid = Grid1D::uniform(IntervalClosed::new(0.0, 1.0), NumIntervals::try_new(4).unwrap());
//!
//! // Uniform refinement: double resolution everywhere
//! let uniform_refinement = base_grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
//! assert_eq!(uniform_refinement.refined_grid().num_intervals().as_ref(), &8);
//!
//! // Selective refinement: refine only specific intervals
//! let selective_plan = BTreeMap::from([
//!     (IntervalId::new(1), PositiveNumPoints1D::try_new(2).unwrap()), // 3 sub-intervals
//!     (IntervalId::new(3), PositiveNumPoints1D::try_new(1).unwrap()), // 2 sub-intervals
//! ]);
//! let selective_refinement = uniform_refinement.into_refined_grid().refine(&selective_plan);
//! ```
//!
//! ### Grid Union for Multi-Physics
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//! use try_create::TryNew;
//!
//! let domain = IntervalClosed::new(0.0, 2.0);
//!
//! // Physics A: coarse global grid
//! let grid_a = Grid1D::uniform(domain, NumIntervals::try_new(4).unwrap());
//!
//! // Physics B: fine local grid
//! let coords_b = SortedSet::from_unsorted(vec![0.0, 0.3, 0.7, 1.1, 1.4, 2.0]);
//! let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords_b).unwrap();
//!
//! // Create unified grid preserving mappings to both original grids
//! let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
//! println!("Unified grid has {} intervals", union.num_refined_intervals().as_ref());
//!
//! // Map data between grids
//! for (refined_id, a_id, b_id) in union.iter_interval_mappings() {
//!     println!("Unified interval {} maps to A[{}] and B[{}]",
//!              refined_id.as_ref(), a_id.as_ref(), b_id.as_ref());
//! }
//! ```
//!
//! ### Domain-Specific Interval Semantics
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//!
//! // Closed domain: [0,2] includes both boundaries
//! let closed_grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
//!     SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])
//! ).unwrap();
//!
//! // Open domain: (0,2) excludes both boundaries  
//! let open_grid = Grid1D::<IntervalOpen<f64>>::try_from_sorted(
//!     SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])
//! ).unwrap();
//!
//! // Different sub-interval types are automatically constructed
//! let closed_interval_0 = closed_grid.interval(&IntervalId::new(0)); // [0,1]
//! let open_interval_0 = open_grid.interval(&IntervalId::new(0));     // (0,1]
//!
//! // Point containment reflects domain semantics
//! use grid1d::intervals::IntervalTrait;
//! // Note: This is a simplified example - actual containment checking requires pattern matching
//! ```
//!
//! ### High-Performance Computations
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! // Large uniform grid for finite difference methods
//! let grid = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 1.0),
//!     NumIntervals::try_new(100_000).unwrap()
//! );
//!
//! // O(1) point location for uniform grids
//! let start = std::time::Instant::now();
//! let test_points: Vec<f64> = (0..10_000).map(|i| i as f64 / 10_000.0).collect();
//! let intervals: Vec<IntervalId> = test_points.iter()
//!     .map(|&p| grid.find_interval_id_of_point(&p))
//!     .collect();
//! let elapsed = start.elapsed();
//! println!("Located {} points in {:?}", test_points.len(), elapsed);
//! ```
//!
//! ## Performance Characteristics
//!
//! ### Time Complexity Summary
//! | Operation | Uniform Grid | Non-Uniform Grid | Grid Union | Notes |
//! |-----------|-------------|------------------|------------|--------|
//! | **Grid Creation** | O(n) | O(n) | O(n+m) | Linear in point count |
//! | **Point Location** | **O(1)** | O(log n) | O(log n) | Uniform grids use analytical formulas |
//! | **Interval Access** | O(1) | O(1) | O(1) | Direct array indexing |
//! | **Intersection** | O(k) | O(log n + k) | O(log n + k) | Where k = result size |
//! | **Union Creation** | - | - | O(n+m) | Optimal merge algorithm |
//! | **Grid Refinement** | O(n×r) | O(n×r) | O(n×r) | Where r = refinement factor |
//!
//! ### Space Complexity
//! - **Coordinates storage**: O(n) for all grid types
//! - **Grid union mappings**: O(n+m) additional space for bidirectional mappings
//! - **Refinement structures**: O(refined intervals) for mapping preservation
//! - **Memory layout**: Cache-friendly contiguous arrays for optimal performance
//!
//! ## Integration with other numerical Ecosystem
//!
//! ### Scalar Type Recommendations
//! Choose scalar types based on your performance and safety requirements:
//!
//! ```rust
//! use grid1d::{Grid1D, intervals::*, scalars::NumIntervals};
//! use num_valid::{RealNative64StrictFiniteInDebug, RealNative64StrictFinite};
//! use try_create::TryNew;
//!
//! // Maximum performance (production)
//! let fast_grid = Grid1D::uniform(
//!     IntervalClosed::new(0.0_f64, 1.0_f64),
//!     NumIntervals::try_new(1000).unwrap()
//! );
//!
//! // Balanced performance + debug safety (recommended)
//! type Real = RealNative64StrictFiniteInDebug;
//! let safe_grid = Grid1D::uniform(
//!     IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap()),
//!     NumIntervals::try_new(1000).unwrap()
//! );
//!
//! // Always validated (safety-critical)
//! type SafeReal = RealNative64StrictFinite;
//! let secure_grid = Grid1D::uniform(
//!     IntervalClosed::new(SafeReal::try_new(0.0).unwrap(), SafeReal::try_new(1.0).unwrap()),
//!     NumIntervals::try_new(1000).unwrap()
//! );
//! ```
//!
//! ### Numerical Method Integration
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! // Finite difference method setup
//! fn setup_finite_difference(domain: IntervalClosed<f64>, resolution: usize)
//!     -> (Grid1D<IntervalClosed<f64>>, Vec<f64>)
//! {
//!     let grid = Grid1D::uniform(domain, NumIntervals::try_new(resolution).unwrap());
//!     let solution = vec![0.0; grid.coords().len()];
//!     (grid, solution)
//! }
//!
//! // Spectral method with periodic boundaries
//! fn setup_spectral_method(n_modes: usize) -> Grid1D<IntervalLowerClosedUpperOpen<f64>> {
//!     let domain = IntervalLowerClosedUpperOpen::new(0.0, 2.0 * std::f64::consts::PI);
//!     Grid1D::uniform(domain, NumIntervals::try_new(n_modes).unwrap())
//! }
//! ```
//!
//! ## Error Handling and Validation
//!
//! ### Construction Error Types
//! The library provides detailed error information for debugging:
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//!
//! // Insufficient points
//! let coords = SortedSet::from_unsorted(vec![0.0]); // Only one point
//! match Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords) {
//!     Err(ErrorsGrid1D::RequiredAtLeastTwoDistinctPoints { num_points_provided, .. }) => {
//!         println!("Need at least 2 points, got {}", num_points_provided);
//!     }
//!     _ => unreachable!(),
//! }
//! ```
//!
//! ### Safe Construction Patterns
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//!
//! fn safe_grid_from_points(
//!     points: Vec<f64>,
//! ) -> Result<Grid1D<IntervalClosed<f64>>, String> {
//!     if points.len() < 2 {
//!         return Err("Need at least 2 points".to_string());
//!     }
//!     
//!     let coords = SortedSet::from_unsorted(points);
//!     Grid1D::try_from_sorted(coords)
//!         .map_err(|e| format!("Grid construction failed: {}", e))
//! }
//! ```
//!
//! ## Common Pitfalls and How to Avoid Them
//!
//! ### Pitfall 1: Assuming All Grids Are Uniform
//! **Problem**: Trying to access uniform-specific methods on non-uniform grids.
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//! use try_create::TryNew;
//!
//! let coords = SortedSet::from_unsorted(vec![0.0, 0.1, 0.5, 1.0]);
//! let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
//!
//! // ERROR: grid is NonUniform, not Uniform!
//! // let delta = grid.delta_points(); // This won't compile
//! ```
//!
//! **Solution**: Use pattern matching or trait methods:
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//! use try_create::TryNew;
//!
//! let coords = SortedSet::from_unsorted(vec![0.0, 0.1, 0.5, 1.0]);
//! let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
//!
//! // Correct: Use pattern matching
//! match &grid {
//!     Grid1D::Uniform(uniform_grid) => {
//!         let delta = uniform_grid.delta_points();
//!         println!("Uniform spacing: {}", delta);
//!     }
//!     Grid1D::NonUniform(_) => {
//!         println!("Non-uniform grid - use general methods");
//!     }
//! }
//!
//! // Or use trait methods that work for all grid types
//! let interval_id = grid.find_interval_id_of_point(&0.3); // Works for both uniform and non-uniform
//! ```
//!
//! ### Pitfall 2: Boundary Point Confusion with Open Intervals
//! **Problem**: Assuming boundary points are always contained in the domain.
//!
//! ```rust,should_panic
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! // Open interval (0, 1) excludes both 0.0 and 1.0
//! let grid = Grid1D::uniform(
//!     IntervalOpen::new(0.0, 1.0),
//!     NumIntervals::try_new(4).unwrap()
//! );
//!
//! // This panics in debug mode - 0.0 is NOT in the open interval!
//! let id = grid.find_interval_id_of_point(&0.0);
//! ```
//!
//! **Solution**: Use `try_` variants or check domain containment:
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(
//!     IntervalOpen::new(0.0, 1.0),
//!     NumIntervals::try_new(4).unwrap()
//! );
//!
//! // Correct: Use try_ variant
//! if let Some(id) = grid.try_find_interval_id_of_point(&0.0) {
//!     println!("Point is in domain at interval {}", id.as_ref());
//! } else {
//!     println!("Point is outside the domain");
//! }
//!
//! // Or check containment first
//! let point = 0.0;
//! if grid.domain().contains_point(&point) {
//!     let id = grid.find_interval_id_of_point(&point);
//!     // Process...
//! }
//! ```
//!
//! ### Pitfall 3: Mixing Different Scalar Types
//! **Problem**: Trying to compare or combine grids with different scalar types.
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use num_valid::RealNative64StrictFinite;
//! use try_create::TryNew;
//!
//! let grid_f64 = Grid1D::uniform(
//!     IntervalClosed::new(0.0_f64, 1.0_f64),
//!     NumIntervals::try_new(10).unwrap()
//! );
//!
//! let grid_validated = Grid1D::uniform(
//!     IntervalClosed::new(
//!         RealNative64StrictFinite::try_new(0.0).unwrap(),
//!         RealNative64StrictFinite::try_new(1.0).unwrap()
//!     ),
//!     NumIntervals::try_new(10).unwrap()
//! );
//!
//! // ERROR: Cannot create union of grids with different scalar types
//! // let union = Grid1DUnion::try_new(&grid_f64, &grid_validated);
//! ```
//!
//! **Solution**: Use consistent scalar types throughout your code:
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use num_valid::RealNative64StrictFiniteInDebug;
//! use try_create::TryNew;
//!
//! // Choose one scalar type and use it consistently
//! type Real = RealNative64StrictFiniteInDebug;
//!
//! let grid1 = Grid1D::uniform(
//!     IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap()),
//!     NumIntervals::try_new(10).unwrap()
//! );
//!
//! let grid2 = Grid1D::uniform(
//!     IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(2.0).unwrap()),
//!     NumIntervals::try_new(20).unwrap()
//! );
//!
//! // Now operations between grids work correctly
//! // let union = Grid1DUnion::try_new(&grid1, &grid2); // Would work if domains matched
//! ```
//!
//! ### Pitfall 4: Inefficient Point Location in Loops
//! **Problem**: Repeatedly calling point location for sequential points.
//!
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 1.0),
//!     NumIntervals::try_new(1000).unwrap()
//! );
//!
//! let points: Vec<f64> = (0..1000).map(|i| i as f64 / 1000.0).collect();
//!
//! // Inefficient: Individual calls in a loop
//! for point in &points {
//!     let _ = grid.find_interval_id_of_point(point);
//! }
//! ```
//!
//! **Solution**: Use batch operations when available:
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! let grid = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 1.0),
//!     NumIntervals::try_new(1000).unwrap()
//! );
//!
//! let points: Vec<f64> = (0..1000).map(|i| i as f64 / 1000.0).collect();
//!
//! // Efficient: Batch operation
//! let interval_ids = grid.find_intervals_for_points(&points);
//! ```
//!
//! ### Pitfall 5: Pattern Matching on Interval Enums
//! **Problem**: Working with heterogeneous interval collections or runtime-determined interval types.
//!
//! The [`Interval`] enum wraps all interval types for uniform handling,
//! but requires pattern matching to access type-specific functionality.
//!
//! ```rust
//! use grid1d::intervals::*;
//!
//! // Creating intervals of different types
//! let closed = IntervalClosed::new(0.0, 1.0);
//! let open = IntervalOpen::new(0.0, 1.0);
//! let half_open = IntervalLowerClosedUpperOpen::new(0.0, 1.0);
//!
//! // Convert to universal Interval enum for uniform handling
//! let intervals: Vec<Interval<f64>> = vec![
//!     closed.into(),
//!     open.into(),
//!     half_open.into(),
//! ];
//!
//! // Pattern matching on Interval enum
//! for interval in &intervals {
//!     match interval {
//!         Interval::FiniteLength(finite) => {
//!             // Further pattern matching on finite length variants
//!             match finite {
//!                 IntervalFiniteLength::PositiveLength(pos_len) => {
//!                     match pos_len {
//!                         IntervalFinitePositiveLength::Closed(closed) => {
//!                             println!("[{}, {}]", closed.lower_bound_value(), closed.upper_bound_value());
//!                         }
//!                         IntervalFinitePositiveLength::Open(open) => {
//!                             println!("({}, {})", open.lower_bound_value(), open.upper_bound_value());
//!                         }
//!                         IntervalFinitePositiveLength::LowerClosedUpperOpen(hco) => {
//!                             println!("[{}, {})", hco.lower_bound_value(), hco.upper_bound_value());
//!                         }
//!                         IntervalFinitePositiveLength::LowerOpenUpperClosed(hoc) => {
//!                             println!("({}, {}]", hoc.lower_bound_value(), hoc.upper_bound_value());
//!                         }
//!                     }
//!                 }
//!                 IntervalFiniteLength::ZeroLength(singleton) => {
//!                     println!("{{{}}}", singleton.value());
//!                 }
//!             }
//!         }
//!         Interval::InfiniteLength(infinite) => {
//!             println!("Unbounded interval: {:?}", infinite);
//!         }
//!     }
//! }
//! ```
//!
//! **Solution for common operations**: Use trait methods that work across all interval types:
//! ```rust
//! use grid1d::intervals::*;
//!
//! let intervals: Vec<Interval<f64>> = vec![
//!     IntervalClosed::new(0.0, 1.0).into(),
//!     IntervalOpen::new(2.0, 3.0).into(),
//!     IntervalLowerClosedUpperOpen::new(4.0, 5.0).into(),
//! ];
//!
//! // Use trait methods for common operations without pattern matching
//! for interval in &intervals {
//!     // contains_point works on all interval types via IntervalTrait
//!     let contains_half = interval.contains_point(&0.5);
//!     println!("Contains 0.5: {}", contains_half);
//!
//!     // Get bounds info using GetLeftValue/GetRightValue traits
//!     // (only for bounded intervals)
//!     if let Interval::FiniteLength(finite) = interval {
//!         match finite {
//!             IntervalFiniteLength::PositiveLength(pos) => {
//!                 let length = pos.length();
//!                 println!("Length: {}", length.as_ref());
//!             }
//!             IntervalFiniteLength::ZeroLength(_) => {
//!                 println!("Length: 0 (singleton)");
//!             }
//!         }
//!     }
//! }
//! ```
//!
//! **Simplified approach for finite intervals only**:
//! ```rust
//! use grid1d::intervals::*;
//!
//! // If you only work with finite positive-length intervals, use IntervalFinitePositiveLength
//! let finite_intervals: Vec<IntervalFinitePositiveLength<f64>> = vec![
//!     IntervalClosed::new(0.0, 1.0).into(),
//!     IntervalOpen::new(2.0, 3.0).into(),
//! ];
//!
//! for interval in &finite_intervals {
//!     // length() is available directly on IntervalFinitePositiveLength
//!     let len = interval.length();
//!     let mid = interval.midpoint();
//!     println!("Length: {}, Midpoint: {}", len.as_ref(), mid);
//! }
//! ```
//!
//! ## Best Practices and Recommendations
//!
//! ### When to Use Each Grid Type
//!
//! **Use [`Grid1DUniform`] when:**
//! - Finite difference methods requiring equal spacing
//! - Spectral methods with uniform point distributions
//! - Maximum performance is critical (O(1) point location)
//! - Simple problems with uniform solution characteristics
//!
//! **Use [`Grid1DNonUniform`] when:**
//! - Adaptive mesh refinement based on solution features
//! - Boundary layer resolution in fluid dynamics
//! - Multi-scale problems requiring variable resolution
//! - Complex geometries or irregular point distributions
//!
//! **Use [`Grid1DUnion`] when:**
//! - Combining grids from different physics solvers
//! - Multi-physics coupling requiring consistent discretization
//! - Grid overlay operations for interpolation
//! - Preserving relationships between multiple grid levels
//!
//! ### Performance Optimization Tips
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! // 1. Choose uniform grids when possible for O(1) point location
//! let uniform_grid = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 1.0),
//!     NumIntervals::try_new(1000).unwrap()
//! );
//!
//! // 2. Use batch operations for multiple points
//! let points = vec![0.1, 0.3, 0.7, 0.9];
//! let intervals = uniform_grid.find_intervals_for_points(&points); // More efficient than individual calls
//!
//! // 3. Pre-allocate vectors for performance-critical loops
//! let mut results = Vec::with_capacity(points.len());
//! for &point in &points {
//!     results.push(uniform_grid.find_interval_id_of_point(&point));
//! }
//!
//! // 4. Use direct coordinate access for sequential processing
//! let coords = uniform_grid.coords().as_ref();
//! for i in 0..coords.len()-1 {
//!     let interval_length = coords[i+1] - coords[i];
//!     // Process interval...
//! }
//! ```
//!
//! ### Memory Management
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! // For large grids, consider the memory footprint
//! let large_grid = Grid1D::uniform(
//!     IntervalClosed::new(0.0, 1.0),
//!     NumIntervals::try_new(1_000_000).unwrap()
//! );
//!
//! // Memory usage: ~8MB for coordinates (1M points × 8 bytes/f64)
//! println!("Grid memory: ~{} MB",
//!          large_grid.coords().len() * std::mem::size_of::<f64>() / 1_000_000);
//!
//! // Use iterators to avoid creating intermediate vectors
//! let sum: f64 = large_grid.coords().iter().sum();
//! let mean = sum / large_grid.coords().len() as f64;
//! ```
//!
//! ## Related Modules and Integration
//!
//! This module integrates with the other numerical ecosystem:
//!
//! - **[`crate::intervals`]**: Provides interval types with boundary semantics
//! - **[`crate::scalars`]**: Type-safe wrappers for numerical quantities
//! - **[`crate::bounds`]**: Low-level boundary representation
//! - **[`crate::topology`]**: Topology types, index spaces, and adjacency queries
//! - **[`num_valid`]**: Generic scalar arithmetic and validation
//! - **[`sorted_vec`]**: Efficient sorted container implementation
//!
//! ### Future Extensions
//! The design supports natural extensions to:
//! - **Multi-dimensional grids**: 2D/3D generalizations using the same principles
//! - **Adaptive algorithms**: Error-driven refinement strategies
//! - **Parallel processing**: Thread-safe grid operations and domain decomposition
//! - **GPU acceleration**: SIMD-friendly data layouts for vectorization
//! - **File I/O**: Serialization for large-scale scientific computing
//!
//! ## Examples and Use Cases
//!
//! ### Computational Fluid Dynamics
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use sorted_vec::partial::SortedSet;
//!
//! // Create boundary layer mesh with fine spacing near walls
//! fn create_boundary_layer_mesh(domain_length: f64, boundary_thickness: f64)
//!     -> Grid1D<IntervalClosed<f64>>
//! {
//!     let mut points = vec![0.0];
//!     
//!     // Fine spacing in boundary layer
//!     for i in 1..=20 {
//!         let eta = (i as f64) / 20.0;
//!         let y = boundary_thickness * (eta * eta); // Quadratic clustering
//!         points.push(y);
//!     }
//!     
//!     // Coarse spacing in outer region
//!     for i in 1..=10 {
//!         let y = boundary_thickness + (domain_length - boundary_thickness) * (i as f64) / 10.0;
//!         points.push(y);
//!     }
//!     
//!     Grid1D::<IntervalClosed<f64>>::try_from_sorted(
//!         SortedSet::from_unsorted(points)
//!     ).unwrap()
//! }
//!
//! let cfd_grid = create_boundary_layer_mesh(1.0, 0.1);
//! println!("CFD grid: {} intervals, min spacing: {:.2e}",
//!          cfd_grid.num_intervals().as_ref(),
//!          cfd_grid.min_interval_length().as_ref());
//! ```
//!
//! ### Adaptive Finite Element Analysis
//! ```rust
//! use grid1d::{*, intervals::*, scalars::*};
//! use try_create::TryNew;
//!
//! // Simulate adaptive refinement based on error indicators
//! fn adaptive_fem_cycle(
//!     mut grid: Grid1D<IntervalClosed<f64>>,
//!     error_threshold: f64,
//!     max_iterations: usize
//! ) -> Grid1D<IntervalClosed<f64>> {
//!     for iteration in 0..max_iterations {
//!         // Compute error indicators (simplified)
//!         let errors: Vec<f64> = (0..*grid.num_intervals().as_ref())
//!             .map(|i| {
//!                 // Simulate error based on interval size (larger intervals = higher error)
//!                 grid.interval_length(&IntervalId::new(i)).as_ref() * 10.0
//!             })
//!             .collect();
//!         
//!         // Mark intervals for refinement
//!         let refinement_plan: std::collections::BTreeMap<IntervalId, PositiveNumPoints1D> = errors
//!             .iter()
//!             .enumerate()
//!             .filter_map(|(i, &error)| {
//!                 if error > error_threshold {
//!                     Some((IntervalId::new(i), PositiveNumPoints1D::try_new(1).unwrap()))
//!                 } else {
//!                     None
//!                 }
//!             })
//!             .collect();
//!         
//!         if refinement_plan.is_empty() {
//!             println!("Converged after {} iterations", iteration);
//!             break;
//!         }
//!         
//!         // Apply refinement
//!         let refinement = grid.refine(&refinement_plan);
//!         grid = Grid1D::NonUniform(refinement.into_refined_grid());
//!         
//!         println!("Iteration {}: refined {} intervals, total: {}",
//!                  iteration, refinement_plan.len(), grid.num_intervals().as_ref());
//!     }
//!     
//!     grid
//! }
//! ```
//!
//! For comprehensive examples, API documentation, and advanced usage patterns,
//! see the individual struct and trait documentation throughout this module.

pub mod bounds;
pub mod coords;
pub use coords::{
    Coords1D, Coords1DInDomain, Coords1DInDomainArc, Coords1DInDomainBorrowed, Coords1DInDomainBox,
    Coords1DInDomainCow, Coords1DInDomainOwned, Coords1DInDomainRc, ErrorsCoords1D,
    ErrorsCoords1DInDomain, linspace, logspace,
};

pub mod intervals;
pub mod operations;
pub use operations::{
    refinement::{Grid1DNonUniformRefinement, Grid1DRefinement, Grid1DUniformRefinement},
    union::{ErrorsGrid1DUnion, Grid1DUnion},
};

pub mod scalars;

pub mod topology;
pub use topology::{
    Adjacency1D, Grid1DIndexSpaces, HasSide, IndexSpace1D, IntervalIndexSpace1D, IntervalSide,
    PointIndexSpace1D, Side, Topology1D,
};

pub mod traits;
pub use traits::{
    BuildIntervalInPartition, HasCoords1D, HasDomain1D, IntervalPartition, TransformedPoints1D,
    TransformedPoints1DArc, TransformedPoints1DBorrowed, TransformedPoints1DOwned,
    TransformedPoints1DRc,
};

use crate::{
    intervals::{
        Contains, ErrorsIntervalConstruction, GetLowerBoundValue, GetUpperBoundValue, Interval,
        IntervalTrait,
        bounded::{IntervalFinitePositiveLengthTrait, SubIntervalInPartition},
    },
    operations::refinement::{
        build_non_uniform_grid_refinement, refine_all_intervals_in_coords1d_same_subdivisions,
    },
    scalars::{IntervalId, NumIntervals, PositiveNumPoints1D},
    topology::sealed::SupportsCircularTopology,
};
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;

//------------------------------------------------------------------------------------------------------------------------------
///  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 [`IntervalPartition`] 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, IntervalPartition,
///     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::{
///     Grid1DUniform, IntervalPartition, HasCoords1D,
///     intervals::*,
///     scalars::{IntervalId, NumIntervals},
/// };
/// 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.);
/// assert_eq!(interval_id, IntervalId::new(0));
///
/// let interval_id = grid.find_interval_id_of_point(&0.0009);
/// assert_eq!(interval_id, IntervalId::new(0));
///
/// let interval_id = grid.find_interval_id_of_point(&0.001);
/// 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);
/// assert_eq!(interval_id, IntervalId::new(1));
///
/// let interval_id = grid.find_interval_id_of_point(&0.002);
/// 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);
/// 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);
/// 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);
/// assert_eq!(interval_id, IntervalId::new(9999));
///
/// let interval_id = grid.find_interval_id_of_point(&10.000);
/// assert_eq!(interval_id, IntervalId::new(9999));
///
/// ```
///
/// ### Different Domain Types
/// ```rust
/// use grid1d::{
///     Grid1DUniform, IntervalPartition,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// 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, IntervalPartition,
///     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::{
///     Grid1DUniform, IntervalPartition,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// 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))
///     .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, IntervalPartition, 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::{
///     Grid1D, Grid1DUniform, IntervalPartition,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// 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 IntervalPartition 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);
/// 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::{
///     Grid1DUniform, IntervalPartition,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// 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::{
///     Grid1DUniform, IntervalPartition,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// 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::{
///     Grid1DUniform, IntervalPartition, HasDomain1D,
///     intervals::*,
///     scalars::{NumIntervals, IntervalId},
/// };
/// 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::{
///     IntervalPartition, 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, IntervalPartition,
///     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
/// - [`IntervalPartition`]: 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`].
    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(num_intervals, Topology1D::RealLine);

        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`].
    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 = Grid1DIndexSpaces::new(num_intervals, Topology1D::Circle);

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

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

    fn domain(&self) -> &Domain1D {
        &self.domain
    }
}

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

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

impl<Domain1D: BuildIntervalInPartition> IntervalPartition for Grid1DUniform<Domain1D> {
    fn index_spaces(&self) -> &Grid1DIndexSpaces {
        &self.index_spaces
    }

    /// 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
    /// [`IntervalPartition`] 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::{
    ///     Grid1DUniform, IntervalPartition,
    ///     intervals::*,
    ///     scalars::{NumIntervals, IntervalId},
    /// };
    /// 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.try_find_interval_id_of_point(&0.0).unwrap(), IntervalId::new(0)); // Domain lower bound
    /// assert_eq!(grid.try_find_interval_id_of_point(&0.5).unwrap(), IntervalId::new(0)); // Interior point
    /// assert_eq!(grid.try_find_interval_id_of_point(&1.0).unwrap(), IntervalId::new(1)); // p₁ → right interval (left-closed)
    /// assert_eq!(grid.try_find_interval_id_of_point(&1.5).unwrap(), IntervalId::new(1)); // Interior point
    /// assert_eq!(grid.try_find_interval_id_of_point(&2.0).unwrap(), IntervalId::new(1)); // Domain upper bound
    /// ```
    fn try_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(grid_idx) => {
                let grid_idx = grid_idx.into_inner();

                // Exact match with a grid point
                if grid_idx == 0 {
                    // This corresponds to the domain lower bound
                    Some(IntervalId::new(0))
                } else if grid_idx == num_intervals {
                    // This corresponds to the domain upper bound
                    Some(IntervalId::new(num_intervals - 1))
                } else {
                    // Interior partition point: for left-closed sub-intervals [pₖ, pₖ₊₁)
                    // it belongs to the right interval (index = grid_idx);
                    // for right-closed sub-intervals (pₖ, pₖ₊₁] it belongs to the left (grid_idx - 1).
                    if Domain1D::LOWER_BOUND_OWNS_BOUNDARY {
                        Some(IntervalId::new(grid_idx))
                    } else {
                        Some(IntervalId::new(grid_idx - 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))
            }
        }
    }

    /// 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
    /// ```
    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
    /// ```
    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)
    }
}
//-------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------------------------------------
///  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 `IntervalPartition` 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 `BuildIntervalInPartition` trait.
///
/// ## Construction and Usage Examples
///
/// ### Basic Construction
/// ```rust
/// use grid1d::{
///     Grid1DNonUniform, HasCoords1D, IntervalPartition,
///     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(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, IntervalPartition, 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(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(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.
/// - [`IntervalPartition`]: 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.
    domain: Domain1D,

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

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

impl<Domain1D: IntervalFinitePositiveLengthTrait> Grid1DNonUniform<Domain1D> {
    fn create_domain_from_coords(
        coords: &Coords1D<Domain1D::RealType>,
    ) -> Result<Domain1D, ErrorsGrid1D<Domain1D>> {
        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(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(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(num_intervals, Topology1D::RealLine);

        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(
        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(num_intervals, Topology1D::Circle);

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

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

    fn domain(&self) -> &Domain1D {
        &self.domain
    }
}

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

    fn coords(&self) -> &Coords1D<Self::Point1DType> {
        &self.coords
    }
}

impl<Domain1D> IntervalPartition for Grid1DNonUniform<Domain1D>
where
    Domain1D: BuildIntervalInPartition,
{
    type UniformlyRefinedGrid1DType = Self;

    fn index_spaces(&self) -> &Grid1DIndexSpaces {
        &self.index_spaces
    }

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

        if <Self::Domain1D as BuildIntervalInPartition>::LOWER_BOUND_OWNS_BOUNDARY {
            // 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 {
                return Some(IntervalId::new(0));
            }
            Some(IntervalId::new(idx - 1))
        }
    }

    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 refined_grid =
            Grid1DNonUniform::try_new(refined_coords).expect("Failed to create a refined grid!");

        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)
    }
}
//-------------------------------------------------------------------------------

///  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 { .. }));
/// ```
#[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>,
    },
}

/// 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 IntervalPartition**: The [`Grid1D`] enum implements the [`IntervalPartition`] trait, providing
///   methods for accessing intervals and performing operations on the grid.
///
/// # Examples
///
/// ## Creating a Uniform Grid
///
/// ```rust
/// use grid1d::{
///     Grid1D, HasCoords1D, IntervalPartition,
///     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 { .. }));
/// ```
#[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, IntervalPartition,
    ///     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.]);
    /// ```
    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>> {
        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>> {
        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 {
            // This unwrap is safe because we've already checked that n_pts >= 2.
            let num_intervals = NumIntervals::try_new(n_pts - 1).unwrap();
            let domain = Domain1D::try_new(
                coords_slice
                    .first()
                    .expect("Coordinates should not be empty")
                    .clone(),
                coords_slice
                    .last()
                    .expect("Coordinates should not be empty")
                    .clone(),
            )?;

            Ok(Grid1D::Uniform(Grid1DUniform::new(domain, num_intervals)))
        } else {
            Ok(Grid1D::NonUniform(Grid1DNonUniform::try_new(coords)?))
        }
    }
}

/// **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;

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

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

    fn coords(&self) -> &Coords1D<Domain1D::RealType> {
        match self {
            Grid1D::Uniform(pts) => pts.coords(),
            Grid1D::NonUniform(pts) => pts.coords(),
        }
    }
}

impl<Domain1D: BuildIntervalInPartition> IntervalPartition for Grid1D<Domain1D> {
    type UniformlyRefinedGrid1DType = Self;

    fn index_spaces(&self) -> &Grid1DIndexSpaces {
        match self {
            Grid1D::Uniform(grid) => grid.index_spaces(),
            Grid1D::NonUniform(grid) => grid.index_spaces(),
        }
    }

    /// 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 try_find_interval_id_of_point(&self, x: &Self::Point1DType) -> Option<IntervalId> {
        match self {
            Grid1D::Uniform(grid) => grid.try_find_interval_id_of_point(x),
            Grid1D::NonUniform(grid) => grid.try_find_interval_id_of_point(x),
        }
    }

    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(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::*, intervals::*, operations::union::Grid1DUnion};
    use std::ops::Deref;

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

    use super::*;

    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::Only(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());
                assert_eq!(i0, IntervalId::new(0));

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

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

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

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

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

                let i6 = grid1d.find_interval_id_of_point(&Real::try_new(3.0).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(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(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(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()),
                    IntervalId::new(0)
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(1.5).unwrap()),
                    IntervalId::new(1)
                );
                assert_eq!(
                    grid.find_interval_id_of_point(&Real::try_new(2.5).unwrap()),
                    IntervalId::new(2)
                );

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

            #[test]
            #[cfg(debug_assertions)]
            #[should_panic(expected = "The value RealValidated(0) is not in the interval")]
            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();

                // This should panic because 0.0 is not in the open interval (0, 2)
                grid.find_interval_id_of_point(&Real::try_new(0.0).unwrap());
            }
        }

        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(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(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.try_find_interval_id_of_point(&Real::try_new(0.0).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(10.0).unwrap()),
                    Some(IntervalId::new(9))
                );

                // Test interior points
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(0.5).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(1.5).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.try_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.try_find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(1))
                ); // p₁ → right interval
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(2.0).unwrap()),
                    Some(IntervalId::new(2))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(9.0).unwrap()),
                    Some(IntervalId::new(9))
                );

                // Test points outside domain
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(-0.1).unwrap()),
                    None
                );
                assert_eq!(
                    grid.try_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.try_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.try_find_interval_id_of_point(&Real::try_new(0.0001).unwrap()),
                    Some(IntervalId::new(1))
                );
                // Point just after first node
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(0.00011).unwrap()),
                    Some(IntervalId::new(1))
                );
                // Point in the middle
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(0.54321).unwrap()),
                    Some(IntervalId::new(5432))
                );
                // Point near end
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(0.99999).unwrap()),
                    Some(IntervalId::new(9999))
                );
                // Point on end
                assert_eq!(
                    grid.try_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.try_find_interval_id_of_point(&Real::try_new(0.1).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(1.9).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(0.0).unwrap()),
                    None
                ); // Outside open domain
                assert_eq!(
                    grid.try_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(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(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.try_find_interval_id_of_point(&Real::try_new(0.0).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(10.0).unwrap()),
                    Some(IntervalId::new(3))
                );

                // Interior points
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(0.5).unwrap()),
                    Some(IntervalId::new(0))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(2.5).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(5.0).unwrap()),
                    Some(IntervalId::new(2))
                );
                assert_eq!(
                    grid.try_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.try_find_interval_id_of_point(&Real::try_new(1.0).unwrap()),
                    Some(IntervalId::new(1))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(3.0).unwrap()),
                    Some(IntervalId::new(2))
                );
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(7.0).unwrap()),
                    Some(IntervalId::new(3))
                );

                // Outside domain
                assert_eq!(
                    grid.try_find_interval_id_of_point(&Real::try_new(-0.1).unwrap()),
                    None
                );
                assert_eq!(
                    grid.try_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
    }

    #[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::Only(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());
                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(),
                );
                assert_eq!(i1, IntervalId::new(0));

                let i2 = grid1d.find_interval_id_of_point(&Real::one());
                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(),
                );
                assert_eq!(i3, IntervalId::new(1));

                let i4 = grid1d.find_interval_id_of_point(&Real::two());
                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(),
                );
                assert_eq!(i5, IntervalId::new(2));

                let i6 = grid1d.find_interval_id_of_point(
                    &Real::try_new(rug::Float::with_val(PRECISION, 3.)).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));
            }
        }
    }
}