grid1d 0.3.6

//! Coordinate collections and domain-validated point sets for 1D grids.
//!
//! This module provides types for representing ordered collections of coordinates
//! and points that are validated to lie within specific domains. These are the
//! fundamental building blocks for grid construction and manipulation.
//!
//! # Core Types
//!
//! | Type | Purpose | Key Feature |
//! |------|---------|-------------|
//! | [`Coords1D`] | Sorted, unique coordinates | O(1) lookup via hash map |
//! | [`Coords1DInDomain`] | Coordinates validated within a domain | Flexible storage (owned, borrowed, Arc, etc.) |
//!
//! # Quick Start
//!
//! ## Creating Coordinates
//!
//! ```rust
//! use grid1d::coords::{Coords1D, Coords1DInDomain};
//! use grid1d::intervals::*;
//! use grid1d::scalars::NumIntervals;
//! use sorted_vec::partial::SortedSet;
//! use try_create::TryNew;
//!
//! // From unsorted values (automatically sorted and deduplicated)
//! let coords = Coords1D::try_from(
//!     SortedSet::from_unsorted(vec![3.0, 1.0, 2.0, 1.0])
//! ).unwrap();
//! assert_eq!(&coords[..], &[1.0, 2.0, 3.0]);
//!
//! // Uniform spacing
//! let domain = IntervalClosed::new(0.0, 10.0);
//! let uniform = Coords1D::new_uniform(&domain, &NumIntervals::try_new(5).unwrap());
//! assert_eq!(uniform.len(), 6); // 5 intervals = 6 points
//! ```
//!
//! ## Domain-Validated Coordinates
//!
//! ```rust
//! use grid1d::coords::{Coords1D, Coords1DInDomain};
//! use grid1d::intervals::*;
//! use sorted_vec::partial::SortedSet;
//!
//! let coords = Coords1D::try_from(
//!     SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])
//! ).unwrap();
//! let domain = IntervalClosed::new(0.0, 1.0);
//!
//! // Validates that all coordinates are within the domain
//! let validated = Coords1DInDomain::try_new(coords, domain).unwrap();
//! ```
//!
//! # Storage Flexibility
//!
//! [`Coords1DInDomain`] supports multiple storage strategies through its generic
//! `CoordStorage` parameter:
//!
//! | Type Alias | Storage | Use Case |
//! |------------|---------|----------|
//! | [`Coords1DInDomainOwned`] | `Coords1D<T>` | Default, full ownership |
//! | [`Coords1DInDomainBorrowed`] | `&Coords1D<T>` | Temporary views |
//! | [`Coords1DInDomainArc`] | `Arc<Coords1D<T>>` | Thread-safe sharing |
//! | [`Coords1DInDomainRc`] | `Rc<Coords1D<T>>` | Single-threaded sharing |
//! | [`Coords1DInDomainBox`] | `Box<Coords1D<T>>` | Heap allocation |
//! | [`Coords1DInDomainCow`] | `Cow<Coords1D<T>>` | Copy-on-write semantics |
//!
//! # Search Operations
//!
//! [`Coords1D`] provides efficient search operations:
//!
//! ```rust
//! use grid1d::coords::Coords1D;
//! use sorted_vec::partial::SortedSet;
//!
//! let coords = Coords1D::try_from(
//!     SortedSet::from_unsorted(vec![10.0, 20.0, 30.0, 40.0])
//! ).unwrap();
//!
//! // O(1) membership test
//! assert!(coords.contains(&20.0));
//! assert!(!coords.contains(&25.0));
//!
//! // O(log n) floor/ceiling operations
//! assert_eq!(coords.find_floor_index(&25.0), Some(1)); // 20.0 at index 1
//! assert_eq!(coords.find_ceil_index(&25.0), Some(2));  // 30.0 at index 2
//!
//! // Find closest coordinate
//! let (closest, idx) = coords.find_closest(&27.0);
//! assert_eq!(*closest, 30.0); // 30.0 is closer than 20.0
//! ```
//!
//! # Mathematical Foundations
//!
//! ## Coordinate Invariants
//!
//! All [`Coords1D`] instances maintain these invariants:
//! - **Non-empty**: At least one coordinate
//! - **Sorted**: Strictly ascending order
//! - **Unique**: No duplicate values
//! - **Finite**: All values are finite (when using validated scalar types)
//!
//! ## Domain Validation
//!
//! [`Coords1DInDomain`] ensures:
//! - **Lower bound**: `min(coords) ≥ domain.lower_bound()` (respecting open/closed)
//! - **Upper bound**: `max(coords) ≤ domain.upper_bound()` (respecting open/closed)
//!
//! # Related Modules
//!
//! - [`crate::intervals`]: Domain interval types
//! - [`crate::scalars`]: Type-safe numeric wrappers
//! - [`crate`]: Grid types that use these coordinates

#![deny(rustdoc::broken_intra_doc_links)]

use crate::{
    intervals::{
        GetLowerBoundValue, GetUpperBoundValue, IntervalBoundsRuntime,
        IntervalFinitePositiveLengthTrait, IntervalTrait, bounded::IntervalClosed,
    },
    scalars::{HashableScalar, NumIntervals, PositiveNumPoints1D},
    traits::{HasCoords1D, HasDomain1D},
};
use derive_more::{AsRef, Index};
use num_valid::{
    RealScalar,
    core::errors::capture_backtrace,
    functions::{Pow, Reciprocal},
};
use serde::{Deserialize, Deserializer, Serialize, Serializer};
use sorted_vec::partial::SortedSet;
use std::{backtrace::Backtrace, collections::HashMap, ops::Deref};
use thiserror::Error;
use try_create::TryNew;

/// **INTERNAL** Merges two sorted coordinate arrays with automatic deduplication.
///
/// This function implements an efficient merge algorithm that combines two sorted
/// arrays while removing duplicate values. It's used internally by grid union
/// operations to create unified coordinate sets.
///
/// ## Algorithm
/// Uses a two-pointer merge technique similar to merge sort, but with deduplication:
/// 1. Compare elements from both arrays
/// 2. Add the smaller element to the result
/// 3. If elements are equal, add only once and advance both pointers
/// 4. Continue until all elements are processed
///
/// ## Performance
/// - **Time complexity**: O(n + m) where n, m are the input array sizes
/// - **Space complexity**: O(n + m) for the output
/// - **Optimality**: This is the optimal algorithm for merging sorted arrays
///
/// ## Parameters
/// - `left`: First sorted array
/// - `right`: Second sorted array
///
/// ## Returns
/// A [`SortedSet`] containing all unique elements from both input arrays.
///
/// ## Panics
/// Panics if any coordinate comparison fails (e.g., due to NaN values).
fn merge_sorted_real_with_dedup<T: RealScalar>(left: &[T], right: &[T]) -> SortedSet<T> {
    let l_length = left.len();
    let r_length = right.len();
    let mut out = Vec::with_capacity(l_length + r_length);
    let mut i = 0;
    let mut j = 0;

    while i < l_length && j < r_length {
        match left[i]
            .partial_cmp(&right[j])
            .expect("Invalid coordinate: comparison failed (NaN or non-finite) in grid coordinates")
        {
            std::cmp::Ordering::Less => {
                out.push(left[i].clone());
                i += 1;
            }
            std::cmp::Ordering::Greater => {
                out.push(right[j].clone());
                j += 1;
            }
            std::cmp::Ordering::Equal => {
                // Equal elements - add only once (deduplication) and advance both
                out.push(left[i].clone());
                i += 1;
                j += 1;
            }
        }
    }

    // Add remaining elements
    out.extend_from_slice(&left[i..]);
    out.extend_from_slice(&right[j..]);

    unsafe { SortedSet::from_sorted(out) }
}
//------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------
/// Generates uniformly spaced points over a finite interval with positive length.
///
/// Creates `n+1` equally spaced points over the given domain, where `n` is the number of intervals.
/// The spacing between consecutive points is `δ = domain_length / n`.
///
/// This function is optimized for creating uniform grids and is used internally by [`crate::Grid1DUniform`].
///
/// ## Mathematical Definition
/// For a domain `[a, b]` with `n` intervals:
/// ```text
/// points[k] = a + k * δ  where δ = (b - a) / n, k = 0, 1, ..., n
/// ```
///
/// ## Performance
/// - **Time complexity**: O(n)
/// - **Space complexity**: O(n)
/// - **Numerical stability**: Uses exact arithmetic for the endpoint calculations
///
/// ## Parameters
/// - `domain`: An interval with finite positive length
/// - `num_intervals`: Number of subintervals to create
///
/// ## Returns
/// A vector of `n+1` coordinates, where the first point equals the domain's lower bound
/// and the last point equals the domain's upper bound.
///
/// ## Example
/// ```rust
/// use grid1d::{linspace, intervals::*, scalars::NumIntervals};
/// use try_create::TryNew;
///
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let coords = linspace(&domain, &NumIntervals::try_new(4).unwrap());
/// assert_eq!(coords, vec![0.0, 0.25, 0.5, 0.75, 1.0]);
/// ```
///
/// ## Numerical Properties
/// - **Endpoint preservation**: First and last points exactly match domain bounds
/// - **Uniform spacing**: All intervals have identical length `δ`
/// - **No accumulated error**: Each point calculated independently to avoid drift
pub fn linspace<IntervalType: IntervalFinitePositiveLengthTrait>(
    domain: &IntervalType,
    num_intervals: &NumIntervals,
) -> Vec<IntervalType::RealType> {
    let n_intervals = *num_intervals.as_ref();
    let n_points = n_intervals + 1;
    let delta_points = domain.length().into_inner()
        * IntervalType::RealType::try_from_f64(n_intervals as f64)
            .unwrap()
            .try_reciprocal()
            .unwrap();

    let lower_bound = domain.lower_bound_value();

    let mut coords = Vec::<_>::with_capacity(n_points);
    coords.push(lower_bound.clone());
    for i in 1..n_intervals {
        let x =
            IntervalType::RealType::try_from_f64(i as f64).unwrap() * &delta_points + lower_bound;
        coords.push(x);
    }
    coords.push(domain.upper_bound_value().clone());

    coords
}

/// Generates logarithmically spaced points in the exponent domain.
///
/// Creates `n+1` points where the logarithms are uniformly distributed over the exponent domain.
/// This is useful for creating grids with logarithmic scaling, commonly used in scientific
/// applications where quantities span multiple orders of magnitude.
///
/// ## Mathematical Definition
/// For an exponent domain `[exp_min, exp_max]` with base `b`:
/// ```text
/// exponents[k] = exp_min + k * δ  where δ = (exp_max - exp_min) / n
/// points[k] = b^exponents[k]
/// ```
///
/// ## Parameters
/// - `base`: The logarithmic base (commonly 10.0 or e)
/// - `exponent_domain`: Interval defining the range of exponents
/// - `num_intervals`: Number of subintervals in the exponent space
///
/// ## Returns
/// A vector of `n+1` coordinates with logarithmic spacing in the original space.
///
/// ## Example
/// ```rust
/// use grid1d::{logspace, intervals::*, scalars::NumIntervals};
/// use try_create::TryNew;
///
/// let exp_domain = IntervalClosed::new(0.0, 2.0); // 10^0 to 10^2
/// let coords = logspace(&10.0, &exp_domain, &NumIntervals::try_new(2).unwrap());
/// assert_eq!(coords, vec![1.0, 10.0, 100.0]);
/// ```
///
/// ## Use Cases
/// - **Frequency analysis**: Creating logarithmic frequency grids
/// - **Scientific data**: Spanning multiple orders of magnitude
/// - **Geometric progression**: When ratios between points should be constant
pub fn logspace<IntervalType: IntervalFinitePositiveLengthTrait>(
    base: &IntervalType::RealType,
    exponent_domain: &IntervalType,
    num_intervals: &NumIntervals,
) -> Vec<IntervalType::RealType> {
    let exponents = linspace(exponent_domain, num_intervals);
    // Convert the exponents to the base raised to the power of each exponent
    let out: Vec<_> = exponents
        .iter()
        .map(|v: &IntervalType::RealType| base.clone().try_pow(v).unwrap())
        .collect();

    out
}
//------------------------------------------------------------------------------------------------

//-------------------------------------------------------------------------------
/// Enum representing errors that can occur when constructing or manipulating one-dimensional coordinates.
///
/// The [`ErrorsCoords1D`] enum provides error types for operations involving the [`Coords1D`] struct, which represents
/// a collection of unique, sorted, one-dimensional coordinates. These errors typically arise when the input data
/// violates the constraints of the `Coords1D` struct, such as being empty.
///
/// # Variants
///
/// - `EmptyObject`: Indicates that the input coordinates are empty.
///
/// # Examples
///
/// ## Handling an Empty Input
///
/// ```rust
/// use grid1d::{Coords1D, ErrorsCoords1D};
/// use sorted_vec::partial::SortedSet;
///
/// let empty_set = SortedSet::<f64>::new();
/// let result = Coords1D::try_from(empty_set);
///
/// match result {
///     Err(ErrorsCoords1D::EmptyObject { .. }) => {
///         println!("Error: The input coordinates are empty!");
///     }
///     _ => panic!("Expected an empty object error."),
/// }
/// ```
#[derive(Debug, Error)]
pub enum ErrorsCoords1D {
    /// Indicates that the input coordinates are empty.
    #[error("The input coordinates are empty!")]
    EmptyObject {
        /// Captured backtrace for debugging the source of the empty coordinate set.
        backtrace: Backtrace,
    },
}

//-------------------------------------------------------------------------------
/// Enum representing errors that can occur when constructing or manipulating a collection of general 1D points.
///
/// The `ErrorsCoords1DInDomain` enum provides error types for operations involving the [`Coords1DInDomain`] struct,
/// which represents a collection of sorted, one-dimensional points within a closed interval. These errors typically
/// arise when the input data violates the constraints of the [`Coords1DInDomain`] struct, such as points lying outside
/// the specified domain.
///
/// # Variants
///
/// - `LowerBoundViolation`: Indicates that the smallest coordinate is less than the lower bound of the domain.
/// - `UpperBoundViolation`: Indicates that the largest coordinate is greater than the upper bound of the domain.
///
/// # Examples
///
/// ## Handling Points Outside the Domain
///
/// ```rust
/// use grid1d::{
///     coords::{Coords1D, Coords1DInDomain, ErrorsCoords1DInDomain,},
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 1.5, 2.0])).unwrap();
/// let domain = IntervalClosed::new(0.0, 1.0);
///
/// let result = Coords1DInDomain::try_new(coords, domain);
///
/// match result {
///     Err(ErrorsCoords1DInDomain::UpperBoundViolation { coord_max, domain, .. }) => {
///         println!("Error: The maximum coordinate ({}) exceeds the domain {:?}", coord_max, domain);
///     }
///     Err(ErrorsCoords1DInDomain::LowerBoundViolation { coord_min, domain, .. }) => {
///         println!("Error: The minimum coordinate ({}) is below the domain {:?}", coord_min, domain);
///     }
///     _ => panic!("Expected a domain violation error."),
/// }
/// ```
#[derive(Debug, Error)]
pub enum ErrorsCoords1DInDomain<IntervalType: IntervalTrait> {
    /// Error indicating that the minimum coordinate violates the domain's lower bound.
    #[error("The minimum coordinate of the points is smaller than the lower bound of the domain.")]
    LowerBoundViolation {
        /// The minimum coordinate value found in the point set.
        coord_min: IntervalType::RealType,
        /// The domain interval that was violated.
        domain: IntervalType,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },

    /// Error indicating that the maximum coordinate violates the domain's upper bound.
    #[error("The maximum coordinate of the points is greater than the upper bound of the domain.")]
    UpperBoundViolation {
        /// The maximum coordinate value found in the point set.
        coord_max: IntervalType::RealType,
        /// The domain interval that was violated.
        domain: IntervalType,
        /// Captured backtrace for debugging.
        backtrace: Backtrace,
    },
}

/// Container for one-dimensional coordinates that are ***unique***, ***sorted***, and ***non-empty***.
///
/// The [`Coords1D`] struct represents a collection of ***unique***, ***sorted*** (in ascending order),
/// one-dimensional coordinates that is guaranteed to never be empty. This container is fundamental
/// to numerical computations, serving as the backbone for defining grids, partitions, interpolation
/// points, and sets of coordinates in one-dimensional space.
///
/// ## Core Mathematical Properties
///
/// - **Uniqueness**: No duplicate coordinates are allowed
/// - **Ordering**: Coordinates are always sorted in strict ascending order
/// - **Non-emptiness**: Always contains at least one coordinate
/// - **Immutability**: Once created, the coordinate set cannot be modified
/// - **Generic**: Works with any type implementing [`RealScalar`]
///
/// ## Generic Design with [`num_valid::RealScalar`]
///
/// **All coordinate operations work seamlessly with any scalar type implementing [`num_valid::RealScalar`]:**
///
/// | Scalar Type | Performance | Validation | Best For |
/// |-------------|-------------|------------|----------|
/// | [`f64`] | ⚡⚡⚡ Maximum | ❌ None | Trusted input, absolute maximum speed |
/// | [`RealNative64StrictFiniteInDebug`](num_valid::RealNative64StrictFiniteInDebug) | ⚡⚡⚡ **Same as f64** | ✅ **Debug only** | **Recommended for most applications** |
/// | [`RealNative64StrictFinite`](num_valid::RealNative64StrictFinite) | ⚡⚡ Small overhead | ✅ Always | Safety-critical applications |
/// | `RealRugStrictFinite` | ⚡ Precision-dependent | ✅ Always | Arbitrary precision needs (available from the [num-valid](https://crates.io/crates/num-valid) crate when compiled with `--features=rug`)|
///
/// ## Design Principles
///
/// ### Mathematical Correctness First
/// - **Invariant preservation**: Coordinates remain unique, sorted, and non-empty throughout lifetime
/// - **Type safety**: Impossible to create invalid coordinate sets through safe constructors
/// - **Numerical stability**: Compatible with validated scalar types to prevent NaN/infinity issues
///
/// ### Zero-Cost Abstractions
/// - **Minimal overhead**: Thin wrapper around [`SortedSet`] with compile-time guarantees
/// - **Direct access**: Zero-cost conversion to slices for performance-critical operations
/// - **Memory efficient**: No additional allocations beyond the underlying array
///
/// ### Integration-Focused Design
/// - **Grid construction**: Primary input for all grid and partition types in the library
/// - **Interpolation**: Direct compatibility with numerical interpolation algorithms
/// - **Conversion**: Seamless interop with standard collections
///
/// ## Performance Characteristics
///
/// | Operation | Time Complexity | Space Complexity | Notes |
/// |-----------|----------------|------------------|-------|
/// | **Construction** | O(n log n) | O(n) | Due to sorting/deduplication in conversion |
/// | **Access by index** | O(1) | O(1) | Direct array indexing |
/// | **Slice access** | O(1) | O(1) | Zero-cost conversion to `&[T]` |
/// | **Search operations** | O(log n) | O(1) | Binary search on sorted data |
/// | **Iteration** | O(n) | O(1) | Cache-friendly sequential access |
/// | **Length query** | O(1) | O(1) | Cached in the structure |
///
/// ## Construction and Usage Examples
///
/// ### Recommended Construction Pattern
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// // Define your optimal type alias
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // From pre-sorted unique coordinates
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![
///         Real::try_new(0.0).unwrap(),
///         Real::try_new(1.0).unwrap(),
///         Real::try_new(2.0).unwrap()
///     ])
/// ).unwrap();
///
/// assert_eq!(coords.len(), 3);
/// assert_eq!(coords.num_points().as_ref(), &3);
/// ```
///
/// ### Uniform Coordinate Generation
/// ```rust
/// use grid1d::{
///     Coords1D,
///     intervals::*,
///     scalars::NumIntervals,
/// };
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// type Real = RealNative64StrictFiniteInDebug;
///
/// let domain = IntervalClosed::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let coords = Coords1D::new_uniform(&domain, &NumIntervals::try_new(4).unwrap());
///
/// // Results in coordinates: [0.0, 0.25, 0.5, 0.75, 1.0]
/// assert_eq!(coords.len(), 5);
/// ```
///
/// ### Different Scalar Types for Different Needs
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
/// use num_valid::{RealNative64StrictFiniteInDebug, RealNative64StrictFinite};
/// use try_create::TryNew;
///
/// // Performance-optimal with debug safety (recommended)
/// type OptimalReal = RealNative64StrictFiniteInDebug;
/// let optimal_coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![
///         OptimalReal::try_new(0.0).unwrap(),
///         OptimalReal::try_new(1.0).unwrap()
///     ])
/// ).unwrap();
///
/// // Always validated for safety-critical applications
/// type SafeReal = RealNative64StrictFinite;
/// let safe_coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![
///         SafeReal::try_new(0.0).unwrap(),
///         SafeReal::try_new(1.0).unwrap()
///     ])
/// ).unwrap();
///
/// // Raw f64 for absolute maximum performance (use with caution)
/// let raw_coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 1.0])
/// ).unwrap();
///
/// // All have identical APIs and mathematical behavior
/// assert_eq!(optimal_coords.num_points(), safe_coords.num_points());
/// assert_eq!(safe_coords.num_points(), raw_coords.num_points());
/// ```
///
/// ## Integration with Grid Construction
///
/// ### Basic Grid Creation
/// ```rust
/// use grid1d::{
///     Coords1D, Grid1D, IntervalPartition,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 2.5, 5.0, 7.5, 10.0])
/// ).unwrap();
///
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_coords(coords).unwrap();
/// assert_eq!(grid.num_intervals().as_ref(), &4);
/// ```
///
/// ### Interpolation and Numerical Methods
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
///
/// let x_values = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 1.0, 4.0, 9.0, 16.0])
/// ).unwrap();
///
/// // Direct access for numerical computations
/// let y_values: Vec<_> = x_values.iter().map(|&x| x.sqrt()).collect();
///
/// // Efficient vectorized operations
/// let sum: f64 = x_values.iter().sum();
/// let mean = sum / x_values.len() as f64;
/// assert!((mean - 6.0).abs() < 1e-10);
/// ```
///
/// ## Access Patterns and Performance
///
/// ### High-Performance Access
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted((0..1000).map(|i| i as f64).collect())
/// ).unwrap();
///
/// // Zero-cost slice access for performance-critical loops
/// let mut sum = 0.0;
/// for value in coords.iter() {
///     sum += value;
/// }
///
/// let sum_2 = coords.into_iter().sum::<f64>();
/// assert!((sum - sum_2).abs() < 1e-10);
/// ```
///
/// ### Search and Lookup Operations
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(
///     unsafe { SortedSet::from_sorted(vec![1.0, 3.0, 5.0, 7.0, 9.0]) }
/// ).unwrap();
///
/// // Using the internal lookup index table (O(1))
/// let search_value = 5.0;
/// let coord_id = match coords.find_index(&search_value) {
///     Some(index) => {
///         println!("Found {} at index {}", search_value, index);
///         Some(index)
///     }
///     None => {
///         println!("Index {} not found", search_value);
///         None
///     }
/// };
/// assert_eq!(coord_id, Some(2));
///
/// // Binary search (O(log n))
/// let search_value = 5.0;
/// let coord_id = match coords.binary_search_by(|probe| probe.partial_cmp(&search_value).unwrap()) {
///     Ok(index) => {
///         println!("Found {} at index {}", search_value, index);
///         Ok(index)
///     }
///     Err(index) => {
///         println!("Would insert {} at index {}", search_value, index);
///         Err(index)
///     }
/// };
/// assert_eq!(coord_id, Ok(2));
///
/// // Find insertion point for maintaining sorted order
/// let new_value = 6.0;
/// let insertion_index = coords.find_insertion_index(&new_value);
/// println!("Insert {} at index {} to maintain order", new_value, insertion_index);
/// assert_eq!(insertion_index, 3);
/// ```
///
/// ## Error Handling and Validation
///
/// ### Construction Error Handling
/// ```rust
/// use grid1d::{Coords1D, ErrorsCoords1D};
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// // Handle empty coordinate sets
/// let empty_set = SortedSet::<f64>::new();
/// match Coords1D::try_from(empty_set) {
///     Ok(_) => unreachable!("Should not succeed with empty coordinates"),
///     Err(ErrorsCoords1D::EmptyObject { .. }) => {
///         println!("Correctly rejected empty coordinate set");
///     }
/// }
///
/// // Successful construction with automatic deduplication
/// let coords_with_duplicates = SortedSet::from_unsorted(vec![1.0, 2.0, 2.0, 3.0]);
/// let coords = Coords1D::try_from(coords_with_duplicates).unwrap();
/// assert_eq!(coords.deref(), &[1.0, 2.0, 3.0]); // Duplicates automatically removed
/// ```
///
/// ### Robust Construction Pattern
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
///
/// fn safe_coords_from_vec(values: Vec<f64>) -> Result<Coords1D<f64>, String> {
///     // Validate input
///     if values.is_empty() {
///         return Err("Cannot create coordinates from empty vector".to_string());
///     }
///     
///     // Check for invalid values
///     for (i, &value) in values.iter().enumerate() {
///         if !value.is_finite() {
///             return Err(format!("Invalid value at index {}: {}", i, value));
///         }
///     }
///     
///     // Create sorted set (handles deduplication and sorting)
///     let sorted_set = SortedSet::from_unsorted(values);
///     
///     // Convert to Coords1D
///     Coords1D::try_from(sorted_set)
///         .map_err(|e| format!("Failed to create coordinates: {}", e))
/// }
///
/// // Usage
/// assert!(safe_coords_from_vec(vec![]).is_err());
/// assert!(safe_coords_from_vec(vec![f64::NAN]).is_err());
/// assert!(safe_coords_from_vec(vec![1.0, 2.0, 3.0]).is_ok());
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Generic Programming Over Coordinate Types
/// ```rust
/// use grid1d::{Coords1D, HasCoords1D};
/// use num_valid::RealScalar;
///
/// fn analyze_coordinates<T>(coords: &Coords1D<T>) -> String
/// where
///     T: RealScalar + std::fmt::Display,
/// {
///     let first = coords.first();
///     let last = coords.last();
///     let count = coords.num_points();
///     
///     format!(
///         "Coordinates: {} points from {} to {} (span: {})",
///         count.as_ref(),
///         first,
///         last,
///         last.clone() - first
///     )
/// }
///
/// let coords = Coords1D::try_from(
///     sorted_vec::partial::SortedSet::from_unsorted(vec![0.0, 2.5, 5.0])
/// ).unwrap();
///
/// println!("{}", analyze_coordinates(&coords));
/// // Output: "Coordinates: 3 points from 0 to 5 (span: 5)"
/// ```
///
/// ### Coordinate Transformation and Mapping
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// let original_coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 1.0, 2.0, 3.0, 4.0])
/// ).unwrap();
///
/// // Transform coordinates (maintaining order is critical)
/// let transformed_values: Vec<f64> = original_coords
///     .iter()
///     .map(|&x| x * x) // Square transformation
///     .collect();
///
/// let transformed_coords = Coords1D::try_from(
///     SortedSet::from_unsorted(transformed_values)
/// ).unwrap();
///
/// assert_eq!(transformed_coords.deref(), &[0.0, 1.0, 4.0, 9.0, 16.0]);
/// ```
///
/// ## Memory Layout and Optimization
///
/// ### Cache-Friendly Iteration
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
///
/// let large_coords = Coords1D::try_from(
///     SortedSet::from_unsorted((0..100_000).map(|i| i as f64).collect())
/// ).unwrap();
///
/// // Sequential access is cache-friendly
/// let mut running_sum = 0.0;
/// for &value in large_coords.iter() {
///     running_sum += value; // Efficient sequential memory access
/// }
///
/// // Random access patterns should be avoided for large coordinate sets
/// // when performance is critical
/// ```
///
/// ## Integration with the other numerical Ecosystem
///
/// ### Grid Construction Workflow
/// ```rust
/// use grid1d::{
///     Coords1D, Grid1D, IntervalPartition,
///     intervals::*,
///     scalars::IntervalId,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // 1. Create coordinates
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 0.3, 0.7, 1.0])
/// ).unwrap();
///
// 2. Build grid from coordinates
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_coords(coords).unwrap();
///
/// // 3. Use grid for interval operations
/// let point = 0.5;
/// let interval_id = grid.find_interval_id_of_point(&point);
/// let interval = grid.interval(&interval_id);
///
/// println!("Point {} is in interval {:?}", point, interval);
/// ```
///
/// ### Multi-Physics Grid Union
/// ```rust
/// use grid1d::{
///     Coords1D, Grid1D, Grid1DUnion, IntervalPartition,
///     intervals::*,
///  };
/// use sorted_vec::partial::SortedSet;
///
/// // Physics A: coarse grid
/// let coords_a = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])
/// ).unwrap();
/// let grid_a = Grid1D::<IntervalClosed<f64>>::try_from_coords(coords_a).unwrap();
///
/// // Physics B: fine grid in specific region
/// let coords_b = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 0.5, 1.0, 1.5, 2.0])
/// ).unwrap();
/// let grid_b = Grid1D::<IntervalClosed<f64>>::try_from_coords(coords_b).unwrap();
///
/// // Union creates refined grid
/// let union = Grid1DUnion::try_new(&grid_a, &grid_b).unwrap();
///
/// println!("Unified grid has {} intervals", union.num_intervals().as_ref());
/// assert_eq!(union.num_intervals().as_ref(), &4);
/// // The unified grid coordinates automatically merge both physics requirements
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`Coords1D`]
/// - **Grid construction**: Primary use case for all 1D grids and partitions
/// - **Interpolation nodes**: Store points for numerical interpolation
/// - **Integration quadrature**: Store quadrature points for numerical integration
/// - **Measurement data**: Store sorted measurement locations or time points
/// - **Boundary conditions**: Store discrete boundary point locations
///
/// ### When NOT to Use [`Coords1D`]
/// - **Unsorted data**: Use regular `Vec<T>` if order doesn't matter
/// - **Data with duplicates**: Use `Vec<T>` if duplicates are meaningful
/// - **Frequently modified sets**: Use `Vec<T>` if you need to add/remove points often
/// - **Multidimensional coordinates**: Use appropriate multidimensional containers
///
/// ### Recommended Patterns
/// ```rust
/// use grid1d::Coords1D;
/// use sorted_vec::partial::SortedSet;
/// use num_valid::{RealNative64StrictFiniteInDebug, RealScalar};
/// use num::Zero;
/// use try_create::TryNew;
///
/// // Define type aliases for your application
/// type Real = RealNative64StrictFiniteInDebug;
/// type Coordinates = Coords1D<Real>;
///
/// // Pattern 1: Validate before construction
/// fn create_validated_coords(values: Vec<f64>) -> Result<Coordinates, String> {
///     // Convert to Real type with validation
///     let real_values: Result<Vec<Real>, _> = values
///         .into_iter()
///         .map(|v| Real::try_new(v))
///         .collect();
///     
///     let real_values = real_values.map_err(|e| format!("Invalid value: {}", e))?;
///     
///     // Create sorted set and coords
///     let sorted_set = SortedSet::from_unsorted(real_values);
///     Coordinates::try_from(sorted_set)
///         .map_err(|e| format!("Coordinate creation failed: {}", e))
/// }
///
/// // Pattern 2: Extract for high-performance computation
/// fn compute_on_coords(coords: &Coordinates) -> Real {
///     let n_as_real = Real::try_from_f64(coords.len() as f64).unwrap();
///
///     coords.iter().cloned().fold(Real::zero(), |acc, x| acc + x) / &n_as_real
/// }
///
/// // Pattern 3: Convert between coordinate types when needed
/// fn convert_coords_precision(coords: Coords1D<f64>) -> Coords1D<Real> {
///     let values: Vec<Real> = coords
///         .iter()
///         .map(|&x| Real::try_new(x).unwrap())
///         .collect();
///     
///     Coords1D::try_from(SortedSet::from_unsorted(values)).unwrap()
/// }
/// ```
///
/// ## Mathematical Properties and Guarantees
///
/// The [`Coords1D`] struct maintains these mathematical properties:
///
/// - **Strict ordering**: `coords[i] < coords[i+1]` for all valid indices `i`
/// - **Uniqueness**: `coords[i] ≠ coords[j]` for all `i ≠ j`
/// - **Non-emptiness**: `coords.len() ≥ 1` always holds
/// - **Finite values**: All coordinates are finite real numbers (when using validated scalar types)
/// - **Index consistency**: `coords.num_points().as_ref() == coords.len()`
/// - **Immutability**: Once created, the coordinate set cannot be modified
///
/// These properties are enforced at construction time and preserved throughout the lifetime
/// of the [`Coords1D`] instance, making it safe to use in mathematical algorithms that
/// depend on these invariants.
#[derive(Debug, Clone, AsRef, Index)]
pub struct Coords1D<RealType: RealScalar> {
    #[as_ref]
    #[index]
    /// The coordinates stored as a [`SortedSet`] object.
    /// These coordinates are stored in ascending order and do not contains duplicated values.
    values: SortedSet<RealType>,

    /// A lookup table mapping each coordinate to its index in the `values` array.
    /// This allows for O(1) retrieval of the index for a given coordinate.
    /// The keys are wrapped in `HashableScalar` to ensure proper hashing and equality checks.
    lookup_index: HashMap<HashableScalar<RealType>, usize>,
}

impl<RealType: RealScalar> PartialEq for Coords1D<RealType> {
    fn eq(&self, other: &Self) -> bool {
        self.values == other.values
    }
}

impl<RealType: RealScalar> Deref for Coords1D<RealType> {
    type Target = [RealType];

    /// Allows `Coords1D` to be treated as a slice `&[RealType]`.
    ///
    /// This provides ergonomic access to all slice methods (`len`, `iter`, `first`, `binary_search`, etc.)
    /// while preventing any mutation, thus preserving the struct's invariants (sorted, unique, non-empty).
    fn deref(&self) -> &Self::Target {
        self.values.as_slice()
    }
}

impl<RealType: RealScalar> Coords1D<RealType> {
    #[inline(always)]
    fn new(coords: SortedSet<RealType>) -> Self {
        let lookup_index = coords
            .iter()
            .cloned()
            .enumerate()
            .map(|(i, v)| (HashableScalar::new(v), i))
            .collect::<HashMap<HashableScalar<RealType>, usize>>();

        Self {
            values: coords,
            lookup_index,
        }
    }

    /// Returns the underlying coordinates as a `Vec<RealType>`, consuming the `Coords1D` instance.
    /// This is a zero-cost operation that transfers ownership of the internal vector.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![2.0, 1.0, 3.0])).unwrap();
    /// let vec = coords.into_vec();
    /// assert_eq!(vec, vec![1.0, 2.0, 3.0]);
    /// ```
    ///
    /// # Note
    /// `to_vec()` is a slice method that is accessible through deref and creates a clone of the data,
    /// while `into_vec()` consumes the [`Coords1D`] instance and transfers ownership of the internal vector.
    #[inline(always)]
    pub fn into_vec(self) -> Vec<RealType> {
        self.values.into_vec()
    }

    /// Returns the number of coordinates in the [`Coords1D`] object.
    #[inline(always)]
    pub fn num_points(&self) -> PositiveNumPoints1D {
        PositiveNumPoints1D::try_new(self.values.len()).unwrap()
    }

    /// Creates a new [`Coords1D`] instance representing uniformly spaced points over a bounded interval.
    ///
    /// This function generates a set of coordinates that are evenly distributed within the specified domain.
    /// The first point coincides with the lower bound of the interval, and the last point coincides with the upper bound.
    /// The number of points generated is equal to `num_intervals + 1`.
    ///
    /// # Parameters
    ///
    /// - `domain`: The interval over which the points are distributed.
    /// - `num_intervals`: The number of uniform subintervals to divide the domain into. The resulting number of points will be `num_intervals + 1`.
    ///
    /// # Returns
    ///
    /// A [`Coords1D`] object containing the uniformly spaced coordinates.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Coords1D, intervals::*,
    ///     scalars::NumIntervals,
    /// };
    /// use std::ops::Deref;
    /// use try_create::TryNew;
    ///
    /// let domain = IntervalClosed::new(0.0, 1.0);
    /// let coords = Coords1D::new_uniform(&domain, &NumIntervals::try_new(4).unwrap());
    ///
    /// assert_eq!(coords.deref(), &[0.0, 0.25, 0.5, 0.75, 1.0]);
    /// ```
    #[inline(always)]
    pub fn new_uniform(
        domain: &impl IntervalFinitePositiveLengthTrait<RealType = RealType>,
        num_intervals: &NumIntervals,
    ) -> Self {
        let coords = linspace(domain, num_intervals);
        let coords = unsafe { SortedSet::from_sorted(coords) };

        Self::new(coords)
    }

    /// Returns the first (minimum) coordinate.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 2.0])).unwrap();
    /// assert_eq!(coords.first(), &1.0);
    /// ```
    #[inline(always)]
    pub fn first(&self) -> &RealType {
        &self[0]
    }

    /// Returns the last (maximum) coordinate.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 2.0])).unwrap();
    /// assert_eq!(coords.last(), &3.0);
    /// ```
    #[inline(always)]
    pub fn last(&self) -> &RealType {
        &self[self.len() - 1]
    }

    /// Returns the span (range) of the coordinates.
    ///
    /// This is the difference between the maximum and minimum coordinates.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 5.0])).unwrap();
    /// assert_eq!(coords.span(), 4.0); // 5.0 - 1.0
    /// ```
    #[inline(always)]
    pub fn span(&self) -> RealType {
        self.last().clone() - self.first()
    }

    /// Checks if the coordinates contain the given value.
    ///
    /// Uses the internal lookup index for O(1) efficiency.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(
    ///     SortedSet::from_unsorted(vec![1.0, 2.0, 3.0])
    /// ).unwrap();
    /// assert!(coords.contains(&2.0));
    /// assert!(!coords.contains(&2.5));
    /// ```
    #[inline(always)]
    pub fn contains(&self, value: &RealType) -> bool {
        self.lookup_index
            .contains_key(&HashableScalar::new(value.clone()))
    }

    /// Finds the index of the given coordinate `value`.
    ///
    /// Uses the internal lookup index for O(1) efficiency.
    ///
    /// # Returns
    /// The index of the coordinate if found, or `None` if not present.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(
    ///     SortedSet::from_unsorted(vec![1.0, 2.0, 3.0])
    /// ).unwrap();
    /// assert_eq!(coords.find_index(&2.0), Some(1));
    /// assert_eq!(coords.find_index(&2.5), None);
    /// ```
    #[inline(always)]
    pub fn find_index(&self, value: &RealType) -> Option<usize> {
        self.lookup_index
            .get(&HashableScalar::new(value.clone()))
            .copied()
    }

    /// Finds the index where the given `value` would be inserted to maintain sorted order.
    ///
    /// If the `value` already exists, returns the index of the existing value.
    /// Uses binary search for O(log n) efficiency.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 5.0])).unwrap();
    /// assert_eq!(coords.find_insertion_index(&2.0), 1); // Insert between 1.0 and 3.0
    /// assert_eq!(coords.find_insertion_index(&3.0), 1); // Existing value
    /// assert_eq!(coords.find_insertion_index(&6.0), 3); // Insert at end
    /// assert_eq!(coords.find_insertion_index(&0.0), 0); // Insert at beginning
    /// ```
    #[inline(always)]
    pub fn find_insertion_index(&self, value: &RealType) -> usize {
        self.partition_point(|coord| coord < value)
    }

    /// Finds the index of the largest coordinate that is less than or equal to the given `value``.
    ///
    /// This is useful for interpolation and interval finding operations.
    ///
    /// # Returns
    /// The index of the largest coordinate ≤ `value`, or `None` if all coordinates are greater than `value`.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 5.0])).unwrap();
    /// assert_eq!(coords.find_floor_index(&1.0), Some(0)); // 1.0 is largest ≤ 1.0
    /// assert_eq!(coords.find_floor_index(&2.5), Some(0)); // 1.0 is largest ≤ 2.5
    /// assert_eq!(coords.find_floor_index(&3.0), Some(1)); // 3.0 is largest ≤ 3.0
    /// assert_eq!(coords.find_floor_index(&4.0), Some(1)); // 3.0 is largest ≤ 4.0
    /// assert_eq!(coords.find_floor_index(&0.5), None);    // No coordinate ≤ 0.5
    /// assert_eq!(coords.find_floor_index(&5.0), Some(2)); // 5.0 is largest ≤ 5.0
    /// assert_eq!(coords.find_floor_index(&6.0), Some(2)); // 5.0 is largest ≤ 6.0
    /// ```
    pub fn find_floor_index(&self, value: &RealType) -> Option<usize> {
        // `partition_point` returns the index of the first element `p` such that `p >= value`.
        let insertion_point = self.find_insertion_index(value);

        // Case 1: The value is smaller than the very first coordinate.
        // This happens if `insertion_point` is 0 and the value at index 0 is strictly greater than `value`.
        // In this case, no floor exists.
        if insertion_point == 0 && &self.values[0] > value {
            return None;
        }

        // Case 2: The value is an exact match with a coordinate.
        // `partition_point` gives us the index of the match.
        if insertion_point < self.values.len() && &self.values[insertion_point] == value {
            return Some(insertion_point);
        }

        // Case 3: The value is between two coordinates or larger than the last coordinate.
        // The floor is the coordinate just before the insertion point.
        // This is safe because if `insertion_point` is 0, we've already handled the case
        // where `value` is smaller than the first element. The only remaining possibility
        // for `insertion_point == 0` is an exact match, which is handled above.
        // Therefore, `insertion_point` must be > 0 here.
        more_asserts::debug_assert_gt!(insertion_point, 0, "insertion_point should be > 0 here");
        Some(insertion_point - 1)
    }

    /// Finds the index of the smallest coordinate that is greater than or equal to the given `value`.
    ///
    /// This is useful for interpolation and interval finding operations.
    ///
    /// # Returns
    /// The index of the smallest coordinate ≥ `value`, or `None` if all coordinates are less than `value`.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 5.0])).unwrap();
    /// assert_eq!(coords.find_ceil_index(&2.5), Some(1)); // 3.0 is smallest ≥ 2.5
    /// assert_eq!(coords.find_ceil_index(&3.0), Some(1)); // 3.0 is smallest ≥ 3.0
    /// assert_eq!(coords.find_ceil_index(&0.5), Some(0)); // 1.0 is smallest ≥ 0.5
    /// assert_eq!(coords.find_ceil_index(&6.0), None);    // No coordinate ≥ 6.0
    /// ```
    #[inline(always)]
    pub fn find_ceil_index(&self, value: &RealType) -> Option<usize> {
        let insertion_point = self.find_insertion_index(value);
        if insertion_point < self.values.len() {
            Some(insertion_point)
        } else {
            None
        }
    }

    /// Finds the indices of coordinates within the given range (inclusive).
    ///
    /// # Parameters
    /// - `min_value`: The minimum value of the range (inclusive)
    /// - `max_value`: The maximum value of the range (inclusive)
    ///
    /// # Returns
    /// A range of indices corresponding to coordinates within [min_value, max_value].
    /// Returns an empty range if no coordinates are in the specified range.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::{Coords1D, intervals::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 2.0, 3.0, 4.0, 5.0])).unwrap();
    /// let range = coords.find_range_indices(&IntervalClosed::try_new(2.5, 4.5).unwrap());
    /// let indices: Vec<usize> = range.collect();
    /// assert_eq!(indices, vec![2, 3]); // Coordinates 3.0 and 4.0
    /// ```
    pub fn find_range_indices(
        &self,
        range_interval: &IntervalClosed<RealType>,
    ) -> std::ops::Range<usize> {
        let start_idx = self
            .find_ceil_index(range_interval.lower_bound_value())
            .unwrap_or(self.values.len());
        let end_idx = match self.find_floor_index(range_interval.upper_bound_value()) {
            Some(idx) => idx + 1,
            None => 0,
        };

        start_idx..end_idx.max(start_idx)
    }

    /// Returns an iterator over coordinate values within the given range (inclusive).
    ///
    /// # Parameters
    /// - `min_value`: The minimum value of the range (inclusive)
    /// - `max_value`: The maximum value of the range (inclusive)
    ///
    /// # Example
    /// ```rust
    /// use grid1d::{Coords1D, intervals::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![2.0, 1.0, 3.0, 4.0, 5.0])).unwrap();
    /// let values: Vec<f64> = coords.range(&IntervalClosed::try_new(2.5, 4.5).unwrap()).cloned().collect();
    /// assert_eq!(values, vec![3.0, 4.0]);
    /// ```
    #[inline(always)]
    pub fn range(
        &self,
        range_interval: &IntervalClosed<RealType>,
    ) -> impl Iterator<Item = &RealType> {
        let range_indices = self.find_range_indices(range_interval);
        self[range_indices].iter()
    }

    /// Finds the closest coordinate to the given value.
    ///
    /// # Returns
    /// A tuple containing the closest coordinate and its index.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 3.0, 5.0])).unwrap();
    /// let (closest, index) = coords.find_closest(&3.8);
    /// assert_eq!(*closest, 3.0);
    /// assert_eq!(index, 1);
    /// ```
    pub fn find_closest(&self, value: &RealType) -> (&RealType, usize) {
        let insertion_point = self.find_insertion_index(value);

        if insertion_point == 0 {
            (&self.values[0], 0)
        } else if insertion_point >= self.values.len() {
            let last_idx = self.values.len() - 1;
            (&self.values[last_idx], last_idx)
        } else {
            let left_idx = insertion_point - 1;
            let right_idx = insertion_point;

            let left_dist = (value.clone() - &self.values[left_idx]).abs();
            let right_dist = (self.values[right_idx].clone() - value).abs();

            if left_dist <= right_dist {
                (&self.values[left_idx], left_idx)
            } else {
                (&self.values[right_idx], right_idx)
            }
        }
    }

    /// Computes the union of two [`Coords1D`] instances.
    ///
    /// The resulting [`Coords1D`] will contain all unique points from both sets, sorted.
    /// This is the core operation for [`crate::Grid1DUnion`].
    ///
    /// # Performance
    /// O(n + m), where n and m are the number of points in the respective sets.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    /// use std::ops::Deref;
    ///
    /// let coords1 = Coords1D::try_from(SortedSet::from_unsorted(vec![2.0, 1.0, 3.0])).unwrap();
    /// let coords2 = Coords1D::try_from(SortedSet::from_unsorted(vec![5.0, 4.0, 3.0])).unwrap();
    /// let union_coords = coords1.union(&coords2);
    /// assert_eq!(union_coords.deref(), &[1.0, 2.0, 3.0, 4.0, 5.0]);
    /// ```
    #[inline(always)]
    pub fn union(&self, other: &Self) -> Self {
        let merged = merge_sorted_real_with_dedup(self, other);
        Self::new(merged)
    }

    /// Computes the intersection of two [`Coords1D`] instances.
    ///
    /// The resulting [`Coords1D`] will contain only the points that are present in both sets.
    /// Returns `None` if the intersection is empty.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    /// use std::ops::Deref;
    ///
    /// let coords1 = Coords1D::try_from(SortedSet::from_unsorted(vec![2.0, 1.0, 3.0])).unwrap();
    /// let coords2 = Coords1D::try_from(SortedSet::from_unsorted(vec![5.0, 4.0, 3.0])).unwrap();
    /// let intersection_coords = coords1.intersection(&coords2).unwrap();
    /// assert_eq!(intersection_coords.deref(), &[3.0]);
    /// ```
    pub fn intersection(&self, other: &Self) -> Option<Self> {
        let len_self = self.len();
        let len_other = other.len();

        let mut intersection_vec = Vec::with_capacity(len_self.min(len_other));
        let mut i = 0;
        let mut j = 0;

        while i < self.len() && j < other.len() {
            match self[i].partial_cmp(&other[j]).unwrap() {
                std::cmp::Ordering::Less => i += 1,
                std::cmp::Ordering::Greater => j += 1,
                std::cmp::Ordering::Equal => {
                    intersection_vec.push(self[i].clone());
                    i += 1;
                    j += 1;
                }
            }
        }

        if intersection_vec.is_empty() {
            None
        } else {
            Some(Self::new(unsafe {
                SortedSet::from_sorted(intersection_vec)
            }))
        }
    }

    /// Checks if `self` is a subset of `other`.
    ///
    /// Returns `true` if every coordinate in `self` is also present in `other`.
    ///
    /// # Example
    /// ```rust
    /// use grid1d::Coords1D;
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords1 = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 2.0])).unwrap();
    /// let coords2 = Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 2.0, 3.0])).unwrap();
    /// assert!(coords1.is_subset(&coords2));
    /// assert!(!coords2.is_subset(&coords1));
    /// ```
    #[inline(always)]
    pub fn is_subset(&self, other: &Self) -> bool {
        if self.len() > other.len() {
            return false;
        }
        self.iter().all(|val| other.contains(val))
    }
}

impl<RealType: RealScalar> TryFrom<SortedSet<RealType>> for Coords1D<RealType> {
    type Error = ErrorsCoords1D;

    #[inline(always)]
    fn try_from(coords: SortedSet<RealType>) -> Result<Self, Self::Error> {
        if coords.is_empty() {
            Err(ErrorsCoords1D::EmptyObject {
                backtrace: capture_backtrace(),
            })
        } else {
            Ok(Self::new(coords))
        }
    }
}

// Serialize implementation: serialize as a Vec<RealType>
impl<RealType: RealScalar> Serialize for Coords1D<RealType> {
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: Serializer,
    {
        self.values.as_slice().serialize(serializer)
    }
}

// Deserialize implementation: deserialize from Vec<RealType> and reconstruct
impl<'de, RealType: RealScalar> Deserialize<'de> for Coords1D<RealType> {
    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
    where
        D: Deserializer<'de>,
    {
        let vec = Vec::<RealType>::deserialize(deserializer)?;
        let sorted_set = SortedSet::from_unsorted(vec);
        Self::try_from(sorted_set)
            .map_err(|e| serde::de::Error::custom(format!("Failed to create Coords1D: {}", e)))
    }
}
//-------------------------------------------------------------------------------

///  A validated collection of one-dimensional coordinates guaranteed to lie within a specific interval domain.
///
/// The [`Coords1DInDomain`] struct represents a set of one-dimensional points that are mathematically
/// guaranteed to be contained within a specified interval domain. Unlike [`Coords1D`] which only
/// ensures uniqueness, sorting, and non-emptiness, [`Coords1DInDomain`] adds the critical constraint
/// that all coordinates must lie within the bounds of its associated interval.
///
/// This struct serves as a **domain-validated coordinate container**, providing **runtime validation**
/// and **post-construction guarantees** that coordinates and their domain are mathematically consistent.
///
/// ## Core Mathematical Properties
///
/// - **Domain Containment**: Every coordinate `x[i]` satisfies `domain.contains_point(&x[i])` (validated at construction)
/// - **Uniqueness**: No duplicate coordinates (inherited from [`Coords1D`])
/// - **Ordering**: Coordinates are sorted in strict ascending order (inherited from [`Coords1D`])
/// - **Non-emptiness**: Always contains at least one coordinate (inherited from [`Coords1D`])
/// - **Type Safety**: The interval type is encoded in the struct's type parameter (compile-time)
/// - **Immutability**: Once created, neither coordinates nor domain can be modified (post-construction guarantee)
///
/// ## Design Philosophy
///
/// ### Mathematical Correctness Through Runtime Validation
/// ```text
/// Invariant: ∀i, domain.contains_point(&coords[i]) = true
/// ```
/// This invariant is **established during construction through runtime validation** and maintained
/// throughout the object's lifetime through immutability, preventing entire classes of domain-related
/// bugs in numerical algorithms.
///
/// ### Domain-Specific Validation
/// The **runtime validation** respects the mathematical semantics of different interval types:
/// - **Closed intervals `[a,b]`**: Both endpoints are valid
/// - **Open intervals `(a,b)`**: Neither endpoint is valid  
/// - **Semi-open intervals**: Only one endpoint is valid
/// - **Unbounded intervals**: Appropriate infinite bounds are handled
///
/// ### Integration with ftl-base Ecosystem
/// - **Chebyshev Transformations**: Primary use case for parametric ↔ physical domain mappings
/// - **Domain Transformations**: Ensures transformed points remain within target domains
/// - **Numerical Method Input**: Provides validated input for interpolation and integration
/// - **Grid Construction**: Can serve as input to grid builders when domain validation is critical
///
/// ## Type Parameters
///
/// - **`IntervalType`**: The interval type implementing [`IntervalTrait`]
///   - Determines the mathematical properties of the domain (bounds, openness/closedness) **at compile time**
///   - Must be consistent with the coordinate type through `IntervalType::RealType` **at compile time**
///   - Examples: [`IntervalClosed<RealNative64StrictFinite>`](crate::intervals::IntervalClosed), [`IntervalOpen<f64>`](crate::intervals::IntervalOpen),
///     [`IntervalLowerClosedUpperUnbounded<f64>`](crate::intervals::IntervalLowerClosedUpperUnbounded)
///
/// - **`CoordsStorage`**: Generic storage type for coordinates (defaults to owned [`Coords1D`])
///   - **Constraint**: Must implement `std::borrow::Borrow<Coords1D<IntervalType::RealType>>`
///   - **Enables flexible memory management** for coordinates without changing the API
///   - **Zero-cost abstraction** with compile-time dispatch to the appropriate storage method
///   - **Supported storage types**:
///     - **Owned**: `Coords1D<T>` (default) - Direct ownership, movable coordinates
///     - **Borrowed**: `&Coords1D<T>` - Zero-cost borrowed reference with lifetime
///     - **Thread-safe shared**: `Arc<Coords1D<T>>` - Concurrent access across threads
///     - **Single-threaded shared**: `Rc<Coords1D<T>>` - Reference-counted sharing
///     - **Unique heap**: `Box<Coords1D<T>>` - Heap allocation for large coordinate sets
///     - **Clone-on-write**: `Cow<Coords1D<T>>` - Efficient lazy cloning semantics
///   - **Use cases**: Memory efficiency for large datasets, thread-safe coordinate sharing,
///     zero-copy operations with borrowed data, flexible ownership patterns
///
/// ## Validation Strategy
///
/// ### Construction-Time Validation (Runtime)
/// The [`try_new`](Self::try_new()) constructor performs **runtime validation**:
/// 1. **Efficiency**: Only checks the minimum and maximum coordinates (O(1) validation)
/// 2. **Completeness**: Since coordinates are sorted, min/max check ensures all points are valid
/// 3. **Fail-Fast**: Invalid inputs are rejected immediately with detailed error information
///
/// ### Post-Construction Guarantees
/// Once constructed successfully:
/// - **Zero Runtime Overhead**: No additional validation needed during use
/// - **Immutable Invariants**: Neither coordinates nor domain can violate constraints
/// - **Type Safety**: The interval type provides compile-time information about domain properties
///
/// ## Performance Characteristics
///
/// | Operation | Time Complexity | Space Complexity | Notes |
/// |-----------|----------------|------------------|--------|
/// | **Construction** | O(1) | O(n) | **Runtime validation** of min/max only |
/// | **Coordinate access** | O(1) | O(1) | Direct delegation to [`Coords1D`] |
/// | **Domain access** | O(1) | O(1) | Direct field access |
/// | **Point validation** | O(1) | O(1) | **Not needed post-construction** |
/// | **Conversion to [`Coords1D`]** | O(1) | O(1) | Zero-cost extraction |
///
/// The **runtime validation** cost is minimal since only the minimum and maximum coordinates need to be checked
/// against the domain bounds, regardless of the total number of points.
///
/// ## Construction and Usage Examples
///
/// ### Basic Construction with Validation
/// ```rust
/// use grid1d::{
///     Coords1D, Coords1DInDomain, HasCoords1D, HasDomain1D,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// // Create coordinates
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])
/// ).unwrap();
///
/// // Define domain
/// let domain = IntervalClosed::new(0.0, 1.0);
///
/// // Create validated points in interval
/// let points = Coords1DInDomain::try_new(coords, domain).unwrap();
///
/// // Access validated coordinates and domain
/// assert_eq!(points.coords().deref(), &[0.1, 0.5, 0.9]);
/// assert_eq!(points.domain().lower_bound_value(), &0.0);
/// assert_eq!(points.domain().upper_bound_value(), &1.0);
/// ```
///
/// ### Domain Type Variations
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, intervals::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])
/// ).unwrap();
///
/// // Closed interval [0,1] - includes both endpoints
/// let closed_points = Coords1DInDomain::try_new(
///     coords.clone(),
///     IntervalClosed::new(0.0, 1.0)
/// ).unwrap();
///
/// // Open interval (0,1) - excludes both endpoints  
/// let open_points = Coords1DInDomain::try_new(
///     coords.clone(),
///     IntervalOpen::new(0.0, 1.0)
/// ).unwrap();
///
/// // Semi-open interval [0,1) - includes left, excludes right
/// let semi_open_points = Coords1DInDomain::try_new(
///     coords,
///     IntervalLowerClosedUpperOpen::new(0.0, 1.0)
/// ).unwrap();
/// ```
///
/// ### Working with Different Scalar Types
/// ```rust
/// use grid1d::{
///     Coords1D, HasCoords1D, Coords1DInDomain,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Create validated coordinates with strict finite type
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
///     Real::try_new(0.1).unwrap(),
///     Real::try_new(0.5).unwrap(),
///     Real::try_new(0.9).unwrap(),
/// ])).unwrap();
///
/// let domain = IntervalClosed::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
///
/// let points = Coords1DInDomain::try_new(coords, domain).unwrap();
/// assert_eq!(points.coords().num_points().as_ref(), &3);
/// ```
///
/// ## Error Handling and Validation
///
/// ### Domain Violation Detection
/// ```rust
/// use grid1d::{
///     coords::{Coords1D, Coords1DInDomain, ErrorsCoords1DInDomain,},
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // Coordinates that exceed the domain bounds
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 0.5, 1.5]) // 1.5 exceeds upper bound
/// ).unwrap();
/// let domain = IntervalClosed::new(0.0, 1.0);
///
/// match Coords1DInDomain::try_new(coords, domain) {
///     Err(ErrorsCoords1DInDomain::UpperBoundViolation { coord_max, domain, .. }) => {
///         println!("Coordinate {} exceeds domain upper bound {}",
///                  coord_max, domain.upper_bound_value());
///         assert_eq!(coord_max, 1.5);
///     }
///     _ => panic!("Expected upper bound violation"),
/// }
/// ```
///
/// ### Lower Bound Violation
/// ```rust
/// use grid1d::{
///     coords::{Coords1D, Coords1DInDomain, ErrorsCoords1DInDomain,},
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // Coordinates that fall below the domain bounds
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![-0.5, 0.5, 1.0]) // -0.5 below lower bound
/// ).unwrap();
/// let domain = IntervalClosed::new(0.0, 1.0);
///
/// match Coords1DInDomain::try_new(coords, domain) {
///     Err(ErrorsCoords1DInDomain::LowerBoundViolation { coord_min, domain, .. }) => {
///         println!("Coordinate {} below domain lower bound {}",
///                  coord_min, domain.lower_bound_value());
///         assert_eq!(coord_min, -0.5);
///     }
///     _ => panic!("Expected lower bound violation"),
/// }
/// ```
///
/// ### Robust Construction Pattern
/// ```rust
/// use grid1d::{
///     coords::{Coords1D, Coords1DInDomain, ErrorsCoords1DInDomain,},
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// fn safe_points_construction<I: IntervalTrait>(
///     values: Vec<I::RealType>,
///     domain: I,
/// ) -> Result<Coords1DInDomain<I>, String>
/// where
///     I::RealType: std::fmt::Display,
/// {
///     // First create valid coordinates
///     let coords = Coords1D::try_from(SortedSet::from_unsorted(values))
///         .map_err(|e| format!("Invalid coordinates: {}", e))?;
///     
///     // Then validate against domain
///     Coords1DInDomain::try_new(coords, domain)
///         .map_err(|e| match e {
///             ErrorsCoords1DInDomain::LowerBoundViolation { coord_min, .. } => {
///                 format!("Coordinate {} is below the domain lower bound", coord_min)
///             }
///             ErrorsCoords1DInDomain::UpperBoundViolation { coord_max, .. } => {
///                 format!("Coordinate {} exceeds the domain upper bound", coord_max)
///             }
///         })
/// }
///
/// // Usage
/// let domain = IntervalClosed::new(0.0, 2.0);
/// let result = safe_points_construction(vec![0.5, 1.0, 1.5], domain.clone());
/// assert!(result.is_ok());
///
/// let bad_result = safe_points_construction(vec![0.5, 1.0, 2.5], domain);
/// assert!(bad_result.is_err());
/// ```
///
/// ### Custom Domain Transformations
/// ```rust
/// use grid1d::{
///     Coords1D, Coords1DInDomain, HasCoords1D, HasDomain1D,
///     intervals::{IntervalClosed, IntervalFinitePositiveLengthTrait, GetLowerBoundValue},
/// };
/// use sorted_vec::partial::SortedSet;
/// use std::ops::Deref;
///
/// fn transform_points_with_validation(
///     source_points: &Coords1DInDomain<IntervalClosed<f64>>,
///     target_domain: IntervalClosed<f64>,
/// ) -> Result<Coords1DInDomain<IntervalClosed<f64>>, String> {
///     // Transform coordinates (example: scale and shift)
///     let source_domain = source_points.domain();
///     let scale = target_domain.length().as_ref() / source_domain.length().as_ref();
///     let shift = target_domain.lower_bound_value() -
///                source_domain.lower_bound_value() * scale;
///     
///     let transformed_coords: Vec<f64> = source_points
///         .coords()
///         .iter()
///         .map(|&x| x * scale + shift)
///         .collect();
///     
///     // Create new validated coordinates
///     let coords = Coords1D::try_from(SortedSet::from_unsorted(transformed_coords))
///         .map_err(|e| format!("Coordinate creation failed: {}", e))?;
///     
///     Coords1DInDomain::try_new(coords, target_domain)
///         .map_err(|e| format!("Domain validation failed: {:?}", e))
/// }
///
/// // Usage
/// let source_coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])
/// ).unwrap();
/// let source_points = Coords1DInDomain::try_new(
///     source_coords,
///     IntervalClosed::new(0.0, 1.0)
/// ).unwrap();
///
/// let target_domain = IntervalClosed::new(-5.0, 5.0);
/// let transformed = transform_points_with_validation(&source_points, target_domain).unwrap();
/// assert_eq!(transformed.coords().deref(), &[-5.0, 0.0, 5.0]);
/// ```
///
/// ## Flexible Storage Types and Memory Management
///
/// The `CoordsStorage` generic parameter enables different storage strategies for coordinates,
/// providing memory efficiency and ownership flexibility without changing the API.
///
/// ### Owned Coordinates (Default)
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, intervals::*};
/// use sorted_vec::partial::SortedSet;
///
/// // Default behavior - owned coordinates
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let points = Coords1DInDomain::try_new(coords, IntervalClosed::new(0.0, 1.0)).unwrap();
/// // `points` owns the coordinate data
/// ```
///
/// ### Borrowed Coordinates (Zero-Cost References)
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, Coords1DInDomainBorrowed, intervals::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let domain = IntervalClosed::new(0.0, 1.0);
///
/// // Zero-cost borrowed reference - no copying
/// let borrowed_points: Coords1DInDomainBorrowed<_> =
///     Coords1DInDomain::try_new(&coords, domain.clone()).unwrap();
///
/// // Or convert from owned to borrowed view
/// let owned_points = Coords1DInDomain::try_new(coords, domain).unwrap();
/// let borrowed_view = owned_points.as_borrowed();
/// ```
///
/// ### Thread-Safe Shared Ownership (Arc)
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, Coords1DInDomainArc, intervals::*};
/// use sorted_vec::partial::SortedSet;
/// use std::sync::Arc;
///
/// // Create large coordinate set (expensive to copy)
/// let large_coords = Coords1D::try_from(SortedSet::from_unsorted(
///     (0..1000).map(|i| i as f64 / 999.0).collect()
/// )).unwrap();
///
/// // Wrap in Arc for thread-safe sharing
/// let arc_coords = Arc::new(large_coords);
/// let points: Coords1DInDomainArc<_> =
///     Coords1DInDomain::try_new(arc_coords, IntervalClosed::new(0.0, 1.0)).unwrap();
///
/// // Can be shared across threads without copying coordinate data
/// // let points_clone = points.clone(); // Cheap Arc clone
/// ```
///
/// ### Single-Threaded Shared Ownership (Rc)
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, Coords1DInDomainRc, intervals::*};
/// use sorted_vec::partial::SortedSet;
/// use std::rc::Rc;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let rc_coords = Rc::new(coords);
///
/// let points: Coords1DInDomainRc<_> =
///     Coords1DInDomain::try_new(rc_coords, IntervalClosed::new(0.0, 1.0)).unwrap();
///
/// // Efficient single-threaded sharing with reference counting
/// let points_clone = points.clone(); // Cheap Rc clone
/// ```
///
/// ### Unique Heap Allocation (Box)
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, Coords1DInDomainBox, intervals::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let boxed_coords = Box::new(coords);
///
/// let points: Coords1DInDomainBox<_> =
///     Coords1DInDomain::try_new(boxed_coords, IntervalClosed::new(0.0, 1.0)).unwrap();
/// // Useful for large coordinate sets to keep them on the heap
/// ```
///
/// ### Clone-on-Write Semantics (Cow)
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, Coords1DInDomainCow, intervals::*};
/// use sorted_vec::partial::SortedSet;
/// use std::borrow::Cow;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let cow_coords = Cow::Borrowed(&coords);
///
/// let points: Coords1DInDomainCow<_> =
///     Coords1DInDomain::try_new(cow_coords, IntervalClosed::new(0.0, 1.0)).unwrap();
/// // Coordinates are borrowed until modification is needed
/// ```
///
/// ### Generic Functions Over Storage Types
/// ```rust
/// use grid1d::{Coords1DInDomain, HasCoords1D, intervals::IntervalTrait};
///
/// // Function works with any storage type
/// fn analyze_coordinates<I, S>(points: &Coords1DInDomain<I, S>) -> usize
/// where
///     I: IntervalTrait,
///     S: std::borrow::Borrow<grid1d::Coords1D<I::RealType>>,
/// {
///     points.coords().len()
/// }
///
/// // Works with owned, borrowed, Arc, Rc, Box, Cow...
/// # use grid1d::{Coords1D, intervals::*};
/// # use sorted_vec::partial::SortedSet;
/// # let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// # let domain = IntervalClosed::new(0.0, 1.0);
/// let owned = Coords1DInDomain::try_new(coords.clone(), domain.clone()).unwrap();
/// let borrowed = Coords1DInDomain::try_new(&coords, domain).unwrap();
///
/// assert_eq!(analyze_coordinates(&owned), 3);
/// assert_eq!(analyze_coordinates(&borrowed), 3);
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Generic Programming Over Domain Types
/// ```rust
/// use grid1d::{
///     Coords1DInDomain, HasCoords1D, HasDomain1D,
///     intervals::*,
/// };
/// use num_valid::RealScalar;
///
/// fn analyze_points<I: IntervalTrait, S>(points: &Coords1DInDomain<I, S>) -> String
/// where
///     I::RealType: std::fmt::Display + Clone,
///     S: std::borrow::Borrow<Coords1D<I::RealType>>,
/// {
///     let coords = points.coords();
///     let domain = points.domain();
///     let span = coords.last().clone() - coords.first();
///     
///     format!(
///         "Points: {} coordinates spanning {} within domain bounds [{}, {}]",
///         coords.len(),
///         span,
///         // Note: These methods would need to be available on IntervalTrait
///         "lower_bound", "upper_bound"
///     )
/// }
///
/// // Works with any interval type
/// # use sorted_vec::partial::SortedSet;
/// # use grid1d::Coords1D;
/// # let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let closed_points = Coords1DInDomain::try_new(coords.clone(), IntervalClosed::new(0.0, 1.0)).unwrap();
/// let open_points = Coords1DInDomain::try_new(coords, IntervalOpen::new(0.0, 1.0)).unwrap();
/// ```
///
/// ### Conversion and Interoperability
/// ```rust
/// use grid1d::{
///     Coords1D, Coords1DInDomain, HasCoords1D,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// let points = Coords1DInDomain::try_new(
///     Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap(),
///     IntervalClosed::new(0.0, 1.0)
/// ).unwrap();
///
/// // Extract coordinates for use with other algorithms
/// let coords: &Coords1D<f64> = points.coords();
///
/// // Convert to vector for external libraries
/// let coord_vec: Vec<f64> = coords.to_vec();
///
/// for &coord in coords.iter() {
///     // High-performance processing...
/// }
/// ```
///
/// ### Batch Validation and Processing
/// ```rust
/// use grid1d::{
///     coords::{Coords1D, Coords1DInDomain, ErrorsCoords1DInDomain,},
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// fn validate_multiple_point_sets(
///     coord_sets: Vec<Vec<f64>>,
///     domain: IntervalClosed<f64>,
/// ) -> (Vec<Coords1DInDomain<IntervalClosed<f64>>>, Vec<String>) {
///     let mut valid_sets = Vec::new();
///     let mut errors = Vec::new();
///     
///     for (i, coords_vec) in coord_sets.into_iter().enumerate() {
///         match Coords1D::try_from(SortedSet::from_unsorted(coords_vec)) {
///             Ok(coords) => {
///                 match Coords1DInDomain::try_new(coords, domain.clone()) {
///                     Ok(points) => valid_sets.push(points),
///                     Err(e) => errors.push(format!("Set {}: {:?}", i, e)),
///                 }
///             }
///             Err(e) => errors.push(format!("Set {} coordinate error: {}", i, e)),
///         }
///     }
///     
///     (valid_sets, errors)
/// }
///
/// // Usage
/// let coord_sets = vec![
///     vec![0.1, 0.5, 0.9],  // Valid
///     vec![0.0, 0.5, 1.0],  // Valid  
///     vec![-0.1, 0.5, 0.9], // Invalid (lower bound violation)
///     vec![0.1, 0.5, 1.1],  // Invalid (upper bound violation)
/// ];
///
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let (valid, errors) = validate_multiple_point_sets(coord_sets, domain);
///
/// assert_eq!(valid.len(), 2);
/// assert_eq!(errors.len(), 2);
/// ```
///
/// ## Memory Layout and Performance
///
/// ### Memory Efficiency
/// ```rust
/// use grid1d::{
///     Coords1D, Coords1DInDomain,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
/// use std::mem;
///
/// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.1, 0.5, 0.9])).unwrap();
/// let domain = IntervalClosed::new(0.0, 1.0);
/// let points = Coords1DInDomain::try_new(coords.clone(), domain.clone()).unwrap();
///
/// // Memory layout is essentially Coords1D + Domain
/// let coords_size = mem::size_of_val(&coords);
/// let domain_size = mem::size_of_val(&domain);
/// let points_size = mem::size_of_val(&points);
///
/// // Coords1DInDomain should have minimal overhead beyond its components
/// assert!(points_size <= coords_size + domain_size + 16); // Allow for some padding
/// ```
///
/// ### Performance-Critical Access
/// ```rust
/// use grid1d::{
///     Coords1DInDomain, HasCoords1D,
///     intervals::*,
/// };
/// # use grid1d::Coords1D;
/// # use sorted_vec::partial::SortedSet;
///
/// # let points = Coords1DInDomain::try_new(
/// #     Coords1D::try_from(SortedSet::from_unsorted((0..1000).map(|i| i as f64 / 1000.0).collect())).unwrap(),
/// #     IntervalClosed::new(0.0, 1.0)
/// # ).unwrap();
///
/// // For performance-critical numerical computations
/// let slice = points.coords();
///
/// // Cache-friendly sequential access
/// let mut sum = 0.0;
/// for &coord in slice.iter() {
///     sum += coord.powi(2); // Example: sum of squares
/// }
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`Coords1DInDomain`]
/// - **Domain transformations**: When mapping between parametric and physical spaces
/// - **Chebyshev methods**: For validated interpolation and quadrature points
/// - **Boundary condition application**: When coordinates must respect domain constraints
/// - **Validated numerical input**: For algorithms that require domain-consistent data
/// - **Multi-domain calculations**: When tracking coordinates across domain boundaries
///
/// ### When to Use [`Coords1D`] Instead
/// - **Pure coordinate storage**: When domain validation is unnecessary
/// - **Grid construction**: When coordinates define the domain rather than being constrained by it
/// - **General numerical arrays**: For coordinates without specific domain requirements
/// - **Performance-critical paths**: When validation overhead is unacceptable
/// - **Coordinate manipulation**: When you need to modify coordinates frequently
///
/// ### Design Patterns
/// ```rust
/// use grid1d::{
///     Coords1D, HasCoords1D, Coords1DInDomain,
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // Pattern 1: Staged validation
/// fn create_validated_interpolation_points(
///     raw_coords: Vec<f64>,
///     domain: IntervalClosed<f64>,
/// ) -> Result<Coords1DInDomain<IntervalClosed<f64>>, String> {
///     // Stage 1: Basic coordinate validation
///     let coords = Coords1D::try_from(SortedSet::from_unsorted(raw_coords))
///         .map_err(|e| format!("Invalid coordinates: {}", e))?;
///     
///     // Stage 2: Domain validation
///     Coords1DInDomain::try_new(coords, domain)
///         .map_err(|e| format!("Domain validation failed: {:?}", e))
/// }
///
/// // Pattern 2: Safe extraction for processing
/// fn process_coordinates_safely(
///     points: &Coords1DInDomain<IntervalClosed<f64>>
/// ) -> Vec<f64> {
///     // Extract coordinates knowing they're domain-valid
///     points.coords()
///         .iter()
///         .map(|&x| x.sqrt()) // Safe because x >= 0 is guaranteed by domain
///         .collect()
/// }
///
/// // Pattern 3: Domain-aware algorithms
/// fn interpolate_with_validation<I: IntervalTrait>(
///     points: &Coords1DInDomain<I>,
///     values: &[I::RealType],
/// ) -> Result<(), String>
/// where
///     I::RealType: Clone + PartialOrd,
/// {
///     if points.coords().len() != values.len() {
///         return Err("Coordinate and value count mismatch".to_string());
///     }
///     
///     // The algorithm can trust that all coordinates are within the domain
///     // No additional bounds checking needed
///     Ok(())
/// }
/// ```
///
/// ## Mathematical Properties and Guarantees
///
/// The [`Coords1DInDomain`] struct maintains these mathematical properties:
///
/// ### Domain Containment Invariant (Runtime-Validated)
/// ```text
/// ∀i ∈ [0, n): domain.contains_point(&coords[i]) = true
/// ```
/// This is the fundamental guarantee that distinguishes [`Coords1DInDomain`] from [`Coords1D`].
/// **This is validated at runtime during construction, not at compile time.**
///
/// ### Inherited Properties from [`Coords1D`] (Construction-Time Guaranteed)
/// - **Strict ordering**: `coords[i] < coords[i+1]` for all valid indices `i`
/// - **Uniqueness**: `coords[i] ≠ coords[j]` for all `i ≠ j`  
/// - **Non-emptiness**: `coords.len() ≥ 1` always holds
/// - **Finite values**: All coordinates are finite real numbers (when using validated scalar types)
///
/// ### Type Safety Properties (Compile-Time)
/// - **Domain-Type Consistency**: The interval type parameter ensures compile-time domain type checking
/// - **Coordinate-Domain Alignment**: Scalar types must match between coordinates and domain **at compile time**
/// - **Immutable Invariants**: Once constructed, neither coordinates nor domain can violate constraints
///
/// ### Validation Efficiency (Runtime)
/// - **O(1) Construction**: Only min/max coordinates are checked, not every point
/// - **Zero Runtime Overhead**: After construction, no additional validation is needed
/// - **Fail-Fast**: Invalid inputs are rejected immediately at construction time **with runtime checks**
///
/// These properties make [`Coords1DInDomain`] ideal for numerical algorithms that depend on
/// domain validity, providing both mathematical correctness and computational efficiency **through
/// a combination of compile-time type safety and runtime validation**.
///
/// ## Error Handling and Validation
///
/// ### Construction Error Handling (Runtime Validation)
/// ```rust
/// use grid1d::{
///     coords::{Coords1D, Coords1DInDomain, ErrorsCoords1DInDomain},
///     intervals::*,
/// };
/// use sorted_vec::partial::SortedSet;
///
/// // Coordinates that exceed the domain bounds
/// let coords = Coords1D::try_from(
///     SortedSet::from_unsorted(vec![0.0, 0.5, 1.5]) // 1.5 exceeds upper bound
/// ).unwrap();
/// let domain = IntervalClosed::new(0.0, 1.0);
///
/// // This will fail at RUNTIME, not compile time
/// match Coords1DInDomain::try_new(coords, domain) {
///     Err(ErrorsCoords1DInDomain::UpperBoundViolation { coord_max, domain, .. }) => {
///         println!("RUNTIME validation caught: Coordinate {} exceeds domain upper bound {}",
///                  coord_max, domain.upper_bound_value());
///         assert_eq!(coord_max, 1.5);
///     }
///     _ => panic!("Expected upper bound violation"),
/// }
/// ```
#[derive(Debug, Clone, PartialEq)]
pub struct Coords1DInDomain<
    IntervalType: IntervalTrait,
    CoordsStorage = Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>,
> where
    CoordsStorage: std::borrow::Borrow<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
{
    /// The coordinates of the points.
    ///
    /// # Note
    /// These values are sorted in ascending order.
    coords: CoordsStorage,

    /// The interval that defines the domain of the points.
    domain: IntervalType,
}

/// Type alias for owned coordinates in domain.
pub type Coords1DInDomainOwned<IntervalType> =
    Coords1DInDomain<IntervalType, Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>;

/// Type alias for borrowed coordinates in domain.
pub type Coords1DInDomainBorrowed<'a, IntervalType> =
    Coords1DInDomain<IntervalType, &'a Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>;

/// Type alias for Arc-wrapped coordinates in domain (thread-safe, shared ownership).
pub type Coords1DInDomainArc<IntervalType> = Coords1DInDomain<
    IntervalType,
    std::sync::Arc<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
>;

/// Type alias for Rc-wrapped coordinates in domain (single-threaded, shared ownership).
pub type Coords1DInDomainRc<IntervalType> = Coords1DInDomain<
    IntervalType,
    std::rc::Rc<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
>;

/// Type alias for Box-wrapped coordinates in domain (unique ownership on heap).
pub type Coords1DInDomainBox<IntervalType> = Coords1DInDomain<
    IntervalType,
    Box<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
>;

/// Type alias for Cow-wrapped coordinates in domain (clone-on-write).
pub type Coords1DInDomainCow<'a, IntervalType> = Coords1DInDomain<
    IntervalType,
    std::borrow::Cow<'a, Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
>;

impl<IntervalType: IntervalTrait, CoordsStorage> Coords1DInDomain<IntervalType, CoordsStorage>
where
    CoordsStorage: std::borrow::Borrow<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
{
    /// Attempts to create a new [`Coords1DInDomain`] instance.
    ///
    /// # Parameters
    ///
    /// - `coords`: A [`Coords1D`] object or reference containing the sorted coordinates of the points.
    /// - `domain`: An interval object defining the domain in which the points must lie.
    ///
    /// # Returns
    ///
    /// - `Ok(Self)`: If all points lie within the specified domain.
    /// - [`Err<ErrorsCoords1DInDomain>`]: If any point lies outside the domain.
    ///
    /// # Errors
    ///
    /// - Returns [`ErrorsCoords1DInDomain::LowerBoundViolation`] if the smallest coordinate is less than the lower bound of the domain.
    /// - Returns [`ErrorsCoords1DInDomain::UpperBoundViolation`] if the largest coordinate is greater than the upper bound of the domain.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Coords1D, Coords1DInDomain, HasCoords1D,
    ///     intervals::*,
    /// };
    /// use sorted_vec::partial::SortedSet;
    /// use std::ops::Deref;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
    /// let domain = IntervalClosed::new(0.0, 1.0);
    ///
    /// // Owned coordinates
    /// let points_owned = Coords1DInDomain::try_new(coords.clone(), domain.clone()).unwrap();
    ///
    /// // Borrowed coordinates  
    /// let points_borrowed = Coords1DInDomain::try_new(&coords, domain).unwrap();
    ///
    /// assert_eq!(points_owned.coords().as_ref().as_slice(), &[0.0, 0.5, 1.0]);
    /// assert_eq!(points_borrowed.coords().as_ref().as_slice(), &[0.0, 0.5, 1.0]);
    ///
    /// // Conversion between owned and borrowed
    /// let borrowed_view = points_owned.as_borrowed();
    /// let owned_copy = points_borrowed.to_owned_coords();
    ///
    /// // Type aliases for convenience
    /// let _owned: grid1d::Coords1DInDomainOwned<IntervalClosed<f64>> = points_owned;
    /// let _borrowed: grid1d::Coords1DInDomainBorrowed<IntervalClosed<f64>> = points_borrowed;
    /// ```
    pub fn try_new(
        coords: CoordsStorage,
        domain: IntervalType,
    ) -> Result<Self, ErrorsCoords1DInDomain<IntervalType>> {
        let coords_ref = coords.borrow();
        let coords_min = coords_ref.first();
        if !domain.contains_point(coords_min) {
            return Err(ErrorsCoords1DInDomain::LowerBoundViolation {
                coord_min: coords_min.clone(),
                domain,
                backtrace: capture_backtrace(),
            });
        }

        let coords_max = coords_ref.last();
        if !domain.contains_point(coords_max) {
            return Err(ErrorsCoords1DInDomain::UpperBoundViolation {
                coord_max: coords_max.clone(),
                domain,
                backtrace: capture_backtrace(),
            });
        };

        Ok(Coords1DInDomain { coords, domain })
    }
}

// Specialized implementations for owned coordinates
impl<IntervalType: IntervalTrait>
    Coords1DInDomain<IntervalType, Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>
{
    /// Creates a borrowed view of this owned `Coords1DInDomain`.
    ///
    /// This is useful when you have owned coordinates but need to pass them to a function
    /// that expects borrowed coordinates without transferring ownership.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{Coords1D, Coords1DInDomain, HasCoords1D, intervals::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
    /// let domain = IntervalClosed::new(0.0, 1.0);
    /// let owned_points = Coords1DInDomain::try_new(coords, domain).unwrap();
    ///
    /// let borrowed_points = owned_points.as_borrowed();
    /// // Both point to the same underlying data
    /// assert_eq!(owned_points.coords(), borrowed_points.coords());
    /// ```
    pub fn as_borrowed(
        &self,
    ) -> Coords1DInDomain<IntervalType, &Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>
    where
        IntervalType: Clone,
    {
        Coords1DInDomain {
            coords: &self.coords,
            domain: self.domain.clone(),
        }
    }
}

// Specialized implementations for borrowed coordinates
impl<IntervalType: IntervalTrait>
    Coords1DInDomain<IntervalType, &Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>
{
    /// Creates an owned copy of this borrowed `Coords1DInDomain`.
    ///
    /// This clones both the coordinates and the domain to create a fully owned instance.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{Coords1D, Coords1DInDomain, HasCoords1D, intervals::*};
    /// use sorted_vec::partial::SortedSet;
    ///
    /// let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
    /// let domain = IntervalClosed::new(0.0, 1.0);
    /// let borrowed_points = Coords1DInDomain::try_new(&coords, domain).unwrap();
    ///
    /// let owned_points = borrowed_points.to_owned_coords();
    /// // Now we have an owned copy
    /// assert_eq!(borrowed_points.coords(), owned_points.coords());
    /// ```
    pub fn to_owned_coords(
        &self,
    ) -> Coords1DInDomain<IntervalType, Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>
    where
        IntervalType: Clone,
        <IntervalType as IntervalBoundsRuntime>::RealType: Clone,
    {
        Coords1DInDomain {
            coords: self.coords.clone(),
            domain: self.domain.clone(),
        }
    }
}

impl<IntervalType: IntervalTrait, CoordsStorage> HasDomain1D
    for Coords1DInDomain<IntervalType, CoordsStorage>
where
    CoordsStorage: std::borrow::Borrow<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
{
    type Domain1D = IntervalType;

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

impl<IntervalType: IntervalTrait, CoordsStorage> HasCoords1D
    for Coords1DInDomain<IntervalType, CoordsStorage>
where
    CoordsStorage: std::borrow::Borrow<Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType>>,
{
    type Point1DType = <IntervalType as IntervalBoundsRuntime>::RealType;

    fn coords(&self) -> &Coords1D<<IntervalType as IntervalBoundsRuntime>::RealType> {
        self.coords.borrow()
    }
}
//-------------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::intervals::*;
    use sorted_vec::partial::SortedSet;
    use std::ops::Deref;

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

        type Real = RealNative64StrictFinite;

        mod coords1d {
            use super::*;

            #[test]
            fn test_coords1d_creation() {
                let sorted_set = SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]);
                let coords = Coords1D::try_from(sorted_set).unwrap();

                assert_eq!(
                    coords.deref(),
                    &[
                        Real::try_new(0.0).unwrap(),
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap()
                    ]
                );
            }

            #[test]
            fn test_coords1d_empty_creation() {
                let empty_set = SortedSet::<Real>::new();
                let coords1d = Coords1D::try_from(empty_set);
                matches!(coords1d, Err(ErrorsCoords1D::EmptyObject { .. }));
            }

            #[test]
            fn test_new_uniform() {
                let domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(10.0).unwrap());
                let coords = Coords1D::new_uniform(&domain, &NumIntervals::try_new(5).unwrap());
                assert_eq!(
                    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!(coords.len(), 6);
            }

            #[test]
            fn test_first_last_span() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(10.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(coords.first(), &Real::try_new(1.0).unwrap());
                assert_eq!(coords.last(), &Real::try_new(10.0).unwrap());
                assert_eq!(coords.span(), Real::try_new(9.0).unwrap());
            }

            #[test]
            fn test_contains_and_find_index() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(7.0).unwrap(),
                ]))
                .unwrap();
                // Test contains
                assert!(coords.contains(&Real::try_new(3.0).unwrap()));
                assert!(coords.contains(&Real::try_new(7.0).unwrap()));
                assert!(!coords.contains(&Real::try_new(4.0).unwrap()));
                assert!(!coords.contains(&Real::try_new(0.0).unwrap()));

                // Test find_index
                assert_eq!(coords.find_index(&Real::try_new(1.0).unwrap()), Some(0));
                assert_eq!(coords.find_index(&Real::try_new(5.0).unwrap()), Some(2));
                assert_eq!(coords.find_index(&Real::try_new(6.0).unwrap()), None);
            }

            #[test]
            fn test_find_insertion_index() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(20.0).unwrap(),
                    Real::try_new(30.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(coords.find_insertion_index(&Real::try_new(5.0).unwrap()), 0); // Before start
                assert_eq!(
                    coords.find_insertion_index(&Real::try_new(10.0).unwrap()),
                    0
                ); // Exact match at start
                assert_eq!(
                    coords.find_insertion_index(&Real::try_new(15.0).unwrap()),
                    1
                ); // Middle
                assert_eq!(
                    coords.find_insertion_index(&Real::try_new(20.0).unwrap()),
                    1
                ); // Exact match in middle
                assert_eq!(
                    coords.find_insertion_index(&Real::try_new(35.0).unwrap()),
                    3
                ); // After end
            }

            #[test]
            fn test_find_floor_index() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(20.0).unwrap(),
                    Real::try_new(30.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(coords.find_floor_index(&Real::try_new(5.0).unwrap()), None);
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(10.0).unwrap()),
                    Some(0)
                );
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(15.0).unwrap()),
                    Some(0)
                );
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(29.9).unwrap()),
                    Some(1)
                );
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(30.0).unwrap()),
                    Some(2)
                );
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(35.0).unwrap()),
                    Some(2)
                );
            }

            #[test]
            fn test_find_ceil_index() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(20.0).unwrap(),
                    Real::try_new(30.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(
                    coords.find_ceil_index(&Real::try_new(5.0).unwrap()),
                    Some(0)
                );
                assert_eq!(
                    coords.find_ceil_index(&Real::try_new(10.0).unwrap()),
                    Some(0)
                );
                assert_eq!(
                    coords.find_ceil_index(&Real::try_new(10.1).unwrap()),
                    Some(1)
                );
                assert_eq!(
                    coords.find_ceil_index(&Real::try_new(30.0).unwrap()),
                    Some(2)
                );
                assert_eq!(coords.find_ceil_index(&Real::try_new(35.0).unwrap()), None);
            }

            #[test]
            fn test_find_range_indices_and_range() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(6.0).unwrap(),
                ]))
                .unwrap();
                let range_interval = IntervalClosed::try_new(
                    Real::try_new(2.5).unwrap(),
                    Real::try_new(5.5).unwrap(),
                )
                .unwrap();

                // Test find_range_indices
                let indices = coords.find_range_indices(&range_interval);
                assert_eq!(indices, 2..5);

                // Test range
                let values: Vec<Real> = coords.range(&range_interval).cloned().collect();
                assert_eq!(
                    values,
                    vec![
                        Real::try_new(3.0).unwrap(),
                        Real::try_new(4.0).unwrap(),
                        Real::try_new(5.0).unwrap()
                    ]
                );

                // Test empty range
                let empty_interval = IntervalClosed::try_new(
                    Real::try_new(7.0).unwrap(),
                    Real::try_new(8.0).unwrap(),
                )
                .unwrap();
                assert_eq!(coords.find_range_indices(&empty_interval), 6..6);
                assert_eq!(coords.range(&empty_interval).count(), 0);
            }

            #[test]
            fn test_find_closest() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(20.0).unwrap(),
                    Real::try_new(30.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(
                    coords.find_closest(&Real::try_new(5.0).unwrap()),
                    (&Real::try_new(10.0).unwrap(), 0)
                ); // Below range
                assert_eq!(
                    coords.find_closest(&Real::try_new(14.0).unwrap()),
                    (&Real::try_new(10.0).unwrap(), 0)
                ); // Closer to left
                assert_eq!(
                    coords.find_closest(&Real::try_new(16.0).unwrap()),
                    (&Real::try_new(20.0).unwrap(), 1)
                ); // Closer to right
                assert_eq!(
                    coords.find_closest(&Real::try_new(20.0).unwrap()),
                    (&Real::try_new(20.0).unwrap(), 1)
                ); // Exact match
                assert_eq!(
                    coords.find_closest(&Real::try_new(35.0).unwrap()),
                    (&Real::try_new(30.0).unwrap(), 2)
                ); // Above range
            }

            #[test]
            fn test_set_operations() {
                let coords1 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                ]))
                .unwrap();
                let coords2 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(6.0).unwrap(),
                ]))
                .unwrap();
                let coords3 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]))
                .unwrap();

                // Test union
                let union = coords1.union(&coords2);
                assert_eq!(
                    union.deref(),
                    &[
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(3.0).unwrap(),
                        Real::try_new(4.0).unwrap(),
                        Real::try_new(5.0).unwrap(),
                        Real::try_new(6.0).unwrap()
                    ]
                );

                // Test intersection
                let intersection = coords1.intersection(&coords2).unwrap();
                assert_eq!(
                    intersection.deref(),
                    &[Real::try_new(3.0).unwrap(), Real::try_new(4.0).unwrap()]
                );
                assert!(coords1.intersection(&coords3).is_some());
                assert!(coords2.intersection(&coords3).is_none());

                // Test is_subset
                assert!(coords3.is_subset(&coords1));
                assert!(!coords1.is_subset(&coords3));
                assert!(!coords1.is_subset(&coords2));
            }

            #[test]
            fn test_into_vec() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(20.0).unwrap(),
                    Real::try_new(30.0).unwrap(),
                ]))
                .unwrap();
                let vec = coords.into_vec();
                assert_eq!(
                    vec,
                    vec![
                        Real::try_new(10.0).unwrap(),
                        Real::try_new(20.0).unwrap(),
                        Real::try_new(30.0).unwrap()
                    ]
                );
            }

            #[test]
            fn test_coords1d_num_points() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(
                    coords.num_points(),
                    PositiveNumPoints1D::try_new(5).unwrap()
                );
            }

            #[test]
            fn test_coords1d_deref() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]))
                .unwrap();
                let slice: &[Real] = &coords;
                assert_eq!(
                    slice,
                    &[
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(3.0).unwrap()
                    ]
                );
            }

            #[test]
            fn test_coords1d_iterator() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                ]))
                .unwrap();
                let collected: Vec<Real> = coords.iter().cloned().collect();
                assert_eq!(
                    collected,
                    vec![
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(3.0).unwrap()
                    ]
                );
            }

            #[test]
            fn test_coords1d_index() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(30.0).unwrap(),
                    Real::try_new(20.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(coords[0], Real::try_new(10.0).unwrap());
                assert_eq!(coords[1], Real::try_new(20.0).unwrap());
                assert_eq!(coords[2], Real::try_new(30.0).unwrap());
            }

            #[test]
            fn test_coords1d_len() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                ]))
                .unwrap();
                assert_eq!(coords.len(), 3);
            }

            #[test]
            fn test_coords1d_edge_cases_for_search() {
                let coords =
                    Coords1D::try_from(SortedSet::from_unsorted(vec![Real::try_new(1.0).unwrap()]))
                        .unwrap();

                // Single element tests
                assert_eq!(coords.find_floor_index(&Real::try_new(0.5).unwrap()), None);
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(1.0).unwrap()),
                    Some(0)
                );
                assert_eq!(
                    coords.find_floor_index(&Real::try_new(1.5).unwrap()),
                    Some(0)
                );

                assert_eq!(
                    coords.find_ceil_index(&Real::try_new(0.5).unwrap()),
                    Some(0)
                );
                assert_eq!(
                    coords.find_ceil_index(&Real::try_new(1.0).unwrap()),
                    Some(0)
                );
                assert_eq!(coords.find_ceil_index(&Real::try_new(1.5).unwrap()), None);
            }

            #[test]
            fn test_coords1d_binary_search() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(7.0).unwrap(),
                    Real::try_new(9.0).unwrap(),
                ]))
                .unwrap();

                // Test exact matches
                assert_eq!(
                    coords.binary_search_by(|probe| probe
                        .partial_cmp(&Real::try_new(3.0).unwrap())
                        .unwrap()),
                    Ok(1)
                );
                assert_eq!(
                    coords.binary_search_by(|probe| probe
                        .partial_cmp(&Real::try_new(7.0).unwrap())
                        .unwrap()),
                    Ok(3)
                );

                // Test non-matches
                assert_eq!(
                    coords.binary_search_by(|probe| probe
                        .partial_cmp(&Real::try_new(4.0).unwrap())
                        .unwrap()),
                    Err(2)
                );
                assert_eq!(
                    coords.binary_search_by(|probe| probe
                        .partial_cmp(&Real::try_new(0.0).unwrap())
                        .unwrap()),
                    Err(0)
                );
                assert_eq!(
                    coords.binary_search_by(|probe| probe
                        .partial_cmp(&Real::try_new(10.0).unwrap())
                        .unwrap()),
                    Err(5)
                );
            }

            #[test]
            fn test_coords1d_find_index() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(7.0).unwrap(),
                    Real::try_new(9.0).unwrap(),
                ]))
                .unwrap();

                // Test exact matches
                assert_eq!(coords.find_index(&Real::try_new(3.0).unwrap()), Some(1));
                assert_eq!(coords.find_index(&Real::try_new(7.0).unwrap()), Some(3));

                // Test non-matches
                assert_eq!(coords.find_index(&Real::try_new(4.0).unwrap()), None);
                assert_eq!(coords.find_index(&Real::try_new(0.0).unwrap()), None);
                assert_eq!(coords.find_index(&Real::try_new(10.0).unwrap()), None);
            }

            #[test]
            fn test_coords1d_partition_point() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                    Real::try_new(6.0).unwrap(),
                    Real::try_new(8.0).unwrap(),
                ]))
                .unwrap();

                assert_eq!(
                    coords.partition_point(|x| x < &Real::try_new(0.0).unwrap()),
                    0
                );
                assert_eq!(
                    coords.partition_point(|x| x < &Real::try_new(3.0).unwrap()),
                    1
                );
                assert_eq!(
                    coords.partition_point(|x| x < &Real::try_new(5.0).unwrap()),
                    2
                );
                assert_eq!(
                    coords.partition_point(|x| x < &Real::try_new(10.0).unwrap()),
                    4
                );
            }

            #[test]
            fn test_coords1d_complex_set_operations() {
                let coords1 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                ]))
                .unwrap();
                let coords2 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                ]))
                .unwrap();
                let coords3 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(5.0).unwrap(),
                    Real::try_new(6.0).unwrap(),
                    Real::try_new(7.0).unwrap(),
                ]))
                .unwrap();

                // Test multiple unions
                let union_12 = coords1.union(&coords2);
                let union_all = union_12.union(&coords3);
                assert_eq!(
                    union_all.deref(),
                    &[
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(3.0).unwrap(),
                        Real::try_new(4.0).unwrap(),
                        Real::try_new(5.0).unwrap(),
                        Real::try_new(6.0).unwrap(),
                        Real::try_new(7.0).unwrap()
                    ]
                );

                // Test subset relationships
                let small =
                    Coords1D::try_from(SortedSet::from_unsorted(vec![Real::try_new(3.0).unwrap()]))
                        .unwrap();
                assert!(small.is_subset(&coords1));
                assert!(small.is_subset(&coords2));
                assert!(!small.is_subset(&coords3));

                // Test empty intersections
                let coords4 = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(10.0).unwrap(),
                    Real::try_new(11.0).unwrap(),
                ]))
                .unwrap();
                assert!(coords1.intersection(&coords4).is_none());
            }

            #[test]
            fn test_coords1d_serialization() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(1.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(5.0).unwrap(),
                ]))
                .unwrap();

                // Serialize to JSON
                let json = serde_json::to_string(&coords).unwrap();

                // Deserialize from JSON
                let deserialized: Coords1D<Real> = serde_json::from_str(&json).unwrap();

                // Verify the deserialized coordinates match the original
                assert_eq!(coords, deserialized);
                assert_eq!(
                    deserialized.deref(),
                    &[
                        Real::try_new(1.0).unwrap(),
                        Real::try_new(2.0).unwrap(),
                        Real::try_new(3.0).unwrap(),
                        Real::try_new(5.0).unwrap()
                    ]
                );
            }

            #[test]
            fn test_coords1d_deserialization_empty_error() {
                // Test that deserializing an empty array fails appropriately
                let json = "[]";
                let result: Result<Coords1D<Real>, _> = serde_json::from_str(json);
                assert!(result.is_err());
            }
        }

        mod points_in_interval {
            use super::*;

            #[test]
            fn test_points_in_interval_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 domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());

                let points = Coords1DInDomain::try_new(coords.clone(), domain.clone()).unwrap();

                assert_eq!(points.coords(), &coords);
                assert_eq!(points.domain(), &domain);
            }

            #[test]
            fn test_points_in_interval_lower_bound_violation() {
                let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![
                    Real::try_new(-0.5).unwrap(),
                    Real::try_new(0.5).unwrap(),
                    Real::try_new(1.0).unwrap(),
                ]))
                .unwrap();
                let domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());

                let pts = Coords1DInDomain::try_new(coords, domain.clone());
                let _c = Real::try_new(-0.5).unwrap();
                matches!(
                    pts,
                    Err(ErrorsCoords1DInDomain::LowerBoundViolation { coord_min: _c, .. })
                );
            }

            #[test]
            fn test_points_in_interval_upper_bound_violation() {
                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.5).unwrap(),
                ]))
                .unwrap();
                let domain =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());

                let pts = Coords1DInDomain::try_new(coords, domain.clone());
                let _c = Real::try_new(1.5).unwrap();
                matches!(
                    pts,
                    Err(ErrorsCoords1DInDomain::UpperBoundViolation { coord_max: _c, .. })
                );
            }
        }
    }

    mod tests_coords1d_in_domain_generic {
        use super::*;

        #[test]
        fn test_owned_coordinates() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            let points_owned = Coords1DInDomain::try_new(coords, domain).unwrap();

            assert_eq!(points_owned.coords().as_ref().as_slice(), &[0.0, 0.5, 1.0]);
            assert_eq!(points_owned.domain(), &IntervalClosed::new(0.0, 1.0));
        }

        #[test]
        fn test_borrowed_coordinates() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            let points_borrowed = Coords1DInDomain::try_new(&coords, domain).unwrap();

            assert_eq!(
                points_borrowed.coords().as_ref().as_slice(),
                &[0.0, 0.5, 1.0]
            );
            assert_eq!(points_borrowed.domain(), &IntervalClosed::new(0.0, 1.0));
        }

        #[test]
        fn test_owned_to_borrowed_conversion() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            let points_owned = Coords1DInDomain::try_new(coords, domain).unwrap();
            let points_borrowed = points_owned.as_borrowed();

            assert_eq!(points_owned.coords(), points_borrowed.coords());
            assert_eq!(points_owned.domain(), points_borrowed.domain());
        }

        #[test]
        fn test_borrowed_to_owned_conversion() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            let points_borrowed = Coords1DInDomain::try_new(&coords, domain).unwrap();
            let points_owned = points_borrowed.to_owned_coords();

            assert_eq!(points_borrowed.coords(), points_owned.coords());
            assert_eq!(points_borrowed.domain(), points_owned.domain());
        }

        #[test]
        fn test_type_aliases() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Test type aliases work correctly
            let _owned: Coords1DInDomainOwned<IntervalClosed<f64>> =
                Coords1DInDomain::try_new(coords.clone(), domain.clone()).unwrap();

            let _borrowed: Coords1DInDomainBorrowed<IntervalClosed<f64>> =
                Coords1DInDomain::try_new(&coords, domain).unwrap();
        }

        #[test]
        fn test_validation_works_with_borrowed() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 2.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Should fail validation for both owned and borrowed
            let result_owned = Coords1DInDomain::try_new(coords.clone(), domain.clone());
            let result_borrowed = Coords1DInDomain::try_new(&coords, domain);

            assert!(result_owned.is_err());
            assert!(result_borrowed.is_err());

            if let Err(ErrorsCoords1DInDomain::UpperBoundViolation { coord_max, .. }) = result_owned
            {
                assert_eq!(coord_max, 2.0);
            } else {
                panic!("Expected UpperBoundViolation");
            }
        }

        #[test]
        fn test_arc_coordinates() {
            use std::sync::Arc;

            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Wrap coordinates in Arc
            let arc_coords = Arc::new(coords);
            let points_arc = Coords1DInDomain::try_new(arc_coords.clone(), domain).unwrap();

            assert_eq!(points_arc.coords().as_ref().as_slice(), &[0.0, 0.5, 1.0]);
            assert_eq!(points_arc.domain(), &IntervalClosed::new(0.0, 1.0));

            // Arc can be shared between multiple Coords1DInDomain instances
            let domain2 = IntervalClosed::new(0.0, 1.0);
            let points_arc2 = Coords1DInDomain::try_new(arc_coords, domain2).unwrap();

            assert_eq!(points_arc.coords(), points_arc2.coords());
        }

        #[test]
        fn test_rc_coordinates() {
            use std::rc::Rc;

            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Wrap coordinates in Rc
            let rc_coords = Rc::new(coords);
            let points_rc = Coords1DInDomain::try_new(rc_coords.clone(), domain).unwrap();

            assert_eq!(points_rc.coords().as_ref().as_slice(), &[0.0, 0.5, 1.0]);
            assert_eq!(points_rc.domain(), &IntervalClosed::new(0.0, 1.0));

            // Rc can be shared between multiple Coords1DInDomain instances
            let domain2 = IntervalClosed::new(0.0, 1.0);
            let points_rc2 = Coords1DInDomain::try_new(rc_coords, domain2).unwrap();

            assert_eq!(points_rc.coords(), points_rc2.coords());
        }

        #[test]
        fn test_box_coordinates() {
            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Wrap coordinates in Box
            let boxed_coords = Box::new(coords);
            let points_box = Coords1DInDomain::try_new(boxed_coords, domain).unwrap();

            assert_eq!(points_box.coords().as_ref().as_slice(), &[0.0, 0.5, 1.0]);
            assert_eq!(points_box.domain(), &IntervalClosed::new(0.0, 1.0));
        }

        #[test]
        fn test_cow_coordinates() {
            use std::borrow::Cow;

            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Cow with owned data
            let cow_owned = Cow::Owned(coords.clone());
            let points_cow_owned = Coords1DInDomain::try_new(cow_owned, domain.clone()).unwrap();

            // Cow with borrowed data
            let cow_borrowed = Cow::Borrowed(&coords);
            let points_cow_borrowed = Coords1DInDomain::try_new(cow_borrowed, domain).unwrap();

            assert_eq!(points_cow_owned.coords(), points_cow_borrowed.coords());
        }

        #[test]
        fn test_smart_pointer_type_aliases() {
            use std::borrow::Cow;
            use std::rc::Rc;
            use std::sync::Arc;

            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            // Test Arc type alias
            let arc_coords = Arc::new(coords.clone());
            let _points_arc: Coords1DInDomainArc<IntervalClosed<f64>> =
                Coords1DInDomain::try_new(arc_coords, domain.clone()).unwrap();

            // Test Rc type alias
            let rc_coords = Rc::new(coords.clone());
            let _points_rc: Coords1DInDomainRc<IntervalClosed<f64>> =
                Coords1DInDomain::try_new(rc_coords, domain.clone()).unwrap();

            // Test Box type alias
            let box_coords = Box::new(coords.clone());
            let _points_box: Coords1DInDomainBox<IntervalClosed<f64>> =
                Coords1DInDomain::try_new(box_coords, domain.clone()).unwrap();

            // Test Cow type alias
            let cow_coords = Cow::Owned(coords);
            let _points_cow: Coords1DInDomainCow<IntervalClosed<f64>> =
                Coords1DInDomain::try_new(cow_coords, domain).unwrap();
        }

        #[test]
        fn test_smart_pointer_trait_methods() {
            use std::sync::Arc;

            let coords = Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();
            let domain = IntervalClosed::new(0.0, 1.0);

            let arc_coords = Arc::new(coords.clone());
            let points_arc = Coords1DInDomain::try_new(arc_coords, domain.clone()).unwrap();
            let points_owned = Coords1DInDomain::try_new(coords.clone(), domain.clone()).unwrap();
            let points_borrowed = Coords1DInDomain::try_new(&coords, domain).unwrap();

            // All should work with HasCoords1D and HasDomain1D traits
            assert_eq!(points_arc.coords(), points_owned.coords());
            assert_eq!(points_arc.coords(), points_borrowed.coords());

            assert_eq!(points_arc.domain(), points_owned.domain());
            assert_eq!(points_arc.domain(), points_borrowed.domain());
        }
    }

    mod tests_transformed_points1d_generic {
        use super::*;
        use crate::traits::{
            TransformedPoints1D, TransformedPoints1DArc, TransformedPoints1DBorrowed,
            TransformedPoints1DOwned, TransformedPoints1DRc,
        };
        use std::sync::Arc;

        // Test struct with owned storage
        struct TestTransformOwned {
            param_points: Coords1DInDomain<IntervalClosed<f64>>,
            phys_points: Coords1DInDomain<IntervalClosed<f64>>,
        }

        impl TransformedPoints1D<Coords1D<f64>, Coords1D<f64>> for TestTransformOwned {
            type ParamDomainType = IntervalClosed<f64>;
            type PhysDomainType = IntervalClosed<f64>;

            fn points_param_interval(
                &self,
            ) -> &Coords1DInDomain<Self::ParamDomainType, Coords1D<f64>> {
                &self.param_points
            }

            fn points_phys_interval(
                &self,
            ) -> &Coords1DInDomain<Self::PhysDomainType, Coords1D<f64>> {
                &self.phys_points
            }
        }

        // Test struct with borrowed storage
        struct TestTransformBorrowed<'a> {
            param_points: Coords1DInDomain<IntervalClosed<f64>, &'a Coords1D<f64>>,
            phys_points: Coords1DInDomain<IntervalClosed<f64>, &'a Coords1D<f64>>,
        }

        impl<'a> TransformedPoints1D<&'a Coords1D<f64>, &'a Coords1D<f64>> for TestTransformBorrowed<'a> {
            type ParamDomainType = IntervalClosed<f64>;
            type PhysDomainType = IntervalClosed<f64>;

            fn points_param_interval(
                &self,
            ) -> &Coords1DInDomain<Self::ParamDomainType, &'a Coords1D<f64>> {
                &self.param_points
            }

            fn points_phys_interval(
                &self,
            ) -> &Coords1DInDomain<Self::PhysDomainType, &'a Coords1D<f64>> {
                &self.phys_points
            }
        }

        // Test struct with Arc storage
        struct TestTransformArc {
            param_points: Coords1DInDomain<IntervalClosed<f64>, Arc<Coords1D<f64>>>,
            phys_points: Coords1DInDomain<IntervalClosed<f64>, Arc<Coords1D<f64>>>,
        }

        impl TransformedPoints1D<Arc<Coords1D<f64>>, Arc<Coords1D<f64>>> for TestTransformArc {
            type ParamDomainType = IntervalClosed<f64>;
            type PhysDomainType = IntervalClosed<f64>;

            fn points_param_interval(
                &self,
            ) -> &Coords1DInDomain<Self::ParamDomainType, Arc<Coords1D<f64>>> {
                &self.param_points
            }

            fn points_phys_interval(
                &self,
            ) -> &Coords1DInDomain<Self::PhysDomainType, Arc<Coords1D<f64>>> {
                &self.phys_points
            }
        }

        #[test]
        fn test_transformed_points_owned_storage() {
            let param_coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![-1.0, 0.0, 1.0])).unwrap();
            let phys_coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();

            let param_domain = IntervalClosed::new(-1.0, 1.0);
            let phys_domain = IntervalClosed::new(0.0, 1.0);

            let transform = TestTransformOwned {
                param_points: Coords1DInDomain::try_new(param_coords, param_domain).unwrap(),
                phys_points: Coords1DInDomain::try_new(phys_coords, phys_domain).unwrap(),
            };

            assert_eq!(
                transform
                    .points_param_interval()
                    .coords()
                    .as_ref()
                    .as_slice(),
                &[-1.0, 0.0, 1.0]
            );
            assert_eq!(
                transform
                    .points_phys_interval()
                    .coords()
                    .as_ref()
                    .as_slice(),
                &[0.0, 0.5, 1.0]
            );
        }

        #[test]
        fn test_transformed_points_borrowed_storage() {
            let param_coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![-1.0, 0.0, 1.0])).unwrap();
            let phys_coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap();

            let param_domain = IntervalClosed::new(-1.0, 1.0);
            let phys_domain = IntervalClosed::new(0.0, 1.0);

            let transform = TestTransformBorrowed {
                param_points: Coords1DInDomain::try_new(&param_coords, param_domain).unwrap(),
                phys_points: Coords1DInDomain::try_new(&phys_coords, phys_domain).unwrap(),
            };

            assert_eq!(
                transform
                    .points_param_interval()
                    .coords()
                    .as_ref()
                    .as_slice(),
                &[-1.0, 0.0, 1.0]
            );
            assert_eq!(
                transform
                    .points_phys_interval()
                    .coords()
                    .as_ref()
                    .as_slice(),
                &[0.0, 0.5, 1.0]
            );
        }

        #[test]
        fn test_transformed_points_arc_storage() {
            let param_coords = Arc::new(
                Coords1D::try_from(SortedSet::from_unsorted(vec![-1.0, 0.0, 1.0])).unwrap(),
            );
            let phys_coords = Arc::new(
                Coords1D::try_from(SortedSet::from_unsorted(vec![0.0, 0.5, 1.0])).unwrap(),
            );

            let param_domain = IntervalClosed::new(-1.0, 1.0);
            let phys_domain = IntervalClosed::new(0.0, 1.0);

            let transform = TestTransformArc {
                param_points: Coords1DInDomain::try_new(param_coords.clone(), param_domain)
                    .unwrap(),
                phys_points: Coords1DInDomain::try_new(phys_coords.clone(), phys_domain).unwrap(),
            };

            assert_eq!(
                transform
                    .points_param_interval()
                    .coords()
                    .as_ref()
                    .as_slice(),
                &[-1.0, 0.0, 1.0]
            );
            assert_eq!(
                transform
                    .points_phys_interval()
                    .coords()
                    .as_ref()
                    .as_slice(),
                &[0.0, 0.5, 1.0]
            );

            // Verify Arc sharing works (the count might be higher due to internal cloning)
            let _param_clone = transform.param_points.clone();
            assert!(Arc::strong_count(&param_coords) >= 2); // At least one in transform, one in clone
        }

        #[test]
        fn test_type_aliases_compile() {
            // These should compile without errors - just testing type aliases
            type _OwnedAlias = TransformedPoints1DOwned<IntervalClosed<f64>, IntervalClosed<f64>>;
            type _BorrowedAlias<'a> =
                TransformedPoints1DBorrowed<'a, IntervalClosed<f64>, IntervalClosed<f64>>;
            type _ArcAlias = TransformedPoints1DArc<IntervalClosed<f64>, IntervalClosed<f64>>;
            type _RcAlias = TransformedPoints1DRc<IntervalClosed<f64>, IntervalClosed<f64>>;

            // Test passes if compilation succeeds
        }
    }

    mod test_logspace {
        use super::*;
        use crate::{logspace, scalars::NumIntervals};
        use try_create::TryNew;

        #[test]
        fn test_logspace_basic() {
            let base = 10.0;
            let exp_domain = IntervalClosed::new(0.0, 2.0); // 10^0 to 10^2
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(2).unwrap());

            assert_eq!(coords.len(), 3);
            assert!((coords[0] - 1.0).abs() < 1e-10); // 10^0 = 1
            assert!((coords[1] - 10.0).abs() < 1e-10); // 10^1 = 10
            assert!((coords[2] - 100.0).abs() < 1e-10); // 10^2 = 100
        }

        #[test]
        fn test_logspace_base_2() {
            let base = 2.0;
            let exp_domain = IntervalClosed::new(0.0, 3.0); // 2^0 to 2^3
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(3).unwrap());

            assert_eq!(coords.len(), 4);
            assert!((coords[0] - 1.0).abs() < 1e-10); // 2^0 = 1
            assert!((coords[1] - 2.0).abs() < 1e-10); // 2^1 = 2
            assert!((coords[2] - 4.0).abs() < 1e-10); // 2^2 = 4
            assert!((coords[3] - 8.0).abs() < 1e-10); // 2^3 = 8
        }

        #[test]
        fn test_logspace_negative_exponents() {
            let base = 10.0;
            let exp_domain = IntervalClosed::new(-2.0, 0.0); // 10^-2 to 10^0
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(2).unwrap());

            assert_eq!(coords.len(), 3);
            assert!((coords[0] - 0.01).abs() < 1e-10); // 10^-2 = 0.01
            assert!((coords[1] - 0.1).abs() < 1e-10); // 10^-1 = 0.1
            assert!((coords[2] - 1.0).abs() < 1e-10); // 10^0 = 1
        }

        #[test]
        fn test_logspace_fractional_base() {
            let base = 0.5;
            let exp_domain = IntervalClosed::new(0.0, 2.0); // 0.5^0 to 0.5^2
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(2).unwrap());

            assert_eq!(coords.len(), 3);
            assert!((coords[0] - 1.0).abs() < 1e-10); // 0.5^0 = 1
            assert!((coords[1] - 0.5).abs() < 1e-10); // 0.5^1 = 0.5
            assert!((coords[2] - 0.25).abs() < 1e-10); // 0.5^2 = 0.25
        }

        #[test]
        fn test_logspace_single_interval() {
            let base = 10.0;
            let exp_domain = IntervalClosed::new(0.0, 1.0);
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(1).unwrap());

            assert_eq!(coords.len(), 2);
            assert!((coords[0] - 1.0).abs() < 1e-10); // 10^0 = 1
            assert!((coords[1] - 10.0).abs() < 1e-10); // 10^1 = 10
        }

        #[test]
        fn test_logspace_many_intervals() {
            let base = 2.0;
            let exp_domain = IntervalClosed::new(0.0, 1.0);
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(10).unwrap());

            assert_eq!(coords.len(), 11);
            assert!((coords[0] - 1.0).abs() < 1e-10); // 2^0 = 1
            assert!((coords[10] - 2.0).abs() < 1e-10); // 2^1 = 2

            // Check that points are increasing
            for i in 0..coords.len() - 1 {
                assert!(coords[i] < coords[i + 1]);
            }
        }

        #[test]
        fn test_logspace_with_half_open_interval() {
            let base = 10.0;
            let exp_domain = IntervalLowerClosedUpperOpen::new(0.0, 2.0);
            let coords = logspace(&base, &exp_domain, &NumIntervals::try_new(2).unwrap());

            assert_eq!(coords.len(), 3);
            // First point is always at lower bound
            more_asserts::assert_lt!((coords[0] - 1.0).abs(), 1e-10);
            // Last point is always at upper bound
            more_asserts::assert_lt!((coords[2] - 100.0).abs(), 1e-10);
        }
    }

    /// Tests for edge cases in find_range_indices
    mod test_find_range_indices_edge_cases {
        use super::*;
        use sorted_vec::partial::SortedSet;

        #[test]
        fn test_find_range_indices_all_coords_above_upper_bound() {
            // Coords: 10, 20, 30, 40, 50
            // Query range: 0 to 5 (all coords are above the range)
            let coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![10.0, 20.0, 30.0, 40.0, 50.0]))
                    .unwrap();
            let range = coords.find_range_indices(&IntervalClosed::new(0.0, 5.0));

            // Should return empty range since no coords are in [0, 5]
            let indices: Vec<usize> = range.collect();
            assert!(indices.is_empty());
        }

        #[test]
        fn test_find_range_indices_all_coords_below_lower_bound() {
            // Coords: 1, 2, 3, 4, 5
            // Query range: 10 to 20 (all coords are below the range)
            let coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 2.0, 3.0, 4.0, 5.0]))
                    .unwrap();
            let range = coords.find_range_indices(&IntervalClosed::new(10.0, 20.0));

            // Should return empty range since no coords are in [10, 20]
            let indices: Vec<usize> = range.collect();
            assert!(indices.is_empty());
        }

        #[test]
        fn test_find_range_indices_partial_overlap() {
            // Coords: 1, 2, 3, 4, 5
            // Query range: 2.5 to 4.5
            let coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 2.0, 3.0, 4.0, 5.0]))
                    .unwrap();
            let range = coords.find_range_indices(&IntervalClosed::new(2.5, 4.5));

            // Should return indices for 3.0 and 4.0
            let indices: Vec<usize> = range.collect();
            assert_eq!(indices, vec![2, 3]);
        }

        #[test]
        fn test_find_range_indices_exact_boundaries() {
            // Coords: 1, 2, 3, 4, 5
            // Query range: 2.0 to 4.0 (exact matches)
            let coords =
                Coords1D::try_from(SortedSet::from_unsorted(vec![1.0, 2.0, 3.0, 4.0, 5.0]))
                    .unwrap();
            let range = coords.find_range_indices(&IntervalClosed::new(2.0, 4.0));

            // Should return indices for 2.0, 3.0, and 4.0
            let indices: Vec<usize> = range.collect();
            assert_eq!(indices, vec![1, 2, 3]);
        }
    }
}