grid1d 0.5.3

#![deny(rustdoc::broken_intra_doc_links)]

//! Bounded interval types with finite, positive-length domains.
//!
//! This module provides all interval types where both bounds are finite. Every type
//! has a well-defined, non-zero length. All items are re-exported at
//! [`crate::intervals`], so `use grid1d::intervals::*` is sufficient.
//!
//! ## Types
//!
//! ### Generic Base Type
//!
//! | Type | Description |
//! |------|------------|
//! | [`IntervalBounded`] | Generic bounded interval parameterised over lower/upper [`crate::bounds::BoundType`] |
//!
//! ### Concrete Type Aliases
//!
//! | Type | Notation | Description |
//! |------|----------|------------|
//! | [`IntervalClosed`] | `[a, b]` | Both bounds included |
//! | [`IntervalOpen`] | `(a, b)` | Both bounds excluded |
//! | [`IntervalLowerClosedUpperOpen`] | `[a, b)` | Lower included, upper excluded |
//! | [`IntervalLowerOpenUpperClosed`] | `(a, b]` | Lower excluded, upper included |
//! | [`IntervalSingleton`] | `{a}` | Single point (zero length) |
//!
//! ### Enum Wrappers
//!
//! | Type | Description |
//! |------|------------|
//! | [`IntervalFinitePositiveLength`] | Enum over the four bounded non-singleton variants |
//! | [`IntervalFiniteLength`] | Extends `IntervalFinitePositiveLength` with [`IntervalSingleton`] |
//! | [`SubIntervalInPartition`] | Enum for grid partition sub-intervals (single/first/middle/last) |
//!
//! ### Traits
//!
//! | Trait | Purpose |
//! |-------|---------|
//! | [`IntervalFinitePositiveLengthTrait`] | Supertrait for all positive-length bounded intervals |
//! | [`IntervalFromBounds`] | Construction from raw bound values: `new()`, `try_new()`, `new_unit_interval()`, `new_symmetric_interval()` |
//!
//! ## Example
//!
//! ```rust
//! use grid1d::intervals::*;
//!
//! let closed = IntervalClosed::new(0.0_f64, 1.0);              // [0, 1]
//! let half_open = IntervalLowerClosedUpperOpen::new(0.0_f64, 1.0); // [0, 1)
//!
//! assert!(closed.contains_point(&1.0));
//! assert!(!half_open.contains_point(&1.0)); // upper bound excluded
//! ```

use crate::{
    GetLowerBoundValue, GetUpperBoundValue,
    bounds::{LowerBoundClosed, UpperBoundClosed},
    grids::traits::Grid1DIntervalBuilder,
    intervals::{
        BoundType, Closed, Contains, ErrorsIntervalConstruction, IntervalBoundRuntime,
        IntervalBoundsRuntime, IntervalHull, IntervalTrait, LowerBound, LowerBoundRuntime, Open,
        UpperBound, UpperBoundRuntime, ValueWithinBound, compute_hull_bounded_intervals,
        operations::IntervalOperations,
    },
};
use duplicate::duplicate_item;
use getset::Getters;
use into_inner::IntoInner;
use num::{One, Zero};
use num_valid::{
    Constants, RealScalar, core::errors::capture_backtrace, scalars::PositiveRealScalar,
};
use serde::{Deserialize, Serialize};
use std::ops::Neg;
use try_create::{New, TryNew};

/// Common interface for intervals with finite, positive, measurable length.
///
/// The [`IntervalFinitePositiveLengthTrait`] trait provides a unified interface for all interval
/// types that have positive length (non-zero measure), enabling generic programming over
/// bounded intervals while excluding both infinite intervals and singleton points.
///
/// ## Mathematical Foundation
///
/// This trait applies to intervals where the upper bound is strictly greater than the lower bound:
/// - **Positive measure**: `μ(I) = |upper - lower| > 0`
/// - **Uncountable cardinality**: Contains infinitely many points
/// - **Non-trivial geometry**: Has well-defined length, midpoint, and interior
/// - **Integration support**: Can support non-zero integrals
///
/// ## Core Operations
///
/// ### Required Methods
/// - [`lower_bound_value`](GetLowerBoundValue::lower_bound_value): Access the lower bound
/// - [`upper_bound_value`](GetUpperBoundValue::upper_bound_value): Access the upper bound
/// - [`into_bounds_pair`](IntervalFinitePositiveLengthTrait::into_bounds_pair): Consume interval to get `(lower, upper)` pair
///
/// ### Provided Methods
///
/// #### Geometric properties
/// - [`length`](IntervalFinitePositiveLengthTrait::length): Compute interval length `|upper - lower|`
/// - [`midpoint`](IntervalFinitePositiveLengthTrait::midpoint): Compute midpoint `(upper + lower) / 2`
/// - [`is_symmetric`](IntervalFinitePositiveLengthTrait::is_symmetric): Test if symmetric around origin (`lower + upper == 0`)
///
/// #### Affine transformations
/// - [`translate`](IntervalFinitePositiveLengthTrait::translate): Shift both bounds by a scalar offset; preserves length and boundary semantics
/// - [`scale`](IntervalFinitePositiveLengthTrait::scale): Scale the half-length around the midpoint by a strictly positive factor
/// - [`expand`](IntervalFinitePositiveLengthTrait::expand): Expand symmetrically by adding a fixed amount to each side; preserves midpoint
///
/// #### Type conversions
/// - [`interior`](IntervalFinitePositiveLengthTrait::interior): Convert to open interval with same bounds
/// - [`closure`](IntervalFinitePositiveLengthTrait::closure): Convert to closed interval with same bounds
///
/// #### Utilities
/// - [`clamp`](IntervalFinitePositiveLengthTrait::clamp): Clamp a scalar value to the interval's `[lower, upper]` range
///
/// ## Usage Examples
///
/// ### Basic Geometric Operations
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// type Real = RealNative64StrictFiniteInDebug;
///
/// let interval = IntervalClosed::new(
///     Real::try_new(-2.0).unwrap(),
///     Real::try_new(3.0).unwrap()
/// );
///
/// // Geometric properties
/// assert_eq!(interval.length().as_ref(), &5.0);           // |3 - (-2)| = 5
/// assert_eq!(interval.midpoint().as_ref(), &0.5);        // (-2 + 3) / 2 = 0.5
/// assert!(!interval.is_symmetric());                      // -2 + 3 ≠ 0
///
/// // Convert to different boundary types
/// let open_version = interval.clone().interior();  // (-2, 3)
/// let closed_version = interval.clone().closure();  // [-2, 3]
/// ```
///
/// ### Generic Programming with Different Interval Types
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealScalar;
///
/// fn analyze_interval<I, T>(interval: &I) -> String
/// where
///     I: IntervalFinitePositiveLengthTrait<RealType = T>,
///     T: RealScalar + std::fmt::Display,
/// {
///     format!(
///         "Interval [{}, {}]: length = {}, symmetric = {}",
///         interval.lower_bound_value(),
///         interval.upper_bound_value(),
///         interval.length(),
///         interval.is_symmetric()
///     )
/// }
///
/// let closed = IntervalClosed::new(0.0, 1.0);
/// let open = IntervalOpen::new(-1.0, 1.0);
///
/// println!("{}", analyze_interval(&closed));
/// // Output: "Interval [0, 1]: length = 1, symmetric = false"
///
/// println!("{}", analyze_interval(&open));
/// // Output: "Interval [-1, 1]: length = 2, symmetric = true"
/// ```
pub trait IntervalFinitePositiveLengthTrait:
    IntervalTrait
    + GetLowerBoundValue<LowerBoundValue = Self::RealType>
    + GetUpperBoundValue<UpperBoundValue = Self::RealType>
{
    /// Return the length of the interval.
    ///
    /// Computes the Euclidean distance between the upper and lower bounds:
    /// `length = upper_bound - lower_bound`. This represents the measure
    /// (Lebesgue measure) of the interval.
    ///
    /// # Mathematical Properties
    ///
    /// - Always positive: `length() > 0` for all intervals implementing this trait
    /// - Scale invariant under positive scaling
    /// - Translation invariant
    /// - Additive for adjacent intervals
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let interval = IntervalClosed::new(-2.0, 3.0);
    /// assert_eq!(interval.length().into_inner(), 5.0);  // |3 - (-2)| = 5
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), requiring only a single subtraction.
    #[inline(always)]
    fn length(&self) -> PositiveRealScalar<Self::RealType> {
        PositiveRealScalar::try_new(self.upper_bound_value().clone() - self.lower_bound_value())
            .unwrap()
    }

    /// Return the midpoint of the interval.
    ///
    /// Computes the center point of the interval using the numerically stable formula
    /// `midpoint = lower_bound + (upper_bound - lower_bound) / 2`. This avoids potential
    /// overflow that would occur with the naive formula `(lower_bound + upper_bound) / 2`
    /// when both bounds are large in magnitude.
    ///
    /// # Mathematical Properties
    ///
    /// - Containment depends on boundary type: closed intervals contain their midpoint,
    ///   while open intervals may or may not (depends on whether midpoint touches boundaries)
    /// - Unique center point: only one point equidistant from both bounds
    /// - Translation equivariant: `translate(interval, c).midpoint() = interval.midpoint() + c`
    /// - Scale equivariant: `scale(interval, k).midpoint() = k * interval.midpoint()`
    ///
    /// # Numerical Stability
    ///
    /// The formula `lower + (upper - lower) / 2` is preferred over `(lower + upper) / 2`
    /// because:
    /// - `upper - lower` is always positive and bounded by the interval length, avoiding
    ///   overflow even when both bounds are large in magnitude (e.g. near `±f64::MAX`)
    /// - The only remaining edge case is an interval of length close to `2 · f64::MAX`,
    ///   which cannot arise when the `RealScalar` type enforces finite bounds
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let closed = IntervalClosed::new(-2.0, 3.0);
    /// assert_eq!(closed.midpoint(), 0.5);  // -2 + (3 - (-2)) / 2 = 0.5
    ///
    /// // Midpoint containment depends on boundary type
    /// assert!(closed.contains_point(&closed.midpoint()));  // Closed contains midpoint
    ///
    /// let open = IntervalOpen::new(-2.0, 3.0);
    /// assert_eq!(open.midpoint(), 0.5);   // Same midpoint value
    /// assert!(open.contains_point(&open.midpoint()));      // Open also contains interior midpoint
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), requiring one subtraction, one multiplication, and one addition.
    #[inline(always)]
    fn midpoint(&self) -> Self::RealType {
        let l = self.lower_bound_value();
        let u = self.upper_bound_value();
        let delta = u.clone() - l;
        delta * Self::RealType::one_div_2() + l
    }

    /// Return `true` if the interval is symmetric about the origin.
    ///
    /// An interval is symmetric if its lower and upper bounds are negatives of each other,
    /// meaning `lower_bound + upper_bound = 0`. This is equivalent to the midpoint being
    /// at the origin.
    ///
    /// # Mathematical Properties
    ///
    /// - Equivalent condition: `midpoint() == 0`
    /// - Preserved under uniform scaling: `scale(symmetric_interval, k)` is symmetric
    /// - Lost under translation: `translate(symmetric_interval, c)` is not symmetric (unless `c = 0`)
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let symmetric = IntervalClosed::new(-2.0, 2.0);
    /// assert!(symmetric.is_symmetric());      // -2 + 2 = 0
    ///
    /// let asymmetric = IntervalClosed::new(1.0, 3.0);
    /// assert!(!asymmetric.is_symmetric());    // 1 + 3 ≠ 0
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), requiring one clone, one sign change and one comparison.
    #[inline]
    fn is_symmetric(&self) -> bool {
        &self.lower_bound_value().clone().neg() == self.upper_bound_value()
    }

    /// Consumes the interval and returns its lower and upper bounds as a pair.
    ///
    /// This is an efficient way to extract both bound values when the interval
    /// object itself is no longer needed, avoiding clones.
    ///
    /// # Returns
    ///
    /// A tuple `(lower_bound, upper_bound)`.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let interval = IntervalClosed::new(1.5, 2.5);
    /// let (lower, upper) = interval.into_bounds_pair();
    /// assert_eq!(lower, 1.5);
    /// assert_eq!(upper, 2.5);
    /// // interval is no longer accessible here
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1) and typically compiles to efficient move operations.
    fn into_bounds_pair(self) -> (Self::RealType, Self::RealType);

    /// Convert this interval to an open interval with the same bounds.
    ///
    /// Creates a new [`IntervalOpen`] with the same lower and upper bound values,
    /// but with both boundaries excluded. This represents the interior of the interval.
    ///
    /// # Mathematical Properties
    ///
    /// - **Bounds preservation**: The numeric values of the bounds remain unchanged
    /// - **Containment**: `interval.interior() ⊆ interval` for all interval types
    /// - **Idempotence**: `interval.interior().interior() = interval.interior()`
    /// - **Interior property**: For any interval I, `interior(I)` is the largest open set contained in I
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let closed = IntervalClosed::new(0.0, 1.0); // [0, 1]
    /// let interior = closed.clone().interior();   // (0, 1)
    ///
    /// // Same bounds, different boundary inclusion
    /// assert_eq!(interior.lower_bound_value(), &0.0);
    /// assert_eq!(interior.upper_bound_value(), &1.0);
    ///
    /// // Different containment behavior
    /// assert!(closed.contains_point(&0.0));       // [0,1] contains 0
    /// assert!(!interior.contains_point(&0.0));    // (0,1) excludes 0
    /// assert!(interior.contains_point(&0.5));     // Both contain interior points
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), involving only a move of the bound values and construction
    /// of a new interval type.
    #[inline(always)]
    fn interior(self) -> IntervalOpen<Self::RealType> {
        let (lower_bound, upper_bound) = self.into_bounds_pair();
        IntervalOpen::new(lower_bound, upper_bound)
    }

    /// Convert this interval to a closed interval with the same bounds.
    ///
    /// Creates a new [`IntervalClosed`] with the same lower and upper bound values,
    /// but with both boundaries included. This represents the closure of the interval.
    ///
    /// # Mathematical Properties
    ///
    /// - **Bounds preservation**: The numeric values of the bounds remain unchanged
    /// - **Containment**: `interval ⊆ interval.closure()` for all interval types
    /// - **Idempotence**: `interval.closure().closure() = interval.closure()`
    /// - **Closure property**: For any interval I, `closure(I)` is the smallest closed set containing I
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let open = IntervalOpen::new(0.0, 1.0);    // (0, 1)
    /// let closure = open.clone().closure();      // [0, 1]
    ///
    /// // Same bounds, different boundary inclusion
    /// assert_eq!(closure.lower_bound_value(), &0.0);
    /// assert_eq!(closure.upper_bound_value(), &1.0);
    ///
    /// // Different containment behavior
    /// assert!(!open.contains_point(&0.0));       // (0,1) excludes 0
    /// assert!(closure.contains_point(&0.0));     // [0,1] contains 0
    /// assert!(closure.contains_point(&0.5));     // Both contain interior points
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), involving only a move of the bound values and construction
    /// of a new interval type.
    #[inline(always)]
    fn closure(self) -> IntervalClosed<Self::RealType> {
        let (lower_bound, upper_bound) = self.into_bounds_pair();
        IntervalClosed::new(lower_bound, upper_bound)
    }

    /// Clamp the value `x` by the lower bound and the upper bound of the current interval.
    ///
    /// If the input value `x` is less than `self.lower_bound()` then this returns `self.lower_bound().
    /// If the input value `x` is greater than `self.upper_bound()` then this returns `self.upper_bound().
    /// Otherwise this returns the value `x`.
    #[inline(always)]
    fn clamp(&self, x: Self::RealType) -> Self::RealType {
        let lower_bound = self.lower_bound_value();
        if &x < lower_bound {
            return lower_bound.clone();
        }
        let upper_bound = self.upper_bound_value();
        if &x > upper_bound {
            return upper_bound.clone();
        }
        x
    }

    /// Translate the interval by the given `amount`.
    ///
    /// Returns a new interval of the same type whose lower and upper bounds are each shifted
    /// by `amount`:
    ///
    /// ```text
    /// translate([a, b], c) = [a + c, b + c]
    /// translate((a, b), c) = (a + c, b + c)
    /// translate([a, b), c) = [a + c, b + c)
    /// translate((a, b], c) = (a + c, b + c]
    /// ```
    ///
    /// # Mathematical Properties
    ///
    /// - **Length preservation**: `translate(I, c).length() == I.length()` for all `c`
    /// - **Midpoint shift**: `translate(I, c).midpoint() == I.midpoint() + c`
    /// - **Boundary semantics preserved**: the open/closed nature of each bound is unchanged
    /// - **Invertibility**: `translate(translate(I, c), -c) == I`
    /// - **Additivity**: `translate(translate(I, a), b) == translate(I, a + b)`
    ///
    /// # Arguments
    ///
    /// - `amount`: The scalar value by which both bounds are shifted.
    ///   Positive values shift right; negative values shift left.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// // Translate a closed interval to the right
    /// let interval = IntervalClosed::new(1.0_f64, 3.0_f64);
    /// let shifted = interval.translate(2.0);
    /// assert_eq!(shifted.lower_bound_value(), &3.0);
    /// assert_eq!(shifted.upper_bound_value(), &5.0);
    /// assert_eq!(shifted.length().into_inner(), interval.length().into_inner()); // length preserved
    ///
    /// // Translate to the left (negative amount)
    /// let left = interval.translate(-1.5);
    /// assert_eq!(left.lower_bound_value(), &-0.5);
    /// assert_eq!(left.upper_bound_value(), &1.5);
    ///
    /// // Boundary semantics are preserved
    /// let half_open = IntervalLowerClosedUpperOpen::new(0.0_f64, 4.0_f64);
    /// let shifted_half = half_open.translate(1.0);
    /// assert_eq!(shifted_half.lower_bound_value(), &1.0);
    /// assert_eq!(shifted_half.upper_bound_value(), &5.0);
    /// assert!(shifted_half.is_lower_bound_closed());
    /// assert!(!shifted_half.is_upper_bound_closed());
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), requiring two additions and one interval construction.
    fn translate(&self, amount: Self::RealType) -> Self;

    /// Scale the interval about its midpoint by a strictly positive `factor`.
    ///
    /// Returns a new interval of the same type whose midpoint is unchanged, but whose
    /// half-length is multiplied by `factor`:
    ///
    /// ```text
    /// let center = midpoint(I)
    /// let new_half_len = (length(I) / 2) * factor
    ///
    /// scale(I, k) = [center - new_half_len,  center + new_half_len]
    /// ```
    ///
    /// The factor must be strictly positive (encoded by [`PositiveRealScalar`]):
    /// - `factor > 1` expands the interval
    /// - `0 < factor < 1` contracts the interval
    /// - `factor = 1` is the identity
    ///
    /// # Mathematical Properties
    ///
    /// - **Midpoint invariance**: `scale(I, k).midpoint() == I.midpoint()` for all `k > 0`
    /// - **Length scaling**: `scale(I, k).length() == I.length() * k`
    /// - **Boundary semantics preserved**: the open/closed nature of each bound is unchanged
    /// - **Identity**: `scale(I, 1)` returns an interval equal to `I`
    /// - **Invertibility**: `scale(scale(I, k), 1/k) == I`
    /// - **Multiplicativity**: `scale(scale(I, a), b) == scale(I, a * b)`
    /// - **Symmetry preservation**: if `I` is symmetric (midpoint = 0), so is `scale(I, k)`
    ///
    /// # Arguments
    ///
    /// - `factor`: A strictly positive scalar by which the half-length is multiplied.
    ///   Construct with [`PositiveRealScalar::try_new`].
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    /// use num_valid::scalars::PositiveRealScalar;
    /// use try_create::TryNew;
    ///
    /// // Expand a closed interval by factor 2 around its midpoint
    /// let interval = IntervalClosed::new(0.0_f64, 4.0_f64); // midpoint = 2, length = 4
    /// let factor = PositiveRealScalar::try_new(2.0_f64).unwrap();
    /// let expanded = interval.scale(&factor);
    /// assert_eq!(expanded.lower_bound_value(), &-2.0); // midpoint - 2*half_len = 2 - 4 = -2
    /// assert_eq!(expanded.upper_bound_value(), &6.0);  // midpoint + 2*half_len = 2 + 4 = 6
    /// assert_eq!(expanded.length().into_inner(), 8.0); // 4 * 2 = 8
    ///
    /// // Contract by factor 0.5
    /// let factor_half = PositiveRealScalar::try_new(0.5_f64).unwrap();
    /// let contracted = interval.scale(&factor_half);
    /// assert_eq!(contracted.lower_bound_value(), &1.0); // 2 - 1 = 1
    /// assert_eq!(contracted.upper_bound_value(), &3.0); // 2 + 1 = 3
    /// assert_eq!(contracted.length().into_inner(), 2.0); // 4 * 0.5 = 2
    ///
    /// // Boundary semantics are preserved
    /// let half_open = IntervalLowerClosedUpperOpen::new(0.0_f64, 4.0_f64);
    /// let scaled = half_open.scale(&factor);
    /// assert!(scaled.is_lower_bound_closed());
    /// assert!(!scaled.is_upper_bound_closed());
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), requiring one midpoint computation, two multiplications,
    /// two additions, and one interval construction.
    fn scale(&self, factor: &PositiveRealScalar<Self::RealType>) -> Self;

    /// Expand the interval symmetrically by a strictly positive `amount` on each side.
    ///
    /// Returns a new interval of the same type whose lower bound is shifted down by `amount`
    /// and whose upper bound is shifted up by `amount`:
    ///
    /// ```text
    /// expand([a, b], ε) = [a - ε,  b + ε]
    /// expand((a, b), ε) = (a - ε,  b + ε)
    /// expand([a, b), ε) = [a - ε,  b + ε)
    /// expand((a, b], ε) = (a - ε,  b + ε]
    /// ```
    ///
    /// Unlike [`scale`](IntervalFinitePositiveLengthTrait::scale), which multiplies the half-length,
    /// `expand` adds a fixed absolute amount to each side.
    ///
    /// # Mathematical Properties
    ///
    /// - **Midpoint invariance**: `expand(I, ε).midpoint() == I.midpoint()` for all `ε > 0`
    /// - **Length increase**: `expand(I, ε).length() == I.length() + 2 * ε`
    /// - **Strict containment**: `I ⊂ expand(I, ε)` (every point of `I` is in the expansion)
    /// - **Boundary semantics preserved**: the open/closed nature of each bound is unchanged
    /// - **Additivity**: `expand(expand(I, a), b) == expand(I, a + b)`
    ///
    /// # Arguments
    ///
    /// - `amount`: A strictly positive scalar added to the upper bound and subtracted from
    ///   the lower bound. Construct with [`PositiveRealScalar::try_new`].
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    /// use num_valid::scalars::PositiveRealScalar;
    /// use try_create::TryNew;
    ///
    /// // Expand a closed interval by 1.0 on each side
    /// let interval = IntervalClosed::new(1.0_f64, 5.0_f64); // midpoint=3, length=4
    /// let eps = PositiveRealScalar::try_new(1.0_f64).unwrap();
    /// let expanded = interval.expand(&eps);
    /// assert_eq!(expanded.lower_bound_value(), &0.0);  // 1 - 1 = 0
    /// assert_eq!(expanded.upper_bound_value(), &6.0);  // 5 + 1 = 6
    /// assert_eq!(expanded.length().into_inner(), 6.0); // 4 + 2*1 = 6
    ///
    /// // Midpoint is unchanged
    /// assert_eq!(expanded.midpoint(), interval.midpoint());
    ///
    /// // Boundary semantics are preserved
    /// let half_open = IntervalLowerClosedUpperOpen::new(1.0_f64, 5.0_f64);
    /// let expanded_half = half_open.expand(&eps);
    /// assert!(expanded_half.is_lower_bound_closed());
    /// assert!(!expanded_half.is_upper_bound_closed());
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1), requiring two additions and one interval construction.
    fn expand(&self, amount: &PositiveRealScalar<Self::RealType>) -> Self;
}
//--------------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------------
/// Trait for constructing a concrete interval type from raw lower and upper bound values.
///
/// Implemented by all bounded finite-positive-length interval types
/// ([`IntervalClosed`], [`IntervalOpen`], [`IntervalLowerClosedUpperOpen`],
/// [`IntervalLowerOpenUpperClosed`]).
/// Use this as an additional bound when generic code needs to *materialise* an interval
/// from raw scalar bounds, rather than only inspect an existing one.
///
/// # Methods
///
/// | Method | Description |
/// |--------|-------------|
/// | [`try_new`](IntervalFromBounds::try_new) | Fallible construction from raw bounds; returns `Err` if `lower >= upper` |
/// | [`new`](IntervalFromBounds::new) | Infallible construction (panics in debug mode on invalid bounds) |
/// | [`new_unit_interval`](IntervalFromBounds::new_unit_interval) | Convenience constructor for `[0, 1]` (or open/half-open equivalent) |
/// | [`new_symmetric_interval`](IntervalFromBounds::new_symmetric_interval) | Convenience constructor for `[-extent, +extent]` centred at zero |
pub trait IntervalFromBounds: IntervalFinitePositiveLengthTrait {
    /// Build and return a new instance of the interval from its lower and upped bound.
    ///
    /// # Arguments
    ///
    /// - `lower_bound`: The lower bound of the interval.
    /// - `upper_bound`: The upper bound of the interval.
    ///
    /// # Errors
    ///
    /// If the lower bound is greater than or equal to the upper bound, this function returns an error
    /// (wrapped by the [`Err`] type).
    /// An error is also returned if any of the bounds is not a finite floating-point number.
    fn try_new(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Result<Self, ErrorsIntervalConstruction<Self::RealType>>;

    /// Build and return a new instance of the interval from its lower and upped bound.
    ///
    /// # Arguments
    ///
    /// - `lower_bound`: The lower bound of the interval.
    /// - `upper_bound`: The upper bound of the interval.
    ///
    /// # Panics
    ///
    /// In debug mode this function internally calls the [`Self::try_new`] method and will panic if
    /// this method returns an error. In release mode the interval is build without any check on the input bounds.
    fn new(lower_bound: Self::RealType, upper_bound: Self::RealType) -> Self;

    /// Build and return the unit interval `[0, 1]` with the boundary semantics of the concrete
    /// implementing type.
    ///
    /// This is a convenience constructor equivalent to `Self::new(T::zero(), T::one())`, where
    /// `T` is `Self::RealType`. It can never fail, since `0 < 1` holds for all supported scalar
    /// types.
    ///
    /// The boundary type (open, closed, half-open) is determined by the concrete implementing type:
    ///
    /// | Type | Result |
    /// |------|--------|
    /// | [`IntervalClosed`] | `[0, 1]` |
    /// | [`IntervalOpen`] | `(0, 1)` |
    /// | [`IntervalLowerClosedUpperOpen`] | `[0, 1)` |
    /// | [`IntervalLowerOpenUpperClosed`] | `(0, 1]` |
    ///
    /// # Examples
    ///
    /// ```
    /// use grid1d::intervals::*;
    ///
    /// // Standard mathematical unit interval [0, 1]
    /// let unit: IntervalClosed<f64> = IntervalClosed::new_unit_interval();
    /// assert_eq!(unit.lower_bound_value(), &0.0);
    /// assert_eq!(unit.upper_bound_value(), &1.0);
    /// assert_eq!(unit.length().into_inner(), 1.0);
    /// assert_eq!(unit.midpoint(), 0.5);
    ///
    /// // Half-open variant: lower included, upper excluded
    /// let half_open: IntervalLowerClosedUpperOpen<f64> =
    ///     IntervalLowerClosedUpperOpen::new_unit_interval();
    /// assert!(half_open.contains_point(&0.0));
    /// assert!(!half_open.contains_point(&1.0));
    /// ```
    fn new_unit_interval() -> Self {
        Self::new(Self::RealType::zero(), Self::RealType::one())
    }

    /// Build and return a symmetric interval `[−extent, +extent]` centred at zero, with the
    /// boundary semantics of the concrete implementing type.
    ///
    /// # Arguments
    ///
    /// - `extent`: the half-length of the interval. Must be strictly positive (enforced by the
    ///   [`PositiveRealScalar`] wrapper).
    ///
    /// This is a convenience constructor equivalent to `Self::new(-extent, extent)`. It can never
    /// fail because `extent > 0` guarantees `-extent < extent`.
    ///
    /// ## Mathematical Properties
    ///
    /// | Property | Value |
    /// |----------|-------|
    /// | Lower bound | `−extent` |
    /// | Upper bound | `+extent` |
    /// | Midpoint | `0` |
    /// | Length | `2 · extent` |
    /// | [`is_symmetric`](IntervalFinitePositiveLengthTrait::is_symmetric) | `true` |
    ///
    /// The boundary type (open, closed, half-open) is determined by the concrete implementing type:
    ///
    /// | Type | Result |
    /// |------|--------|
    /// | [`IntervalClosed`] | `[−extent, extent]` |
    /// | [`IntervalOpen`] | `(−extent, extent)` |
    /// | [`IntervalLowerClosedUpperOpen`] | `[−extent, extent)` |
    /// | [`IntervalLowerOpenUpperClosed`] | `(−extent, extent]` |
    ///
    /// # Examples
    ///
    /// ```
    /// use grid1d::intervals::*;
    /// use num_valid::scalars::PositiveRealScalar;
    /// use try_create::TryNew;
    ///
    /// let extent = PositiveRealScalar::try_new(3.0_f64).unwrap();
    ///
    /// let sym: IntervalClosed<f64> = IntervalClosed::new_symmetric_interval(&extent);
    /// assert_eq!(sym.lower_bound_value(), &-3.0);
    /// assert_eq!(sym.upper_bound_value(), &3.0);
    /// assert_eq!(sym.midpoint(), 0.0);
    /// assert_eq!(sym.length().into_inner(), 6.0);
    /// assert!(sym.is_symmetric());
    /// ```
    fn new_symmetric_interval(extent: &PositiveRealScalar<Self::RealType>) -> Self {
        let extent = extent.as_ref();
        Self::new(extent.clone().neg(), extent.clone())
    }
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType> IntervalFromBounds
    for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
where
    Self: IntervalBoundsRuntime<RealType = RealType>,
    LowerBound<RealType, LowerBoundType>: ValueWithinBound<RealType = RealType>,
    UpperBound<RealType, UpperBoundType>: ValueWithinBound<RealType = RealType>,
{
    /// Build and return a new instance of the interval from its lower and upped bound.
    ///
    /// # Arguments
    ///
    /// - `lower_bound`: The lower bound of the interval.
    /// - `upper_bound`: The upper bound of the interval.
    ///
    /// # Errors
    ///
    /// If the lower bound is greater than or equal to the upper bound, this function returns an error
    /// (wrapped by the [`Err`] type).
    /// An error is also returned if any of the bounds is not a finite floating-point number.
    fn try_new(
        lower_bound: RealType,
        upper_bound: RealType,
    ) -> Result<Self, ErrorsIntervalConstruction<RealType>> {
        if lower_bound < upper_bound {
            let lower_bound = LowerBound::<RealType, LowerBoundType>::new(lower_bound);
            let upper_bound = UpperBound::<RealType, UpperBoundType>::new(upper_bound);

            Ok(Self {
                lower_bound,
                upper_bound,
            })
        } else {
            Err(
                ErrorsIntervalConstruction::LowerBoundGreaterOrEqualThanUpperBound {
                    lower_bound,
                    upper_bound,
                    backtrace: capture_backtrace(),
                },
            )
        }
    }
    fn new(lower_bound: Self::RealType, upper_bound: Self::RealType) -> Self {
        #[cfg(debug_assertions)]
        {
            Self::try_new(lower_bound, upper_bound)
                .expect("Invalid interval bounds for the closed interval!")
        }
        #[cfg(not(debug_assertions))]
        {
            Self {
                lower_bound: LowerBound::<Self::RealType, LowerBoundType>::new(lower_bound),
                upper_bound: UpperBound::<Self::RealType, UpperBoundType>::new(upper_bound),
            }
        }
    }
}
//--------------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------------
///  Container for all intervals with finite, measurable length.
///
/// The [`IntervalFiniteLength`] enum represents the fundamental mathematical distinction
/// between intervals that have positive length (contain infinitely many points) and
/// those with zero length (contain exactly one point). This categorization is essential
/// for mathematical analysis, numerical integration, and measure theory applications.
///
/// ## Mathematical Foundation
///
/// In mathematical analysis, finite intervals fall into two distinct categories:
/// - **Positive measure**: Intervals like `[a,b]`, `(a,b)`, `[a,b)`, `(a,b]` where `a < b`
/// - **Zero measure**: Singleton sets `{a}` or `[a]` containing exactly one point
///
/// This distinction is crucial because:
/// - Positive-length intervals have **non-zero measure** in Lebesgue measure theory
/// - Zero-length intervals (singletons) have **zero measure** but are still measurable sets
/// - Integration over zero-length intervals always yields zero
/// - Positive-length intervals can support non-trivial probability distributions
///
/// ## Generic Over Any [`num_valid::RealScalar`] Type
///
/// **The [`IntervalFiniteLength`] enum works seamlessly with any scalar type implementing [`num_valid::RealScalar`]:**
///
/// | Scalar Type | Performance | Validation | Best For |
/// |-------------|-------------|------------|----------|
/// | [`f64`] | ⚡⚡⚡ Maximum | ❌ None | Trusted input, maximum speed |
/// | [`RealNative64StrictFiniteInDebug`](num_valid::RealNative64StrictFiniteInDebug) | ⚡⚡⚡ **Same as f64** | ✅ **Debug only** | **Recommended for most uses** |
/// | [`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`)|
///
/// ## Enum Variants
///
/// ### [`IntervalFiniteLength::PositiveLength`]
/// Contains intervals where the upper bound is strictly greater than the lower bound:
/// - **Closed**: `[a, b]` where `a < b` - includes both endpoints
/// - **Open**: `(a, b)` where `a < b` - excludes both endpoints  
/// - **Half-open**: `[a, b)` or `(a, b]` where `a < b` - includes one endpoint
/// - **Measure**: Always `> 0` (specifically `|b - a|`)
/// - **Cardinality**: Uncountably infinite points
///
/// ### [`IntervalFiniteLength::ZeroLength`]
/// Contains singleton intervals with exactly one point:
/// - **Notation**: `[a]` or `{a}` - a single point
/// - **Measure**: Always `0`
/// - **Cardinality**: Exactly one point
/// - **Use cases**: Point masses, discrete events, boundary conditions
///
/// ## Construction Patterns
///
/// ### From Concrete Interval Types (Recommended)
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::{New, TryNew};
///
/// // Use the recommended performance-optimal type
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Create positive-length intervals
/// let closed = IntervalClosed::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let finite_positive: IntervalFiniteLength<Real> = closed.into();
///
/// // Create zero-length intervals  
/// let singleton = IntervalSingleton::new(Real::try_new(0.5).unwrap());
/// let finite_zero: IntervalFiniteLength<Real> = singleton.into();
/// ```
///
/// ### Direct Construction
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// // Build positive-length interval directly
/// let positive = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(
///         IntervalClosed::try_new(0.0, 1.0).unwrap()
///     )
/// );
///
/// // Build zero-length interval directly
/// let zero = IntervalFiniteLength::ZeroLength(
///     IntervalSingleton::new(42.0)
/// );
/// ```
///
/// ## Core Operations
///
/// All [`IntervalFiniteLength`] instances support the complete [`IntervalTrait`] interface:
///
/// ### Point Containment Testing
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let positive = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(
///         IntervalClosed::try_new(0.0, 1.0).unwrap())
/// );
/// let zero = IntervalFiniteLength::ZeroLength(IntervalSingleton::new(0.5));
///
/// // Positive-length interval contains many points
/// assert!(positive.contains_point(&0.0));    // Lower bound
/// assert!(positive.contains_point(&0.5));    // Interior point  
/// assert!(positive.contains_point(&1.0));    // Upper bound
/// assert!(!positive.contains_point(&2.0));   // Outside interval
///
/// // Zero-length interval contains only one specific point
/// assert!(zero.contains_point(&0.5));        // Exact match
/// assert!(!zero.contains_point(&0.0));       // Any other point
/// assert!(!zero.contains_point(&1.0));       // Any other point
/// ```
///
/// ### Interval Containment Testing
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let outer = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(
///         IntervalClosed::try_new(0.0, 2.0).unwrap())
/// );
/// let inner = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Open(
///         IntervalOpen::try_new(0.5, 1.5).unwrap())
/// );
/// let singleton = IntervalFiniteLength::ZeroLength(IntervalSingleton::new(1.0));
///
/// assert!(outer.contains_interval(&inner));           // [0,2] contains (0.5,1.5)
/// assert!(outer.contains_interval(&singleton));       // [0,2] contains {1}
/// assert!(inner.contains_interval(&singleton));       // (0.5,1.5) contains {1}
/// assert!(!inner.contains_interval(&outer));          // (0.5,1.5) doesn't contain [0,2]
/// assert!(!singleton.contains_interval(&outer));      // {1} doesn't contain [0,2]
/// ```
///
/// ### Intersection Operations
/// ```rust
/// use grid1d::intervals::*;
///
/// let a = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 2.0))
/// );
/// let b = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(1.0, 3.0))
/// );
///
/// if let Some(intersection) = a.intersection(&b) {
///     // Result: [1.0, 2.0] as general Interval enum
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(positive)) => {
///             println!("Intersection has positive length: {}", positive.length());
///         }
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(singleton)) => {
///             println!("Intersection is a single point: {}", singleton.value());
///         }
///         _ => unreachable!("Finite intervals always have finite intersections"),
///     }
/// }
///
/// // Intersection creating a singleton
/// let c = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 1.0))
/// );
/// let d = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(1.0, 2.0))
/// );
///
/// if let Some(intersection) = c.intersection(&d) {
///     // Result: singleton at point 1.0
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(singleton)) => {
///             assert_eq!(singleton.value(), &1.0);
///         }
///         _ => panic!("Expected singleton intersection"),
///     }
/// }
/// ```
///
/// ## Pattern Matching and Type Analysis
///
/// ### Distinguishing Between Positive and Zero Length
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// fn analyze_finite_interval(interval: &IntervalFiniteLength<f64>) -> String {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             format!(
///                 "Positive-length interval: [{}, {}], length = {}, midpoint = {}",
///                 positive.lower_bound_value(),
///                 positive.upper_bound_value(),
///                 positive.length(),
///                 positive.midpoint()
///             )
///         }
///         IntervalFiniteLength::ZeroLength(singleton) => {
///             format!(
///                 "Zero-length interval (singleton): [{}], measure = 0",
///                 singleton.value()
///             )
///         }
///     }
/// }
///
/// let positive = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(1.0, 3.0))
/// );
/// let zero = IntervalFiniteLength::ZeroLength(IntervalSingleton::new(2.0));
///
/// println!("{}", analyze_finite_interval(&positive));
/// // Output: "Positive-length interval: [1, 3], length = 2, midpoint = 2"
///
/// println!("{}", analyze_finite_interval(&zero));  
/// // Output: "Zero-length interval (singleton): [2], measure = 0"
/// ```
///
/// ### Integration and Measure Theory Applications
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealScalar;
/// use try_create::New;
///
/// /// Compute the measure (length) of any finite interval
/// fn measure<T: RealScalar>(interval: &IntervalFiniteLength<T>) -> T {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => positive.length().into_inner(),
///         IntervalFiniteLength::ZeroLength(_) => T::zero(), // Measure is always 0
///     }
/// }
///
/// /// Integrate a constant function over any finite interval
/// fn integrate_constant<T: RealScalar + Clone>(
///     interval: &IntervalFiniteLength<T>,
///     constant_value: T
/// ) -> T {
///     // Integral = constant × measure of interval
///     constant_value * measure(interval)
/// }
///
/// let positive = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 2.0))
/// );
/// let zero = IntervalFiniteLength::ZeroLength(IntervalSingleton::new(1.0));
///
/// assert_eq!(measure(&positive), 2.0);        // Length of [0,2] is 2
/// assert_eq!(measure(&zero), 0.0);            // Length of {1} is 0
///
/// assert_eq!(integrate_constant(&positive, 5.0), 10.0);  // 5 × 2 = 10
/// assert_eq!(integrate_constant(&zero, 5.0), 0.0);       // 5 × 0 = 0
/// ```
///
/// ## Type Conversion and Recovery
///
/// ### Converting to Specific Types
/// ```rust
/// use grid1d::intervals::*;
///
/// let finite = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 1.0))
/// );
///
/// // Convert to general interval enum
/// let general: Interval<f64> = finite.clone().into();
///
/// // Try to recover positive-length interval
/// match finite {
///     IntervalFiniteLength::PositiveLength(positive) => {
///         println!("Successfully extracted positive-length interval");
///         // Can access all positive-length specific methods
///         println!("Length: {}", positive.length());
///         println!("Midpoint: {}", positive.midpoint());
///     }
///     IntervalFiniteLength::ZeroLength(_) => {
///         println!("This is a singleton, not positive-length");
///     }
/// }
/// ```
///
/// ### Working with Results from Operations
/// ```rust
/// use grid1d::intervals::*;
///
/// fn process_intersection_result(
///     intersection: Option<Interval<f64>>
/// ) -> String {
///     match intersection {
///         Some(Interval::FiniteLength(finite_length)) => {
///             match finite_length {
///                 IntervalFiniteLength::PositiveLength(positive) => {
///                     format!("Non-trivial intersection with length {}", positive.length())
///                 }
///                 IntervalFiniteLength::ZeroLength(singleton) => {
///                     format!("Intervals meet at single point: {}", singleton.value())
///                 }
///             }
///         }
///         Some(Interval::InfiniteLength(_)) => {
///             "Infinite intersection (shouldn't happen with finite inputs)".to_string()
///         }
///         None => "No intersection (disjoint intervals)".to_string(),
///     }
/// }
/// ```
///
/// ## Performance Characteristics
///
/// ### Memory Layout
/// - **Positive-length variants**: Two scalar values + boundary type information (16-24 bytes typical)
/// - **Zero-length variants**: One scalar value (8 bytes typical for f64)
/// - **Enum overhead**: Small discriminant tag (1-8 bytes depending on alignment)
/// - **Total size**: Typically 24-32 bytes depending on scalar type and alignment
///
/// ### Operation Complexity
/// | Operation | Positive Length | Zero Length | Notes |
/// |-----------|----------------|-------------|-------|
/// | **Point queries** | O(1) - 2 comparisons | O(1) - 1 comparison | Constant time regardless of variant |
/// | **Containment tests** | O(1) - boundary logic | O(1) - equality check | No significant difference |
/// | **Intersections** | O(1) - bound arithmetic | O(1) - equality/containment | May allocate result |
/// | **Length calculation** | O(1) - subtraction | O(1) - returns zero | Zero-length always returns 0 |
/// | **Pattern matching** | O(1) - compile-time optimized | O(1) - compile-time optimized | No runtime dispatch cost |
///
/// ### Scalar Type Performance Impact
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::{RealNative64StrictFiniteInDebug, RealNative64StrictFinite, RealScalar};
///
/// // These have identical performance in release builds:
/// type FastFiniteInterval = IntervalFiniteLength<f64>;
/// type OptimalFiniteInterval = IntervalFiniteLength<RealNative64StrictFiniteInDebug>;
///
/// // This has small validation overhead:
/// type SafeFiniteInterval = IntervalFiniteLength<RealNative64StrictFinite>;
///
/// // Performance test function (generic over scalar type)
/// fn benchmark_finite_operations<T: RealScalar + Clone>(
///     interval1: &IntervalFiniteLength<T>,
///     interval2: &IntervalFiniteLength<T>
/// ) {
///     // All operations compile to identical assembly for f64 and RealNative64StrictFiniteInDebug
///     let _contains = interval1.contains_interval(interval2);
///     let _intersection = interval1.intersection(interval2);
///     // Pattern matching is zero-cost in all cases
///     match interval1 {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             let _length = positive.length().into_inner(); // O(1) subtraction
///         }
///         IntervalFiniteLength::ZeroLength(_) => {
///             let _length = T::zero(); // O(1) constant
///         }
///     }
/// }
/// ```
///
/// ## Mathematical Properties and Guarantees
///
/// ### Measure Theory Properties
/// - **Positive-length intervals** always have **positive measure**: `μ([a,b]) = b - a > 0`
/// - **Zero-length intervals** always have **zero measure**: `μ({a}) = 0`
/// - **Additivity**: For disjoint finite intervals, `μ(A ∪ B) = μ(A) + μ(B)`
/// - **Translation invariance**: `μ(A + c) = μ(A)` for any constant `c`
///
/// ### Set Theory Properties
/// - **Positive-length intervals** are **uncountably infinite** sets
/// - **Zero-length intervals** are **finite** sets (cardinality 1)
/// - **Closure**: Intersection of finite intervals is always finite (or empty)
/// - **Boundedness**: All finite intervals are bounded in the real line
///
/// ### Integration Properties
/// ```rust
/// use grid1d::intervals::*;
///
/// // Fundamental theorem: integral over zero-length interval is always zero
/// fn integral_over_singleton<F>(f: F, singleton: &IntervalSingleton<f64>) -> f64
/// where F: Fn(f64) -> f64
/// {
///     0.0 // Always zero regardless of function f
/// }
///
/// // For positive-length intervals, integral depends on function and interval length
/// fn can_have_nonzero_integral<F>(
///     f: F,
///     positive: &IntervalFinitePositiveLength<f64>
/// ) -> bool
/// where F: Fn(f64) -> f64
/// {
///     positive.length().into_inner() > 0.0 // Can be non-zero if function is non-zero somewhere
/// }
/// ```
///
/// ## Integration with Library Ecosystem
///
/// ### With Numerical Integration
/// ```rust
/// use grid1d::intervals::*;
///
/// fn adaptive_integrate(
///     f: impl Fn(f64) -> f64,
///     domain: IntervalFiniteLength<f64>
/// ) -> f64 {
///     match domain {
///         IntervalFiniteLength::PositiveLength(interval) => {
///             // Use sophisticated numerical integration
///             // (Gaussian quadrature, adaptive subdivision, etc.)
///             integrate_positive_length(f, &interval)
///         }
///         IntervalFiniteLength::ZeroLength(_) => {
///             // Integral over single point is always zero
///             0.0
///         }
///     }
/// }
///
/// fn integrate_positive_length(
///     f: impl Fn(f64) -> f64,
///     interval: &IntervalFinitePositiveLength<f64>
/// ) -> f64 {
///     // Placeholder for actual integration algorithm
///     let midpoint = interval.midpoint();
///     f(midpoint) * interval.length().into_inner() // Simple midpoint rule
/// }
/// ```
///
/// ### With Probability Theory
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Represents a probability distribution support
/// struct ProbabilitySupport {
///     support: IntervalFiniteLength<f64>
/// }
///
/// impl ProbabilitySupport {
///     fn new(support: IntervalFiniteLength<f64>) -> Self {
///         Self { support }
///     }
///     
///     fn is_discrete(&self) -> bool {
///         matches!(self.support, IntervalFiniteLength::ZeroLength(_))
///     }
///     
///     fn is_continuous(&self) -> bool {
///         matches!(self.support, IntervalFiniteLength::PositiveLength(_))
///     }
///     
///     fn measure(&self) -> f64 {
///         match &self.support {
///             IntervalFiniteLength::PositiveLength(interval) => interval.length().into_inner(),
///             IntervalFiniteLength::ZeroLength(_) => 0.0,
///         }
///     }
/// }
/// ```
///
/// ## Common Use Cases and Examples
///
/// ### Computational Geometry
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// /// Project a point onto a finite interval
/// fn project_point_onto_interval(
///     point: f64,
///     interval: &IntervalFiniteLength<f64>
/// ) -> f64 {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             // Clamp point to interval bounds
///             point.max(*positive.lower_bound_value())
///                  .min(*positive.upper_bound_value())
///         }
///         IntervalFiniteLength::ZeroLength(singleton) => {
///             // Always project to the single point
///             *singleton.value()
///         }
///     }
/// }
///
/// let interval = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 1.0))
/// );
/// let singleton = IntervalFiniteLength::ZeroLength(IntervalSingleton::new(0.5));
///
/// assert_eq!(project_point_onto_interval(-0.5, &interval), 0.0);  // Clamp to lower bound
/// assert_eq!(project_point_onto_interval(0.7, &interval), 0.7);   // Inside interval
/// assert_eq!(project_point_onto_interval(1.5, &interval), 1.0);   // Clamp to upper bound
///
/// assert_eq!(project_point_onto_interval(-100.0, &singleton), 0.5); // Always singleton value
/// assert_eq!(project_point_onto_interval(100.0, &singleton), 0.5);  // Always singleton value
/// ```
///
/// ### Scientific Computing
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// /// Compute the uniform grid points within a finite interval
/// fn uniform_grid(interval: &IntervalFiniteLength<f64>, num_points: usize) -> Vec<f64> {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             if num_points <= 1 {
///                 vec![positive.midpoint()]
///             } else {
///                 let lower = *positive.lower_bound_value();
///                 let upper = *positive.upper_bound_value();
///                 let step = (upper - lower) / (num_points - 1) as f64;
///                 (0..num_points)
///                     .map(|i| lower + i as f64 * step)
///                     .collect()
///             }
///         }
///         IntervalFiniteLength::ZeroLength(singleton) => {
///             // Grid of a singleton is just the point itself
///             vec![*singleton.value(); num_points.max(1)]
///         }
///     }
/// }
///
/// let interval = IntervalFiniteLength::PositiveLength(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 10.0))
/// );
/// let grid = uniform_grid(&interval, 5);
/// assert_eq!(grid, vec![0.0, 2.5, 5.0, 7.5, 10.0]);
///
/// let singleton = IntervalFiniteLength::ZeroLength(IntervalSingleton::new(42.0));
/// let singleton_grid = uniform_grid(&singleton, 3);
/// assert_eq!(singleton_grid, vec![42.0, 42.0, 42.0]);
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`IntervalFiniteLength`]
/// - **API Design**: When you need to accept any finite interval but want to exclude infinite ones
/// - **Numerical Integration**: Finite intervals have well-defined integrals
/// - **Measure Theory**: When you need to distinguish between zero and positive measure
/// - **Optimization**: When working with bounded feasible regions
///
/// ### When to Use Specific Variants
/// - **PositiveLength**: For continuous domains, integration regions, probability supports
/// - **ZeroLength**: For discrete events, point constraints, boundary conditions
///
/// ### Recommended Patterns
/// ```rust
/// use grid1d::intervals::*;
/// use num::Zero;
/// use num_valid::RealNative64StrictFiniteInDebug;
///
/// // Define type aliases for your application domain
/// type Real = RealNative64StrictFiniteInDebug;
/// type FiniteInterval = IntervalFiniteLength<Real>;
/// type PositiveInterval = IntervalFinitePositiveLength<Real>;
/// type Singleton = IntervalSingleton<Real>;
///
/// // Pattern: Handle both cases explicitly
/// fn compute_weighted_measure(
///     interval: &FiniteInterval,
///     weight: Real
/// ) -> Real {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             weight * positive.length().into_inner() // Non-zero contribution possible
///         }
///         IntervalFiniteLength::ZeroLength(_) => {
///             Real::zero() // Always zero regardless of weight
///         }
///     }
/// }
///
/// // Pattern: Convert to most specific type when needed
/// fn extract_bounds(interval: FiniteInterval) -> Result<(Real, Real), String> {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             Ok(positive.into_bounds_pair()) // Get (lower, upper) bounds
///         }
///         IntervalFiniteLength::ZeroLength(singleton) => {
///             let value = singleton.into_inner();
///             Ok((value.clone(), value)) // Lower and upper are the same
///         }
///     }
/// }
/// ```
///
/// ## Error Handling and Edge Cases
///
/// ### Common Pitfalls and Solutions
/// ```rust
/// use grid1d::intervals::*;
///
/// // WRONG: Assuming all finite intervals have positive length
/// fn naive_midpoint(interval: &IntervalFiniteLength<f64>) -> f64 {
///     // This panics for singletons!
///     // interval.lower_bound_value() + interval.upper_bound_value()) / 2.0
///     unimplemented!("This approach is incorrect")
/// }
///
/// // CORRECT: Handle both cases appropriately  
/// fn robust_midpoint(interval: &IntervalFiniteLength<f64>) -> f64 {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             positive.midpoint() // Proper midpoint calculation
///         }
///         IntervalFiniteLength::ZeroLength(singleton) => {
///             *singleton.value() // The single point is the "midpoint"
///         }
///     }
/// }
///
/// // WRONG: Dividing by length without checking for zero
/// fn naive_normalize(interval: &IntervalFiniteLength<f64>, value: f64) -> f64 {
///     // This can divide by zero for singletons!
///     // value / interval.length()
///     unimplemented!("This approach is incorrect")
/// }
///
/// // CORRECT: Handle zero-length case explicitly
/// fn robust_normalize(
///     interval: &IntervalFiniteLength<f64>,
///     value: f64
/// ) -> Result<f64, &'static str> {
///     match interval {
///         IntervalFiniteLength::PositiveLength(positive) => {
///             Ok(value / positive.length().into_inner())
///         }
///         IntervalFiniteLength::ZeroLength(_) => {
///             if value == 0.0 {
///                 Ok(0.0) // 0/0 case - define as 0
///             } else {
///                 Err("Cannot normalize non-zero value over zero-length interval")
///             }
///         }
///     }
/// }
/// ```
///
/// ## Mathematical Correctness Guarantees
///
/// The [`IntervalFiniteLength`] enum maintains all essential mathematical properties:
///
/// - **Measure consistency**: The measure is always non-negative and finite
/// - **Containment transitivity**: If A ⊆ B and B ⊆ C, then A ⊆ C  
/// - **Intersection commutativity**: A ∩ B = B ∩ A
/// - **Intersection associativity**: (A ∩ B) ∩ C = A ∩ (B ∩ C)
/// - **Monotonicity**: If A ⊆ B, then μ(A) ≤ μ(B)
/// - **Finite additivity**: For disjoint A and B, μ(A ∪ B) = μ(A) + μ(B)
///
/// These properties are preserved regardless of the scalar type used and are guaranteed
/// by the underlying implementations of the [`IntervalTrait`] and associated methods.
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
#[serde(bound(deserialize = "RealType: for<'a> Deserialize<'a>"))]
pub enum IntervalFiniteLength<RealType: RealScalar> {
    /// An interval of ***finite*** and ***positive length***.
    ///
    /// This variant represents intervals where the upper bound is strictly greater than
    /// the lower bound, resulting in a set with uncountably infinite points and positive
    /// measure. Examples include `[a,b]`, `(a,b)`, `[a,b)`, and `(a,b]` where `a < b`.
    ///
    /// **Properties:**
    /// - **Measure**: `> 0` (specifically `|upper - lower|`)
    /// - **Cardinality**: Uncountably infinite
    /// - **Topology**: Connected set with non-empty interior
    /// - **Integration**: Can support non-trivial integrals
    PositiveLength(IntervalFinitePositiveLength<RealType>),

    /// An interval of ***finite*** and ***zero length***.
    ///
    /// This variant represents singleton intervals containing exactly one point.
    /// These are degenerate intervals denoted as `[a]` or `{a}` in mathematical notation.
    ///
    /// **Properties:**
    /// - **Measure**: Always `0`
    /// - **Cardinality**: Exactly `1`
    /// - **Topology**: Isolated point with empty interior
    /// - **Integration**: All integrals evaluate to `0`
    ZeroLength(IntervalSingleton<RealType>),
}

impl<RealType: RealScalar> IntervalFiniteLength<RealType> {
    /// Constructs a finite-length interval from runtime-typed bounds.
    ///
    /// This function handles all combinations of open and closed bounds for finite intervals,
    /// automatically selecting the appropriate concrete interval type. It also handles the special
    /// case of singleton intervals (when both bounds are closed and equal).
    ///
    /// # Arguments
    ///
    /// * `lower_bound` - Runtime-typed lower bound (closed or open)
    /// * `upper_bound` - Runtime-typed upper bound (closed or open)
    ///
    /// # Bound Type Selection
    ///
    /// | Lower | Upper | Result Type | Notation |
    /// |-------|-------|-------------|----------|
    /// | Closed | Closed | [`IntervalClosed`] or [`IntervalSingleton`] | `[a, b]` or `{a}` |
    /// | Closed | Open | [`IntervalLowerClosedUpperOpen`] | `[a, b)` |
    /// | Open | Closed | [`IntervalLowerOpenUpperClosed`] | `(a, b]` |
    /// | Open | Open | [`IntervalOpen`] | `(a, b)` |
    ///
    /// # Returns
    ///
    /// * `Ok(IntervalFiniteLength)` - Successfully constructed interval
    /// * `Err(ErrorsIntervalConstruction)` - If `lower_bound ≥ upper_bound` (except for singleton case)
    ///
    /// # Examples
    ///
    /// ```ignore
    /// // Note: This is a private function used internally
    /// use grid1d::intervals::*;
    /// use grid1d::bounds::*;
    /// use try_create::New;
    ///
    /// // Closed interval [0, 1]
    /// let lower: LowerBoundRuntime<f64> = LowerBoundClosed::new(0.0).into();
    /// let upper: UpperBoundRuntime<f64> = UpperBoundClosed::new(1.0).into();
    /// let interval = IntervalFiniteLength::try_from_runtime_bounds(lower, upper).unwrap();
    /// ```
    ///
    /// # Singleton Detection
    ///
    /// When both bounds are closed and have equal values, automatically creates an
    /// [`IntervalSingleton`] instead of a zero-length [`IntervalClosed`].
    ///
    /// # Design Notes
    ///
    /// This function is marked `#[allow(dead_code)]` as it serves as infrastructure for
    /// the universal `Interval::try_new` constructor. It may be promoted to `pub(crate)`
    /// if needed elsewhere in the codebase.
    #[allow(dead_code)]
    pub(crate) fn try_from_runtime_bounds(
        lower_bound: LowerBoundRuntime<RealType>,
        upper_bound: UpperBoundRuntime<RealType>,
    ) -> Result<Self, ErrorsIntervalConstruction<RealType>> {
        // Delegate to the appropriate concrete interval type's try_new() method
        // This ensures consistent validation and maintains a single source of truth
        match (lower_bound, upper_bound) {
            (LowerBoundRuntime::Closed(lower_bound), UpperBoundRuntime::Closed(upper_bound)) => {
                if lower_bound.as_ref() == upper_bound.as_ref() {
                    Ok(IntervalSingleton::new(lower_bound.into_inner()).into())
                } else {
                    Ok(IntervalClosed::try_from_static_bounds(lower_bound, upper_bound)?.into())
                }
            }
            (LowerBoundRuntime::Closed(lower_bound), UpperBoundRuntime::Open(upper_bound)) => Ok(
                IntervalLowerClosedUpperOpen::try_from_static_bounds(lower_bound, upper_bound)?
                    .into(),
            ),
            (LowerBoundRuntime::Open(lower_bound), UpperBoundRuntime::Closed(upper_bound)) => Ok(
                IntervalLowerOpenUpperClosed::try_from_static_bounds(lower_bound, upper_bound)?
                    .into(),
            ),
            (LowerBoundRuntime::Open(lower_bound), UpperBoundRuntime::Open(upper_bound)) => {
                Ok(IntervalOpen::try_from_static_bounds(lower_bound, upper_bound)?.into())
            }
        }
    }
}

impl<RealType: RealScalar> IntervalOperations for IntervalFiniteLength<RealType> {}

impl<RealType: RealScalar> IntervalBoundsRuntime for IntervalFiniteLength<RealType> {
    type RealType = RealType;

    #[inline]
    fn lower_bound_runtime(&self) -> Option<LowerBoundRuntime<RealType>> {
        match self {
            IntervalFiniteLength::PositiveLength(interval) => interval.lower_bound_runtime(),
            IntervalFiniteLength::ZeroLength(interval) => interval.lower_bound_runtime(),
        }
    }

    #[inline]
    fn upper_bound_runtime(&self) -> Option<UpperBoundRuntime<RealType>> {
        match self {
            IntervalFiniteLength::PositiveLength(interval) => interval.upper_bound_runtime(),
            IntervalFiniteLength::ZeroLength(interval) => interval.upper_bound_runtime(),
        }
    }
}

impl<RealType: RealScalar> Contains for IntervalFiniteLength<RealType> {
    /// Returns `true` if the 1D point/coordinate `x` is contained in `self`.
    #[inline(always)]
    fn contains_point(&self, x: &RealType) -> bool {
        match self {
            IntervalFiniteLength::PositiveLength(interval) => interval.contains_point(x),
            IntervalFiniteLength::ZeroLength(interval) => interval.contains_point(x),
        }
    }
}

impl<RealType: RealScalar> IntervalTrait for IntervalFiniteLength<RealType> {}
//--------------------------------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------------------------------
/// Generic bounded interval structure for intervals with two finite endpoints.
///
/// The [`IntervalBounded`] struct represents the most general form of finite intervals,
/// where both the lower and upper bounds are finite real numbers. This generic structure
/// can represent all possible boundary combinations (open/closed) through its type parameters,
/// providing maximum flexibility while maintaining type safety.
///
/// ## Generic Parameters
///
/// - `RealType`: Must implement [`RealScalar`] for mathematical operations and comparisons
/// - `LowerBoundType`: Either [`Closed`] or [`Open`] to specify lower boundary inclusion
/// - `UpperBoundType`: Either [`Closed`] or [`Open`] to specify upper boundary inclusion
///
/// ## Boundary Type Combinations
///
/// The type parameters allow representing all possible finite interval types:
///
/// | Lower | Upper | Notation | Mathematical Set |
/// |-------|-------|----------|------------------|
/// | `Closed` | `Closed` | `[a, b]` | `{x : a ≤ x ≤ b}` |
/// | `Open` | `Open` | `(a, b)` | `{x : a < x < b}` |
/// | `Closed` | `Open` | `[a, b)` | `{x : a ≤ x < b}` |
/// | `Open` | `Closed` | `(a, b]` | `{x : a < x ≤ b}` |
///
/// ## Mathematical Properties
///
/// All bounded intervals share these fundamental properties:
/// - **Finite measure**: `μ([a,b]) = b - a > 0` (for valid intervals)
/// - **Compactness**: Only closed intervals `[a,b]` are compact
/// - **Connectedness**: All bounded intervals are connected sets
/// - **Boundedness**: Both bounds are finite real numbers
///
/// ## Type Aliases
///
/// The library provides convenient type aliases for common cases:
/// - [`IntervalClosed`] = `IntervalBounded<RealType, Closed, Closed>`
/// - [`IntervalOpen`] = `IntervalBounded<RealType, Open, Open>`
/// - [`IntervalLowerClosedUpperOpen`] = `IntervalBounded<RealType, Closed, Open>`
/// - [`IntervalLowerOpenUpperClosed`] = `IntervalBounded<RealType, Open, Closed>`
///
/// ## Construction
///
/// ```rust
/// use grid1d::{
///     bounds::{Closed, Open},
///     intervals::bounded::*,
/// };
///
/// // Direct construction (requires explicit type parameters)
/// let closed: IntervalBounded<f64, Closed, Closed> =
///     IntervalBounded::try_new(0.0, 1.0).unwrap();
///
/// // Preferred: use type aliases
/// let closed = IntervalClosed::try_new(0.0, 1.0).unwrap();
/// let open = IntervalOpen::try_new(0.0, 1.0).unwrap();
/// ```
///
/// ## Performance Characteristics
///
/// - **Memory**: Two scalar values + bound type information (compile-time)
/// - **Operations**: O(1) for all basic operations (point queries, intersections)
/// - **Type safety**: Boundary semantics encoded at compile time
/// - **Zero-cost abstractions**: Boundary types have no runtime overhead
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize, Getters)]
#[serde(bound(deserialize = "RealType: for<'a> Deserialize<'a>"))]
pub struct IntervalBounded<
    RealType: RealScalar,
    LowerBoundType: BoundType,
    UpperBoundType: BoundType,
> {
    /// The lower bound of the interval, included in the interval.
    ///
    /// This value represents the smallest point contained in the interval.
    /// For a closed lower bound `[a`, any point `x` in the interval satisfies `x ≥ a`.
    /// For an open lower bound `(a`, any point `x` in the interval satisfies `x > a`.
    #[getset(get = "pub")]
    lower_bound: LowerBound<RealType, LowerBoundType>,

    /// The upper bound of the interval, included in the interval.
    ///
    /// This value represents the largest point contained in the interval.
    /// For a closed upper bound `b]`, any point `x` in the interval satisfies `x ≤ b`.
    /// For an open upper bound `b)`, any point `x` in the interval satisfies `x < b`.
    #[getset(get = "pub")]
    upper_bound: UpperBound<RealType, UpperBoundType>,
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType> GetLowerBoundValue
    for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
{
    type LowerBoundValue = RealType;

    #[inline(always)]
    fn lower_bound_value(&self) -> &RealType {
        &self.lower_bound.value
    }

    #[inline(always)]
    fn is_lower_bound_closed(&self) -> bool {
        LowerBoundType::is_closed()
    }
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType> GetUpperBoundValue
    for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
{
    type UpperBoundValue = RealType;

    #[inline(always)]
    fn upper_bound_value(&self) -> &RealType {
        &self.upper_bound.value
    }

    #[inline(always)]
    fn is_upper_bound_closed(&self) -> bool {
        UpperBoundType::is_closed()
    }
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType>
    IntervalFinitePositiveLengthTrait for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
where
    Self: IntervalBoundsRuntime<RealType = RealType>,
    LowerBound<RealType, LowerBoundType>: ValueWithinBound<RealType = RealType>,
    UpperBound<RealType, UpperBoundType>: ValueWithinBound<RealType = RealType>,
{
    /// Return the inner bounds of the interval as a pair in which the first entry of the pair
    /// is the lower bound and the second entry of the pair is the upper bound.
    ///
    /// # Note
    /// This function consumes `Self` (no cloning of the bounds is needed).
    #[inline(always)]
    fn into_bounds_pair(self) -> (Self::RealType, Self::RealType) {
        (self.lower_bound.into_inner(), self.upper_bound.into_inner())
    }

    fn translate(&self, amount: Self::RealType) -> Self {
        let lower_bound = self.lower_bound_value().clone() + amount.clone();
        let upper_bound = self.upper_bound_value().clone() + amount;
        Self::new(lower_bound, upper_bound)
    }

    fn scale(&self, factor: &PositiveRealScalar<Self::RealType>) -> Self {
        let center = self.midpoint();
        let half_len = Self::RealType::one_div_2() * self.length().as_ref() * factor.as_ref();
        Self::new(center.clone() - &half_len, center.clone() + half_len)
    }

    fn expand(&self, amount: &PositiveRealScalar<Self::RealType>) -> Self {
        let lower_bound = self.lower_bound_value().clone() - amount.as_ref();
        let upper_bound = self.upper_bound_value().clone() + amount.as_ref();
        Self::new(lower_bound, upper_bound)
    }
}

impl<RealType: RealScalar> GetLowerBoundValue for IntervalFinitePositiveLength<RealType> {
    type LowerBoundValue = RealType;

    #[inline(always)]
    fn lower_bound_value(&self) -> &RealType {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.lower_bound_value(),
            IntervalFinitePositiveLength::Open(interval) => interval.lower_bound_value(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.lower_bound_value()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.lower_bound_value()
            }
        }
    }

    #[inline(always)]
    fn is_lower_bound_closed(&self) -> bool {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.is_lower_bound_closed(),
            IntervalFinitePositiveLength::Open(interval) => interval.is_lower_bound_closed(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.is_lower_bound_closed()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.is_lower_bound_closed()
            }
        }
    }
}

impl<RealType: RealScalar> GetUpperBoundValue for IntervalFinitePositiveLength<RealType> {
    type UpperBoundValue = RealType;

    #[inline(always)]
    fn upper_bound_value(&self) -> &RealType {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.upper_bound_value(),
            IntervalFinitePositiveLength::Open(interval) => interval.upper_bound_value(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.upper_bound_value()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.upper_bound_value()
            }
        }
    }

    #[inline(always)]
    fn is_upper_bound_closed(&self) -> bool {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.is_upper_bound_closed(),
            IntervalFinitePositiveLength::Open(interval) => interval.is_upper_bound_closed(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.is_upper_bound_closed()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.is_upper_bound_closed()
            }
        }
    }
}

impl<RealType: RealScalar> IntervalFinitePositiveLengthTrait
    for IntervalFinitePositiveLength<RealType>
{
    /// Return the inner bounds of the interval as a pair in which the first entry of the pair
    /// is the lower bound and the second entry of the pair is the upper bound.
    ///
    /// # Note
    /// This function consumes `Self` (no cloning of the bounds is needed).
    #[inline(always)]
    fn into_bounds_pair(self) -> (Self::RealType, Self::RealType) {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.into_bounds_pair(),
            IntervalFinitePositiveLength::Open(interval) => interval.into_bounds_pair(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.into_bounds_pair()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.into_bounds_pair()
            }
        }
    }

    fn translate(&self, amount: Self::RealType) -> Self {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => {
                interval.translate(amount.clone()).into()
            }
            IntervalFinitePositiveLength::Open(interval) => {
                interval.translate(amount.clone()).into()
            }
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.translate(amount.clone()).into()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.translate(amount.clone()).into()
            }
        }
    }

    fn scale(&self, factor: &PositiveRealScalar<Self::RealType>) -> Self {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.scale(factor).into(),
            IntervalFinitePositiveLength::Open(interval) => interval.scale(factor).into(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.scale(factor).into()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.scale(factor).into()
            }
        }
    }
    fn expand(&self, amount: &PositiveRealScalar<Self::RealType>) -> Self {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.expand(amount).into(),
            IntervalFinitePositiveLength::Open(interval) => interval.expand(amount).into(),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.expand(amount).into()
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.expand(amount).into()
            }
        }
    }
}

#[duplicate_item(
    interval_type                  lower_bound_type upper_bound_type;
    [IntervalClosed]               [Closed]         [Closed];
    [IntervalOpen]                 [Open]           [Open];
    [IntervalLowerClosedUpperOpen] [Closed]         [Open];
    [IntervalLowerOpenUpperClosed] [Open]           [Closed];
)]
impl<RealType: RealScalar> IntervalBoundsRuntime for interval_type<RealType> {
    type RealType = RealType;

    #[inline(always)]
    fn lower_bound_runtime(&self) -> Option<LowerBoundRuntime<RealType>> {
        Some(LowerBoundRuntime::lower_bound_type(
            self.lower_bound().clone(),
        ))
    }

    #[inline(always)]
    fn upper_bound_runtime(&self) -> Option<UpperBoundRuntime<RealType>> {
        Some(UpperBoundRuntime::upper_bound_type(
            self.upper_bound().clone(),
        ))
    }
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType> IntervalTrait
    for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
where
    Self: IntervalBoundsRuntime<RealType = RealType>,
    LowerBound<RealType, LowerBoundType>: ValueWithinBound<RealType = RealType>,
    UpperBound<RealType, UpperBoundType>: ValueWithinBound<RealType = RealType>,
{
}

//--------------------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
///  A closed interval `[a, b]` where both endpoints are included.
///
/// The [`IntervalClosed`] represents the most common type of interval in mathematics,
/// where both the lower bound `a` and upper bound `b` are part of the interval.
/// Every real number `x` satisfying `a ≤ x ≤ b` belongs to this interval.
///
/// ## Mathematical Properties
///
/// - **Notation**: `[a, b]`
/// - **Definition**: `{x ∈ ℝ : a ≤ x ≤ b}`
/// - **Endpoints**: Both included
/// - **Compactness**: Always compact (closed and bounded)
/// - **Connectedness**: Always connected (no gaps)
///
/// ## Use Cases
///
/// Closed intervals are ideal when:
/// - You need to include boundary values (e.g., `[0., 1.]` for probabilities)
/// - Working with integration bounds where endpoints matter
/// - Defining domains where the boundary is part of the solution
/// - Representing physical quantities with inclusive limits
///
/// ## Construction
///
/// ```rust
/// use grid1d::intervals::*;
///
/// // Infallible construction (panics in debug mode if invalid)
/// let interval = IntervalClosed::new(0.0, 1.0);
///
/// // Fallible construction (returns Result)
/// let result = IntervalClosed::try_new(0.0, 1.0).unwrap();
/// ```
///
/// ## Key Operations
///
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalClosed::new(-1.0, 1.0);
///
/// // Point containment (includes endpoints)
/// assert!(interval.contains_point(&-1.0));  // Lower bound included
/// assert!(interval.contains_point(&0.0));   // Interior point
/// assert!(interval.contains_point(&1.0));   // Upper bound included
/// assert!(!interval.contains_point(&2.0));  // Outside interval
///
/// // Geometric properties
/// assert_eq!(interval.length().into_inner(), 2.0);       // |b - a|
/// assert_eq!(interval.midpoint(), 0.0);    // (a + b) / 2
/// assert!(interval.is_symmetric());         // True when a = -b
///
/// // Clamping values to interval bounds
/// assert_eq!(interval.clamp(-2.0), -1.0);   // Below range
/// assert_eq!(interval.clamp(0.5), 0.5);     // Within range
/// assert_eq!(interval.clamp(2.0), 1.0);     // Above range
/// ```
///
/// ## Comparison with Other Interval Types
///
/// | Property | Closed `[a,b]` | Open `(a,b)` | Half-open `[a,b)` |
/// |----------|----------------|--------------|-------------------|
/// | Contains `a` | ✓ | ✗ | ✓ |
/// | Contains `b` | ✓ | ✗ | ✗ |
/// | Length | `b - a` | `b - a` | `b - a` |
/// | Compactness | Always | Never | Never |
/// | Supremum | `b` (achieved) | `b` (not achieved) | `b` (not achieved) |
/// | Infimum | `a` (achieved) | `a` (not achieved) | `a` (achieved) |
///
/// ## Performance Characteristics
///
/// - **Memory**: Two values of type `RealType` plus bound type information
/// - **Point queries**: O(1) - two comparisons
/// - **Intersection**: O(1) - bound arithmetic
/// - **Construction**: O(1) - single comparison in debug mode
///
/// ## Error Conditions
///
/// Construction fails when:
/// - `lower_bound > upper_bound` (invalid interval)
/// - Either bound is not finite (NaN or infinite)
///
/// ```rust
/// use grid1d::intervals::*;
///
/// // These will fail
/// assert!(IntervalClosed::try_new(1.0, 0.0).is_err());     // Inverted bounds
/// assert!(IntervalClosed::try_new(f64::NAN, 1.0).is_err()); // NaN bound
/// ```
///
/// ## Advanced Usage
///
/// ### Generic Programming
///
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealScalar;
///
/// fn process_closed_interval<T: RealScalar>(interval: &IntervalClosed<T>) {
///     // Work with any numeric type that implements RealScalar
///     println!("Length: {}", interval.length());
/// }
/// ```
///
/// ### Conversion to Other Types
///
/// ```rust
/// use grid1d::intervals::*;
///
/// let closed = IntervalClosed::new(0.0, 1.0);
///
/// // Convert to general interval enum
/// let general: Interval<f64> = closed.clone().into();
///
/// // Convert to finite positive length category
/// let finite: IntervalFinitePositiveLength<f64> = closed.into();
/// ```
///
/// ## Mathematical Examples
///
/// ### Unit Interval
/// ```rust
/// use grid1d::intervals::*;
///
/// let unit = IntervalClosed::new(0.0, 1.0); // [0., 1.]
/// ```
///
/// ### Symmetric Interval
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::{RealNative64StrictFinite, Constants};
///
/// let pi = RealNative64StrictFinite::pi();
///
/// let symmetric = IntervalClosed::new(-pi.clone(), pi); // [-π, π]
/// assert!(symmetric.is_symmetric());
/// ```
pub type IntervalClosed<RealType> = IntervalBounded<RealType, Closed, Closed>;

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType>
    IntervalBounded<RealType, LowerBoundType, UpperBoundType>
{
    /// Build and return a new instance of the interval from two `IntervalBound` objects.
    ///
    /// This constructor accepts pre-constructed [`LowerBound`] and [`UpperBound`] objects,
    /// allowing for more flexible interval construction when bounds are already available
    /// in their typed form.
    ///
    /// # Arguments
    ///
    /// - `lower_bound`: A [`LowerBound`] object representing the interval's lower boundary.
    /// - `upper_bound`: An [`UpperBound`] object representing the interval's upper boundary.
    ///
    /// # Errors
    ///
    /// Returns an error if the lower bound value is greater than or equal to the upper bound value.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    /// use grid1d::bounds::*;
    /// use try_create::New;
    ///
    /// // Create bounds separately
    /// let lower = LowerBoundClosed::new(0.0);
    /// let upper = UpperBoundClosed::new(1.0);
    ///
    /// // Construct interval from bounds
    /// let interval = IntervalClosed::try_from_static_bounds(lower, upper).unwrap();
    /// assert_eq!(interval.lower_bound_value(), &0.0);
    /// assert_eq!(interval.upper_bound_value(), &1.0);
    /// ```
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, bounds::*};
    /// use try_create::New;
    ///
    /// // Create half-open interval [0, 1) from bounds
    /// let lower = LowerBoundClosed::new(0.0);
    /// let upper = UpperBoundOpen::new(1.0);
    ///
    /// let interval = IntervalLowerClosedUpperOpen::try_from_static_bounds(lower, upper).unwrap();
    /// assert!(interval.contains_point(&0.0));  // Lower bound included
    /// assert!(!interval.contains_point(&1.0)); // Upper bound excluded
    /// ```
    ///
    /// # Type Safety
    ///
    /// This method maintains type safety by requiring that the bound types match the interval's
    /// generic parameters. For example, an [`IntervalClosed`] requires both bounds to be closed:
    ///
    /// ```compile_fail
    /// use grid1d::{*, bounds::*, intervals::*};
    /// use try_create::New;
    ///
    /// let lower = LowerBoundClosed::new(0.0);
    /// let upper = UpperBoundOpen::new(1.0);  // Wrong type!
    ///
    /// // This will not compile because IntervalClosed expects closed bounds
    /// let interval = IntervalClosed::try_from_static_bounds(lower, upper);
    /// ```
    pub fn try_from_static_bounds(
        lower_bound: LowerBound<RealType, LowerBoundType>,
        upper_bound: UpperBound<RealType, UpperBoundType>,
    ) -> Result<Self, ErrorsIntervalConstruction<RealType>> {
        if lower_bound.as_ref() < upper_bound.as_ref() {
            Ok(Self {
                lower_bound,
                upper_bound,
            })
        } else {
            Err(
                ErrorsIntervalConstruction::LowerBoundGreaterOrEqualThanUpperBound {
                    lower_bound: lower_bound.as_ref().clone(),
                    upper_bound: upper_bound.as_ref().clone(),
                    backtrace: capture_backtrace(),
                },
            )
        }
    }

    /// Build and return a new instance of the interval from two `IntervalBound` objects.
    ///
    /// This is the infallible version of [`Self::try_new`]. It constructs an interval
    /// from pre-existing bound objects without validation.
    ///
    /// # Arguments
    ///
    /// - `lower_bound`: A [`LowerBound`] object representing the interval's lower boundary.
    /// - `upper_bound`: An [`UpperBound`] object representing the interval's upper boundary.
    ///
    /// # Panics
    ///
    /// In debug mode, this function internally calls [`Self::try_new`] and will panic
    /// if that method returns an error. In release mode, the interval is constructed without
    /// any validation.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    /// use grid1d::bounds::*;
    /// use try_create::New;
    ///
    /// let lower = LowerBoundOpen::new(0.0);
    /// let upper = UpperBoundOpen::new(1.0);
    ///
    /// let interval = IntervalOpen::from_static_bounds(lower, upper);
    /// assert!(!interval.contains_point(&0.0)); // Both bounds excluded
    /// assert!(!interval.contains_point(&1.0));
    /// assert!(interval.contains_point(&0.5));  // Interior point included
    /// ```
    #[inline(always)]
    pub fn from_static_bounds(
        lower_bound: LowerBound<RealType, LowerBoundType>,
        upper_bound: UpperBound<RealType, UpperBoundType>,
    ) -> Self {
        #[cfg(debug_assertions)]
        {
            Self::try_from_static_bounds(lower_bound, upper_bound)
                .expect("Invalid interval bounds!")
        }
        #[cfg(not(debug_assertions))]
        {
            Self {
                lower_bound,
                upper_bound,
            }
        }
    }
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType> IntervalOperations
    for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
where
    Self: IntervalBoundsRuntime<RealType = RealType>,
    LowerBound<RealType, LowerBoundType>: ValueWithinBound<RealType = RealType>,
    UpperBound<RealType, UpperBoundType>: ValueWithinBound<RealType = RealType>,
{
}

impl<RealType: RealScalar, LowerBoundType: BoundType, UpperBoundType: BoundType> Contains
    for IntervalBounded<RealType, LowerBoundType, UpperBoundType>
where
    Self: IntervalBoundsRuntime<RealType = RealType>,
    LowerBound<RealType, LowerBoundType>: ValueWithinBound<RealType = RealType>,
    UpperBound<RealType, UpperBoundType>: ValueWithinBound<RealType = RealType>,
{
    #[inline(always)]
    fn contains_point(&self, x: &RealType) -> bool {
        self.lower_bound.value_within_bound(x) && self.upper_bound.value_within_bound(x)
    }
}

//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
///  An open interval `(a, b)` where both endpoints are excluded.
///
/// The [`IntervalOpen`] represents intervals where neither the lower bound `a` nor
/// the upper bound `b` are part of the interval. Every real number `x` satisfying
/// `a < x < b` belongs to this interval, but the boundary points themselves are excluded.
///
/// ## Mathematical Properties
///
/// - **Notation**: `(a, b)`
/// - **Definition**: `{x ∈ ℝ : a < x < b}`
/// - **Endpoints**: Both excluded
/// - **Topology**: Open set in the real line topology
/// - **Compactness**: Never compact (open sets are not compact)
/// - **Connectedness**: Always connected (no gaps)
/// - **Supremum/Infimum**: `sup = b` and `inf = a` (exist but not achieved)
///
/// ## Use Cases
///
/// Open intervals are ideal when:
/// - You need to exclude boundary values (e.g., `(0, 1)` for strict inequalities)
/// - Working with limit domains where endpoints represent asymptotic behavior
/// - Defining convergence regions where the limit points are approached but not reached
/// - Modeling strict constraints in optimization problems
/// - Representing neighborhoods in topological analysis
///
/// ## Comparison with Other Interval Types
///
/// | Property | Open `(a,b)` | Closed `[a,b]` | Half-open `[a,b)` |
/// |----------|--------------|----------------|-------------------|
/// | Contains `a` | ✗ | ✓ | ✓ |
/// | Contains `b` | ✗ | ✓ | ✗ |
/// | Length | `b - a` | `b - a` | `b - a` |
/// | Compactness | Never | Always | Never |
/// | Supremum | `b` (not achieved) | `b` (achieved) | `b` (not achieved) |
/// | Infimum | `a` (not achieved) | `a` (achieved) | `a` (achieved) |
/// | Closure | `[a, b]` | `[a, b]` | `[a, b]` |
/// | Interior | `(a, b)` | `(a, b)` | `(a, b)` |
///
/// ## Construction
///
/// ### Safe Construction (Recommended)
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// // Use the recommended performance-optimal type
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Fallible construction (returns Result)
/// let result = IntervalOpen::try_new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let interval = result.unwrap();
///
/// // Infallible construction (panics in debug mode if invalid)
/// let interval = IntervalOpen::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// ```
///
/// ### Raw Construction (Use with Caution)
/// ```rust
/// use grid1d::intervals::*;
///
/// // Direct construction with raw f64 - only when you trust your data
/// let interval = IntervalOpen::new(0.0, 1.0);
/// ```
///
/// ## Key Operations
///
/// ### Point Containment (Excludes Endpoints)
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalOpen::new(-1.0, 1.0);
///
/// // Boundary points are excluded
/// assert!(!interval.contains_point(&-1.0));  // Lower bound excluded
/// assert!(!interval.contains_point(&1.0));   // Upper bound excluded
///
/// // Interior points are included
/// assert!(interval.contains_point(&0.0));    // Interior point
/// assert!(interval.contains_point(&-0.5));   // Interior point
/// assert!(interval.contains_point(&0.999));  // Close to boundary but inside
///
/// // Points outside the interval
/// assert!(!interval.contains_point(&-2.0));  // Below interval
/// assert!(!interval.contains_point(&2.0));   // Above interval
/// ```
///
/// ### Geometric Properties
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalOpen::new(-2.0, 3.0);
///
/// // Length is the same as for closed intervals
/// assert_eq!(interval.length().into_inner(), 5.0);           // |b - a| = |3 - (-2)| = 5
/// assert_eq!(interval.midpoint(), 0.5);        // (a + b) / 2 = (-2 + 3) / 2 = 0.5
/// assert!(!interval.is_symmetric());            // a + b ≠ 0
///
/// // Bounds access
/// assert_eq!(interval.lower_bound_value(), &-2.0);
/// assert_eq!(interval.upper_bound_value(), &3.0);
///
/// // Consume to get bounds
/// let (lower, upper) = interval.into_bounds_pair();
/// assert_eq!((lower, upper), (-2.0, 3.0));
/// ```
///
/// ### Interval Containment Testing
/// ```rust
/// use grid1d::intervals::*;
///
/// let outer_open = IntervalOpen::new(0.0, 3.0);
/// let inner_open = IntervalOpen::new(0.5, 2.5);
/// let inner_closed = IntervalClosed::new(0.5, 2.5);
/// let boundary_closed = IntervalClosed::new(0.0, 3.0);
///
/// // Open intervals can contain smaller open intervals
/// assert!(outer_open.contains_interval(&inner_open));     // (0,3) ⊇ (0.5,2.5) ✓
///
/// // Open intervals can contain closed intervals (if boundaries don't touch)
/// assert!(outer_open.contains_interval(&inner_closed));   // (0,3) ⊇ [0.5,2.5] ✓
///
/// // Open intervals cannot contain intervals that include their excluded boundaries
/// assert!(!outer_open.contains_interval(&boundary_closed)); // (0,3) ⊉ [0,3] ✗
/// assert!(!inner_closed.contains_interval(&outer_open));    // [0.5,2.5] ⊉ (0,3) ✗
/// ```
///
/// ## Intersection Examples
///
/// ### Open with Open
/// ```rust
/// use grid1d::intervals::*;
///
/// let a = IntervalOpen::new(0.0, 3.0);    // (0, 3)
/// let b = IntervalOpen::new(1.0, 4.0);    // (1, 4)
///
/// if let Some(intersection) = a.intersection(&b) {
///     // Result: (1, 3) - intersection preserves open boundaries
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
///             IntervalFinitePositiveLength::Open(open_interval)
///         )) => {
///             assert_eq!(open_interval.lower_bound_value(), &1.0);
///             assert_eq!(open_interval.upper_bound_value(), &3.0);
///             assert_eq!(open_interval.length().into_inner(), 2.0);
///         }
///         _ => unreachable!("Expected open interval"),
///     }
/// }
/// ```
///
/// ### Open with Closed
/// ```rust
/// use grid1d::intervals::*;
///
/// let open = IntervalOpen::new(0.0, 3.0);      // (0, 3)
/// let closed = IntervalClosed::new(1.0, 4.0);  // [1, 4]
///
/// if let Some(intersection) = open.intersection(&closed) {
///     // Result: [1, 3) - takes most restrictive boundary combination
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
///             IntervalFinitePositiveLength::LowerClosedUpperOpen(half_open)
///         )) => {
///             assert_eq!(half_open.lower_bound_value(), &1.0);  // Closed from [1,4]
///             assert_eq!(half_open.upper_bound_value(), &3.0);  // Open from (0,3)
///         }
///         _ => unreachable!("Expected half-open interval"),
///     }
/// }
/// ```
///
/// ### Edge Cases: Touching Boundaries
/// ```rust
/// use grid1d::intervals::*;
///
/// let left_open = IntervalOpen::new(0.0, 2.0);   // (0, 2)
/// let right_open = IntervalOpen::new(2.0, 4.0);  // (2, 4)
///
/// // Open intervals that meet at a point have no intersection
/// assert!(left_open.intersection(&right_open).is_none()); // No overlap at excluded point 2
///
/// let left_open = IntervalOpen::new(0.0, 2.0);     // (0, 2)
/// let right_closed = IntervalClosed::new(2.0, 4.0); // [2, 4]
///
/// // Open and closed intervals meeting at boundary also have no intersection
/// assert!(left_open.intersection(&right_closed).is_none()); // Point 2 excluded from (0,2)
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Generic Programming with Different Scalar Types
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::{RealScalar, RealNative64StrictFinite};
/// use try_create::TryNew;
///
/// // Generic function works with any RealScalar type
/// fn analyze_open_interval<T: RealScalar + std::fmt::Display>(
///     interval: &IntervalOpen<T>
/// ) -> String {
///     format!(
///         "Open interval ({}, {}) with length {} and midpoint {}",
///         interval.lower_bound_value(),
///         interval.upper_bound_value(),
///         interval.length(),
///         interval.midpoint()
///     )
/// }
///
/// // Works with validated types
/// let safe_interval = IntervalOpen::new(
///     RealNative64StrictFinite::try_new(0.0).unwrap(),
///     RealNative64StrictFinite::try_new(1.0).unwrap()
/// );
/// println!("{}", analyze_open_interval(&safe_interval));
///
/// // Also works with raw f64
/// let raw_interval = IntervalOpen::new(0.0, 1.0);
/// println!("{}", analyze_open_interval(&raw_interval));
/// ```
///
/// ### Conversion to Other Interval Types
/// ```rust
/// use grid1d::intervals::*;
///
/// let open = IntervalOpen::new(0.0, 1.0);
///
/// // Convert to category enums
/// let finite_positive: IntervalFinitePositiveLength<f64> = open.clone().into();
/// let finite: IntervalFiniteLength<f64> = open.clone().into();
/// let general: Interval<f64> = open.into();
///
/// // Pattern match on the general enum
/// match general {
///     Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
///         IntervalFinitePositiveLength::Open(recovered_open)
///     )) => {
///         assert_eq!(recovered_open.lower_bound_value(), &0.0);
///         assert_eq!(recovered_open.upper_bound_value(), &1.0);
///     }
///     _ => unreachable!("Should be open interval"),
/// }
/// ```
///
/// ## Mathematical Applications
///
/// ### Limit Analysis and Convergence
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Model a convergence neighborhood around a limit point
/// struct ConvergenceNeighborhood {
///     /// Punctured neighborhood: points near the limit but excluding the limit itself
///     neighborhood: IntervalOpen<f64>,
///     limit_point: f64,
/// }
///
/// impl ConvergenceNeighborhood {
///     fn new(limit_point: f64, epsilon: f64) -> Self {
///         Self {
///             neighborhood: IntervalOpen::new(
///                 limit_point - epsilon,
///                 limit_point + epsilon
///             ),
///             limit_point,
///         }
///     }
///     
///     /// Check if a point is in the punctured neighborhood
///     fn is_in_neighborhood(&self, x: f64) -> bool {
///         self.neighborhood.contains_point(&x) && x != self.limit_point
///     }
///     
///     /// Get the epsilon (half-width) of the neighborhood
///     fn epsilon(&self) -> f64 {
///         self.neighborhood.length().into_inner() / 2.0
///     }
/// }
///
/// let neighborhood = ConvergenceNeighborhood::new(2.0, 0.1);
///
/// assert!(neighborhood.is_in_neighborhood(1.95));  // Close to limit
/// assert!(neighborhood.is_in_neighborhood(2.05));  // Close to limit
/// assert!(!neighborhood.is_in_neighborhood(2.0));  // Limit point excluded
/// assert!(!neighborhood.is_in_neighborhood(1.8));  // Outside neighborhood
/// ```
///
/// ### Optimization with Strict Constraints
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Optimization problem with strict inequality constraints
/// struct StrictConstraint {
///     feasible_region: IntervalOpen<f64>
/// }
///
/// impl StrictConstraint {
///     /// Create constraint: a < x < b
///     fn new_strict_bounds(lower: f64, upper: f64) -> Self {
///         Self {
///             feasible_region: IntervalOpen::new(lower, upper)
///         }
///     }
///     
///     /// Check if a point satisfies the strict constraint
///     fn is_feasible(&self, x: f64) -> bool {
///         self.feasible_region.contains_point(&x)
///     }
///     
///     /// Find interior point (useful for interior-point methods)
///     fn interior_point(&self) -> f64 {
///         self.feasible_region.midpoint()
///     }
///     
///     /// Get constraint violation (positive = infeasible)
///     fn violation(&self, x: f64) -> f64 {
///         let lower = *self.feasible_region.lower_bound_value();
///         let upper = *self.feasible_region.upper_bound_value();
///         
///         if x <= lower {
///             lower - x + f64::EPSILON  // Add small epsilon for strict inequality
///         } else if x >= upper {
///             x - upper + f64::EPSILON  // Add small epsilon for strict inequality
///         } else {
///             0.0  // Feasible
///         }
///     }
/// }
///
/// let constraint = StrictConstraint::new_strict_bounds(0.0, 1.0);  // 0 < x < 1
///
/// assert!(constraint.is_feasible(0.5));     // Interior point
/// assert!(!constraint.is_feasible(0.0));    // Boundary excluded
/// assert!(!constraint.is_feasible(1.0));    // Boundary excluded
/// assert!(!constraint.is_feasible(-0.1));   // Outside region
///
/// assert_eq!(constraint.interior_point(), 0.5);
/// assert!(constraint.violation(0.5) == 0.0);
/// assert!(constraint.violation(0.0) > 0.0);
/// assert!(constraint.violation(1.0) > 0.0);
/// ```
///
/// ### Numerical Integration (Special Handling Required)
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Integrate over open interval using transformation
/// fn integrate_over_open_interval(
///     f: impl Fn(f64) -> f64,
///     domain: &IntervalOpen<f64>,
///     num_points: usize
/// ) -> f64 {
///     let a = *domain.lower_bound_value();
///     let b = *domain.upper_bound_value();
///     
///     // Use Gauss-Legendre quadrature which naturally avoids endpoints
///     // Transform from [-1,1] to (a,b)
///     let transform = |t: f64| (a + b) / 2.0 + (b - a) / 2.0 * t;
///     let jacobian = (b - a) / 2.0;
///     
///     // Gauss-Legendre nodes and weights (simplified for 2 points)
///     let nodes = [-1.0 / f64::sqrt(3.0), 1.0 / f64::sqrt(3.0)];
///     let weights = [1.0, 1.0];
///     
///     let mut integral = 0.0;
///     for i in 0..nodes.len() {
///         let x = transform(nodes[i]);
///         integral += weights[i] * f(x);
///     }
///     
///     integral * jacobian
/// }
///
/// // Example: integrate x^2 over (0, 1)
/// let domain = IntervalOpen::new(0.0, 1.0);
/// let integral = integrate_over_open_interval(|x| x * x, &domain, 2);
///
/// // Should approximate ∫₀¹ x² dx = 1/3 ≈ 0.333...
/// // (Note: for open interval, the result is the same as for closed interval
/// //  since the boundary has measure zero)
/// assert!((integral - 1.0/3.0).abs() < 0.01);
/// ```
///
/// ## Performance Characteristics
///
/// ### Memory Layout
/// - **Storage**: Two scalar values plus boundary type information
/// - **Size**: Identical to [`IntervalClosed`] (typically 16-24 bytes for f64)
/// - **Alignment**: Natural alignment of the scalar type
/// - **No heap allocation**: All data stored inline
///
/// ### Operation Complexity
/// | Operation | Time Complexity | Notes |
/// |-----------|-----------------|-------|
/// | **Point queries** | O(1) | Two strict comparisons (`a < x < b`) |
/// | **Containment tests** | O(1) | Boundary logic comparison |
/// | **Intersections** | O(1) | Pure computation, may allocate result |
/// | **Geometric calculations** | O(1) | Length, midpoint, symmetry |
/// | **Construction** | O(1) | Single validation check in debug mode |
///
/// ### Comparison with Closed Intervals
/// ```rust
/// use grid1d::intervals::*;
///
/// // Both have identical computational complexity
/// let open = IntervalOpen::new(0.0, 1.0);
/// let closed = IntervalClosed::new(0.0, 1.0);
///
/// // Same geometric properties
/// assert_eq!(open.length(), closed.length());
/// assert_eq!(open.midpoint(), closed.midpoint());
/// assert_eq!(open.is_symmetric(), closed.is_symmetric());
///
/// // Different boundary behavior
/// assert!(!open.contains_point(&0.0));     // Open excludes boundaries
/// assert!(closed.contains_point(&0.0));    // Closed includes boundaries
///
/// // Same performance characteristics
/// let start = std::time::Instant::now();
/// for i in 0..1000000 {
///     let x = i as f64 / 1000000.0;
///     let _ = open.contains_point(&x);
/// }
/// let open_time = start.elapsed();
///
/// let start = std::time::Instant::now();
/// for i in 0..1000000 {
///     let x = i as f64 / 1000000.0;
///     let _ = closed.contains_point(&x);
/// }
/// let closed_time = start.elapsed();
///
/// // Performance should be virtually identical
/// // (The difference is just `<` vs `<=` comparisons)
/// ```
///
/// ## Error Handling and Edge Cases
///
/// ### Construction Validation
/// ```rust,ignore
/// use grid1d::intervals::*;
///
/// // Same validation as other interval types
/// assert!(IntervalOpen::try_new(1.0, 0.0).is_err());  // Invalid bounds
/// assert!(IntervalOpen::try_new(1.0, 1.0).is_err());  // Zero length
/// assert!(IntervalOpen::try_new(f64::NAN, 1.0).is_err()); // NaN bound
///
/// // Valid construction
/// assert!(IntervalOpen::try_new(0.0, 1.0).is_ok());
///
/// #[cfg(debug_assertions)]
/// {
///     // Panics in debug mode with invalid bounds
///     std::panic::catch_unwind(|| {
///         IntervalOpen::new(1.0, 0.0);
///     }).expect_err("Should panic on invalid bounds");
/// }
/// ```
///
/// ### Numerical Precision Considerations
/// ```rust
/// use grid1d::intervals::*;
///
/// // Very small open intervals
/// let tiny = IntervalOpen::new(1.0, 1.0 + f64::EPSILON);
/// let length = tiny.length().into_inner();
/// assert!(length > 0.0);
/// assert!(length < 1e-15);
/// assert_eq!(length, f64::EPSILON);
///
/// // Points very close to boundary
/// let interval = IntervalOpen::new(0.0, 1.0);
/// assert!(!interval.contains_point(&0.0));           // Exact boundary
/// assert!(interval.contains_point(&f64::EPSILON));   // Just inside
/// assert!(!interval.contains_point(&1.0));           // Exact boundary
/// assert!(interval.contains_point(&(1.0 - f64::EPSILON))); // Just inside
///
/// // Midpoint is always strictly inside
/// let midpoint = interval.midpoint();
/// assert!(interval.contains_point(&midpoint));
/// ```
///
/// ### Common Pitfalls and Solutions
/// ```rust
/// use grid1d::intervals::*;
///
/// // WRONG: Assuming open intervals behave like closed intervals at boundaries
/// fn naive_boundary_check(interval: &IntervalOpen<f64>) -> bool {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     // This is incorrect for open intervals!
///     interval.contains_point(&lower) && interval.contains_point(&upper)
/// }
///
/// // CORRECT: Understanding that open intervals exclude boundaries
/// fn correct_boundary_check(interval: &IntervalOpen<f64>) -> bool {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     // Open intervals never contain their boundaries
///     !interval.contains_point(&lower) && !interval.contains_point(&upper)
/// }
///
/// let open = IntervalOpen::new(0.0, 1.0);
/// assert!(!naive_boundary_check(&open));     // Always false for open intervals
/// assert!(correct_boundary_check(&open));    // Always true for open intervals
///
/// // WRONG: Trying to find points at the exact boundary
/// fn naive_sample_endpoints(interval: &IntervalOpen<f64>) -> Vec<f64> {
///     vec![
///         *interval.lower_bound_value(),  // Not in the interval!
///         *interval.upper_bound_value(),  // Not in the interval!
///     ]
/// }
///
/// // CORRECT: Sample points strictly inside the interval
/// fn correct_sample_interior(interval: &IntervalOpen<f64>, count: usize) -> Vec<f64> {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     let length = interval.length().into_inner();
///     
///     (0..count)
///         .map(|i| {
///             let t = (i + 1) as f64 / (count + 1) as f64;  // Avoid 0 and 1
///             lower + t * length
///         })
///         .collect()
/// }
///
/// let samples = correct_sample_interior(&open, 3);
/// for point in &samples {
///     assert!(open.contains_point(point));  // All points are valid
/// }
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`IntervalOpen`]
/// - **Strict inequalities**: When boundary points must be excluded
/// - **Limit analysis**: For punctured neighborhoods and convergence studies
/// - **Optimization**: Interior-point methods and strict constraints
/// - **Topology**: When you need open sets for topological properties
/// - **Numerical analysis**: When boundary behavior is singular or undefined
///
/// ### When NOT to Use [`IntervalOpen`]
/// - **Integration bounds**: Usually want closed intervals for definite integrals
/// - **Physical measurements**: Real-world ranges typically include boundaries
/// - **Probability ranges**: Cumulative distributions usually include endpoints
/// - **Array indexing**: Discrete domains need inclusive bounds
///
/// ### Recommended Patterns
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// // Define type aliases for your application
/// type Real = RealNative64StrictFiniteInDebug;
/// type OpenInterval = IntervalOpen<Real>;
///
/// // Pattern: Always check containment instead of assuming boundary behavior
/// fn safe_evaluate_over_interval<F>(
///     f: F,
///     domain: &OpenInterval,
///     x: Real
/// ) -> Option<Real>
/// where
///     F: Fn(Real) -> Real,
/// {
///     if domain.contains_point(&x) {
///         Some(f(x))
///     } else {
///         None  // Point not in domain
///     }
/// }
///
/// // Pattern: Use midpoint for guaranteed interior point
/// fn get_safe_interior_point(interval: &OpenInterval) -> Real {
///     interval.midpoint()  // Always strictly inside for positive-length intervals
/// }
///
/// // Pattern: Convert to closed interval when you need inclusive boundaries
/// fn extend_to_closure(interval: OpenInterval) -> IntervalClosed<Real> {
///     let (lower, upper) = interval.into_bounds_pair();
///     IntervalClosed::new(lower, upper)
/// }
/// ```
///
/// ## Mathematical Correctness Guarantees
///
/// The [`IntervalOpen`] struct maintains all essential mathematical properties:
///
/// - **Open set property**: Never contains its boundary points
/// - **Connectedness**: Always connected (no gaps) for positive-length intervals  
/// - **Length invariance**: `length() = upper - lower` regardless of boundary type
/// - **Intersection correctness**: Results follow standard interval arithmetic
/// - **Containment transitivity**: If A ⊆ B and B ⊆ C, then A ⊆ C
/// - **Supremum/Infimum**: Upper and lower bounds are supremum and infimum respectively
/// - **Interior**: The interval equals its own interior
/// - **Closure**: The closure is the corresponding closed interval `[a, b]`
///
/// These properties are preserved regardless of the scalar type used and are guaranteed
/// by the underlying implementations of the [`IntervalTrait`] and [`IntervalFinitePositiveLengthTrait`].
pub type IntervalOpen<RealType> = IntervalBounded<RealType, Open, Open>;
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
///  A left half-open interval `(a, b]` where the lower endpoint is excluded and the upper endpoint is included.
///
/// The [`IntervalLowerOpenUpperClosed`] represents intervals of the form `(a, b]` where the lower bound `a`
/// is not part of the interval but the upper bound `b` is included. This type of interval is fundamental
/// in mathematical analysis, particularly for right-continuous functions, cumulative distribution functions,
/// and mathematical constructions where you need to exclude the starting point but include the ending point.
///
/// ## Mathematical Properties
///
/// - **Notation**: `(a, b]` or `]a, b]`
/// - **Definition**: `{x ∈ ℝ : a < x ≤ b}` where `a < b`
/// - **Lower boundary**: Excluded (open)
/// - **Upper boundary**: Included (closed)
/// - **Measure**: `|b - a|` (positive length)
/// - **Topology**: Neither open nor closed (half-open set)
/// - **Cardinality**: Uncountably infinite points
/// - **Compactness**: Never compact (not closed)
///
/// ## Use Cases
///
/// Left half-open intervals are ideal when:
/// - **Right-continuous functions**: Domains where functions are continuous from the right
/// - **Cumulative distributions**: CDF values where the upper bound represents an inclusive threshold
/// - **Interval unions**: Creating continuous coverage without overlaps when combined with `[a,b)` intervals
/// - **Limit processes**: Where the target value is achieved but the starting point represents a limit
/// - **Probability theory**: Events where the maximum value is included but minimum is approached
/// - **Numerical analysis**: Discretization schemes requiring upper boundary inclusion
///
/// ## Comparison with Other Interval Types
///
/// | Property | `(a,b]` | `[a,b]` | `(a,b)` | `[a,b)` |
/// |----------|---------|---------|---------|---------|
/// | Contains `a` | ✗ | ✓ | ✗ | ✓ |
/// | Contains `b` | ✓ | ✓ | ✗ | ✗ |
/// | Length | `b - a` | `b - a` | `b - a` | `b - a` |
/// | Compactness | Never | Always | Never | Never |
/// | Supremum | `b` (achieved) | `b` (achieved) | `b` (not achieved) | `b` (not achieved) |
/// | Infimum | `a` (not achieved) | `a` (achieved) | `a` (not achieved) | `a` (achieved) |
/// | Closure | `[a, b]` | `[a, b]` | `[a, b]` | `[a, b]` |
/// | Interior | `(a, b)` | `(a, b)` | `(a, b)` | `(a, b)` |
///
/// ## Construction
///
/// ### Safe Construction (Recommended)
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// // Use the recommended performance-optimal type
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Fallible construction (returns Result)
/// let result = IntervalLowerOpenUpperClosed::try_new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let interval = result.unwrap();
///
/// // Infallible construction (panics in debug mode if invalid)
/// let interval = IntervalLowerOpenUpperClosed::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// ```
///
/// ### Raw Construction (Use with Caution)
/// ```rust
/// use grid1d::intervals::*;
///
/// // Direct construction with raw f64 - only when you trust your data
/// let interval = IntervalLowerOpenUpperClosed::new(0.0, 1.0);
/// ```
///
/// ## Key Operations
///
/// ### Point Containment (Lower Excluded, Upper Included)
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalLowerOpenUpperClosed::new(-1.0, 1.0);  // (-1, 1]
///
/// // Lower boundary is excluded
/// assert!(!interval.contains_point(&-1.0));  // Lower bound excluded
///
/// // Upper boundary is included
/// assert!(interval.contains_point(&1.0));    // Upper bound included
///
/// // Interior points are included
/// assert!(interval.contains_point(&0.0));    // Interior point
/// assert!(interval.contains_point(&0.5));    // Interior point
/// assert!(interval.contains_point(&0.999));  // Close to upper boundary but inside
///
/// // Points outside the interval
/// assert!(!interval.contains_point(&-2.0));  // Below interval
/// assert!(!interval.contains_point(&2.0));   // Above interval
/// ```
///
/// ### Geometric Properties
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalLowerOpenUpperClosed::new(-2.0, 3.0);
///
/// // Length is the same as for other interval types
/// assert_eq!(interval.length().into_inner(), 5.0);           // |b - a| = |3 - (-2)| = 5
/// assert_eq!(interval.midpoint(), 0.5);        // (a + b) / 2 = (-2 + 3) / 2 = 0.5
/// assert!(!interval.is_symmetric());            // a + b ≠ 0
///
/// // Bounds access
/// assert_eq!(interval.lower_bound_value(), &-2.0);
/// assert_eq!(interval.upper_bound_value(), &3.0);
///
/// // Consume to get bounds
/// let (lower, upper) = interval.into_bounds_pair();
/// assert_eq!((lower, upper), (-2.0, 3.0));
/// ```
///
/// ### Boundary Behavior Analysis
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalLowerOpenUpperClosed::new(0.0, 1.0);  // (0, 1]
///
/// // Test boundary inclusion
/// assert!(!interval.contains_point(&0.0));    // Lower: excluded
/// assert!(interval.contains_point(&1.0));     // Upper: included
///
/// // Test points near boundaries
/// assert!(interval.contains_point(&0.0001));  // Just above lower bound
/// assert!(interval.contains_point(&0.9999));  // Just below upper bound
/// assert!(!interval.contains_point(&-0.0001)); // Just below lower bound
/// assert!(!interval.contains_point(&1.0001));  // Just above upper bound
/// ```
///
/// ## Interval Containment and Set Operations
///
/// ### Containment Testing
/// ```rust
/// use grid1d::intervals::*;
///
/// let outer = IntervalLowerOpenUpperClosed::new(0.0, 3.0);  // (0, 3]
/// let inner = IntervalLowerOpenUpperClosed::new(0.5, 2.5);  // (0.5, 2.5]
/// let boundary = IntervalLowerOpenUpperClosed::new(1.0, 3.0); // (1, 3]
///
/// // Left half-open intervals can contain smaller left half-open intervals
/// assert!(outer.contains_interval(&inner));           // (0, 3] ⊇ (0.5, 2.5] ✓
/// assert!(outer.contains_interval(&boundary));        // (0, 3] ⊇ (1, 3] ✓
///
/// // Smaller intervals cannot contain larger ones
/// assert!(!inner.contains_interval(&outer));          // (0.5, 2.5] ⊉ (0, 3] ✗
/// assert!(!boundary.contains_interval(&outer));       // (1, 3] ⊉ (0, 3] ✗
///
/// // Cross-type containment
/// let closed = IntervalClosed::new(0.5, 2.5);               // [0.5, 2.5]
/// let open = IntervalOpen::new(0.5, 2.5);                   // (0.5, 2.5)
/// let other_half = IntervalLowerClosedUpperOpen::new(0.5, 2.5); // [0.5, 2.5)
///
/// assert!(outer.contains_interval(&closed));          // (0, 3] ⊇ [0.5, 2.5] ✓
/// assert!(outer.contains_interval(&open));            // (0, 3] ⊇ (0.5, 2.5) ✓
/// assert!(outer.contains_interval(&other_half));      // (0, 3] ⊇ [0.5, 2.5) ✓
/// ```
///
/// ## Intersection Examples
///
/// ### Left Half-Open with Left Half-Open
/// ```rust
/// use grid1d::intervals::*;
///
/// let a = IntervalLowerOpenUpperClosed::new(0.0, 3.0);  // (0, 3]
/// let b = IntervalLowerOpenUpperClosed::new(1.0, 4.0);  // (1, 4]
///
/// if let Some(intersection) = a.intersection(&b) {
///     // Result: (1, 3] - intersection preserves left half-open boundaries
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
///             IntervalFinitePositiveLength::LowerOpenUpperClosed(half_open)
///         )) => {
///             assert_eq!(half_open.lower_bound_value(), &1.0);
///             assert_eq!(half_open.upper_bound_value(), &3.0);
///             assert_eq!(half_open.length().into_inner(), 2.0);
///         }
///         _ => unreachable!("Expected left half-open interval"),
///     }
/// }
/// ```
///
/// ### Left Half-Open with Other Types
/// ```rust
/// use grid1d::intervals::*;
///
/// let half_open = IntervalLowerOpenUpperClosed::new(0.0, 3.0);  // (0, 3]
/// let closed = IntervalClosed::new(1.0, 4.0);                   // [1, 4]
/// let open = IntervalOpen::new(2.0, 5.0);                       // (2, 5)
/// let other_half = IntervalLowerClosedUpperOpen::new(1.5, 2.5); // [1.5, 2.5)
///
/// // Intersection of (0, 3] with [1, 4]
/// // Result: [1, 3]
/// let half_open_intersection_closed: IntervalClosed<_>
///     = half_open.intersection(&closed).unwrap().try_into().unwrap();
/// assert_eq!(half_open_intersection_closed.lower_bound_value(), &1.0); // Closed from [1, 4]
/// assert_eq!(half_open_intersection_closed.upper_bound_value(), &3.0); // Closed from (0, 3]
///
/// // Intersection of (0, 3] with (2, 5)
/// // Result: (2, 3]
/// let half_open_intersection_open: IntervalLowerOpenUpperClosed<_>
///     = half_open.intersection(&open).unwrap().try_into().unwrap();
/// assert_eq!(half_open_intersection_open.lower_bound_value(), &2.0); // Open from (2, 5)
/// assert_eq!(half_open_intersection_open.upper_bound_value(), &3.0); // Closed from (0, 3]
///
/// // Intersection of (0, 3] with [1.5, 2.5)
/// // Result: [1.5, 2.5)
/// let half_open_intersection_other_half: IntervalLowerClosedUpperOpen<_>
///     = half_open.intersection(&other_half).unwrap().try_into().unwrap();
/// assert_eq!(half_open_intersection_other_half.lower_bound_value(), &1.5); // Closed from [1.5, 2.5)
/// assert_eq!(half_open_intersection_other_half.upper_bound_value(), &2.5); // Open from [1.5, 2.5)
/// ```
///
/// ### Edge Cases: Touching Boundaries
/// ```rust
/// use grid1d::intervals::*;
///
/// let left = IntervalLowerClosedUpperOpen::new(0.0, 2.0);   // [0, 2)
/// let right = IntervalLowerOpenUpperClosed::new(2.0, 4.0);  // (2, 4]
///
/// // Intersection of [0, 2) with (2, 4]
/// // No intersection - gap at point 2
/// assert!(left.intersection(&right).is_none()); // No overlap
///
/// let left_closed = IntervalClosed::new(0.0, 2.0);                  // [0, 2]
/// let right_half = IntervalLowerOpenUpperClosed::new(2.0, 4.0);     // (2, 4]
///
/// // Intersection of [0, 2] with (2, 4]
/// // No intersection - only touching at excluded point 2
/// assert!(left_closed.intersection(&right_half).is_none()); // No overlap at excluded point 2
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Right-Continuous Function Domains
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Represents a right-continuous step function
/// struct RightContinuousFunction {
///     domain: IntervalLowerOpenUpperClosed<f64>,
///     value: f64,
/// }
///
/// impl RightContinuousFunction {
///     fn new(domain_start_exclusive: f64, domain_end_inclusive: f64, function_value: f64) -> Self {
///         Self {
///             domain: IntervalLowerOpenUpperClosed::new(domain_start_exclusive, domain_end_inclusive),
///             value: function_value,
///         }
///     }
///     
///     /// Evaluate the function at a point
///     fn evaluate(&self, x: f64) -> Option<f64> {
///         if self.domain.contains_point(&x) {
///             Some(self.value)
///         } else {
///             None
///         }
///     }
///     
///     /// Check if the function is defined at the right endpoint
///     fn is_right_continuous_at_endpoint(&self) -> bool {
///         let endpoint = *self.domain.upper_bound_value();
///         self.domain.contains_point(&endpoint)  // Should be true for (a,b]
///     }
///     
///     /// Get the domain bounds
///     fn domain_bounds(&self) -> (f64, f64) {
///         (*self.domain.lower_bound_value(), *self.domain.upper_bound_value())
///     }
/// }
///
/// /// Piecewise right-continuous function
/// struct PiecewiseRightContinuous {
///     pieces: Vec<RightContinuousFunction>
/// }
///
/// impl PiecewiseRightContinuous {
///     fn new(pieces: Vec<RightContinuousFunction>) -> Self {
///         Self { pieces }
///     }
///     
///     fn evaluate(&self, x: f64) -> Option<f64> {
///         for piece in &self.pieces {
///             if let Some(value) = piece.evaluate(x) {
///                 return Some(value);
///             }
///         }
///         None
///     }
/// }
///
/// let function = PiecewiseRightContinuous::new(vec![
///     RightContinuousFunction::new(0.0, 1.0, 2.0),    // f(x) = 2 for x ∈ (0, 1]
///     RightContinuousFunction::new(1.0, 2.0, 5.0),    // f(x) = 5 for x ∈ (1, 2]
///     RightContinuousFunction::new(2.0, 3.0, 1.0),    // f(x) = 1 for x ∈ (2, 3]
/// ]);
///
/// assert_eq!(function.evaluate(0.0), None);        // Not in domain
/// assert_eq!(function.evaluate(0.5), Some(2.0));   // In (0, 1]
/// assert_eq!(function.evaluate(1.0), Some(2.0));   // At right endpoint of (0, 1]
/// assert_eq!(function.evaluate(1.5), Some(5.0));   // In (1, 2]
/// assert_eq!(function.evaluate(2.0), Some(5.0));   // At right endpoint of (1, 2]
/// assert_eq!(function.evaluate(3.0), Some(1.0));   // At right endpoint of (2, 3]
/// ```
///
/// ### Numerical Integration with Right Endpoints
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Integrate a right-continuous function over (a,b] using right Riemann sums
/// fn integrate_right_continuous(
///     f: impl Fn(f64) -> f64,
///     domain: IntervalLowerOpenUpperClosed<f64>,
///     num_subdivisions: usize
/// ) -> f64 {
///     let start = *domain.lower_bound_value();
///     let end = *domain.upper_bound_value();
///     let step = domain.length().into_inner() / num_subdivisions as f64;
///     
///     let mut integral = 0.0;
///     
///     for i in 0..num_subdivisions {
///         // Use right endpoint for function evaluation (natural for right-continuous functions)
///         let sub_end = start + (i + 1) as f64 * step;
///         let subinterval_length = step;
///         
///         integral += f(sub_end) * subinterval_length;
///     }
///     
///     integral
/// }
///
/// // Example: integrate f(x) = x over (0, 2]
/// let domain = IntervalLowerOpenUpperClosed::new(0.0, 2.0);
/// let integral = integrate_right_continuous(|x| x, domain, 1000);
///
/// // Should approximate ∫₀² x dx = [x²/2]₀² = 2 (note: endpoint behavior doesn't affect continuous integrals significantly)
/// assert!((integral - 2.0).abs() < 0.01);
/// ```
///
/// ## Performance Characteristics
///
/// ### Memory Layout
/// - **Storage**: Two scalar values plus boundary type information
/// - **Size**: Identical to other bounded interval types (typically 16-24 bytes for f64)
/// - **Alignment**: Natural alignment of the scalar type
/// - **No heap allocation**: All data stored inline
///
/// ### Operation Complexity
/// | Operation | Time Complexity | Notes |
/// |-----------|-----------------|-------|
/// | **Point queries** | O(1) | Two comparisons: `a < x ≤ b` |
/// | **Containment tests** | O(1) | Boundary logic comparison |
/// | **Intersections** | O(1) | Pure computation, may allocate result |
/// | **Geometric calculations** | O(1) | Length, midpoint, symmetry |
/// | **Construction** | O(1) | Single validation check in debug mode |
///
/// ### Comparison with Other Interval Types
/// ```rust
/// use grid1d::intervals::*;
///
/// // All bounded interval types have identical computational complexity
/// let left_half_open = IntervalLowerOpenUpperClosed::new(0.0, 1.0);   // (0, 1]
/// let closed = IntervalClosed::new(0.0, 1.0);                         // [0, 1]
/// let open = IntervalOpen::new(0.0, 1.0);                             // (0, 1)
/// let right_half_open = IntervalLowerClosedUpperOpen::new(0.0, 1.0);  // [0, 1)
///
/// // Same geometric properties
/// assert_eq!(left_half_open.length(), closed.length());
/// assert_eq!(left_half_open.midpoint(), closed.midpoint());
/// assert_eq!(left_half_open.is_symmetric(), closed.is_symmetric());
///
/// // Different boundary behavior
/// assert!(!left_half_open.contains_point(&0.0));     // Excludes lower bound
/// assert!(left_half_open.contains_point(&1.0));      // Includes upper bound
/// assert!(closed.contains_point(&0.0));              // Includes both bounds
/// assert!(closed.contains_point(&1.0));              
/// assert!(!open.contains_point(&0.0));               // Excludes both bounds
/// assert!(!open.contains_point(&1.0));               
/// assert!(right_half_open.contains_point(&0.0));     // Includes lower bound
/// assert!(!right_half_open.contains_point(&1.0));    // Excludes upper bound
/// ```
///
/// ## Mathematical Correctness Guarantees
///
/// The [`IntervalLowerOpenUpperClosed`] struct maintains all essential mathematical properties:
///
/// - **Boundary semantics**: Lower bound always excluded, upper bound always included
/// - **Positive measure**: `length() = upper - lower > 0` for all valid instances
/// - **Set operations**: Containment and intersection follow standard interval arithmetic
/// - **Right-continuity**: Natural choice for right-continuous function domains
/// - **CDF compatibility**: Perfect fit for cumulative distribution step functions
/// - **Order preservation**: Respects the natural ordering of real numbers
/// - **Translation invariance**: Shifting preserves the half-open structure
///
/// These properties are preserved regardless of the scalar type used and are guaranteed
/// by the underlying implementations of the [`IntervalTrait`] and [`IntervalFinitePositiveLengthTrait`].
pub type IntervalLowerOpenUpperClosed<RealType> = IntervalBounded<RealType, Open, Closed>;
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
///  A right half-open interval `[a, b)` where the lower endpoint is included and the upper endpoint is excluded.
///
/// The [`IntervalLowerClosedUpperOpen`] represents intervals of the form `[a, b)` where the lower bound `a`
/// is part of the interval but the upper bound `b` is not included. This type of interval is fundamental
/// in mathematical analysis, particularly for partitions, left-continuous functions, and discrete-continuous
/// interfaces where you need to include the starting point but exclude the ending point.
///
/// ## Mathematical Properties
///
/// - **Notation**: `[a, b)` or `[a, b[`
/// - **Definition**: `{x ∈ ℝ : a ≤ x < b}` where `a < b`
/// - **Lower boundary**: Included (closed)
/// - **Upper boundary**: Excluded (open)
/// - **Measure**: `|b - a|` (positive length)
/// - **Topology**: Neither open nor closed (half-open set)
/// - **Cardinality**: Uncountably infinite points
/// - **Compactness**: Never compact (not closed)
///
/// ## Use Cases
///
/// Right half-open intervals are ideal when:
/// - **Array indexing**: Representing ranges like `[0, n)` for array indices
/// - **Left-continuous functions**: Domains where functions are continuous from the left
/// - **Interval partitions**: Non-overlapping subdivision of larger intervals
/// - **Time intervals**: Start time included, end time excluded (e.g., `[9:00, 10:00)`)
/// - **Discrete-continuous boundaries**: Transitioning from discrete to continuous domains
/// - **Cumulative processes**: Where the starting point is included but the endpoint represents a limit
///
/// ## Comparison with Other Interval Types
///
/// | Property | `[a,b)` | `[a,b]` | `(a,b)` | `(a,b]` |
/// |----------|---------|---------|---------|---------|
/// | Contains `a` | ✓ | ✓ | ✗ | ✗ |
/// | Contains `b` | ✗ | ✓ | ✗ | ✓ |
/// | Length | `b - a` | `b - a` | `b - a` | `b - a` |
/// | Compactness | Never | Always | Never | Never |
/// | Supremum | `b` (not achieved) | `b` (achieved) | `b` (not achieved) | `b` (achieved) |
/// | Infimum | `a` (achieved) | `a` (achieved) | `a` (not achieved) | `a` (not achieved) |
/// | Closure | `[a, b]` | `[a, b]` | `[a, b]` | `[a, b]` |
/// | Interior | `(a, b)` | `(a, b)` | `(a, b)` | `(a, b)` |
///
/// ## Construction
///
/// ### Safe Construction (Recommended)
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// // Use the recommended performance-optimal type
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Fallible construction (returns Result)
/// let result = IntervalLowerClosedUpperOpen::try_new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let interval = result.unwrap();
///
/// // Infallible construction (panics in debug mode if invalid)
/// let interval = IntervalLowerClosedUpperOpen::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// ```
///
/// ### Raw Construction (Use with Caution)
/// ```rust
/// use grid1d::intervals::*;
///
/// // Direct construction with raw f64 - only when you trust your data
/// let interval = IntervalLowerClosedUpperOpen::new(0.0, 1.0);
/// ```
///
/// ## Key Operations
///
/// ### Point Containment (Lower Included, Upper Excluded)
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalLowerClosedUpperOpen::new(-1.0, 1.0);  // [-1, 1)
///
/// // Lower boundary is included
/// assert!(interval.contains_point(&-1.0));   // Lower bound included
///
/// // Upper boundary is excluded
/// assert!(!interval.contains_point(&1.0));   // Upper bound excluded
///
/// // Interior points are included
/// assert!(interval.contains_point(&0.0));    // Interior point
/// assert!(interval.contains_point(&-0.5));   // Interior point
/// assert!(interval.contains_point(&0.999));  // Close to upper boundary but inside
///
/// // Points outside the interval
/// assert!(!interval.contains_point(&-2.0));  // Below interval
/// assert!(!interval.contains_point(&2.0));   // Above interval
/// ```
///
/// ### Geometric Properties
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalLowerClosedUpperOpen::new(-2.0, 3.0);
///
/// // Length is the same as for other interval types
/// assert_eq!(interval.length().into_inner(), 5.0);           // |b - a| = |3 - (-2)| = 5
/// assert_eq!(interval.midpoint(), 0.5);        // (a + b) / 2 = (-2 + 3) / 2 = 0.5
/// assert!(!interval.is_symmetric());            // a + b ≠ 0
///
/// // Bounds access
/// assert_eq!(interval.lower_bound_value(), &-2.0);
/// assert_eq!(interval.upper_bound_value(), &3.0);
///
/// // Consume to get bounds
/// let (lower, upper) = interval.into_bounds_pair();
/// assert_eq!((lower, upper), (-2.0, 3.0));
/// ```
///
/// ### Boundary Behavior Analysis
/// ```rust
/// use grid1d::intervals::*;
///
/// let interval = IntervalLowerClosedUpperOpen::new(0.0, 1.0);  // [0, 1)
///
/// // Test boundary inclusion
/// assert!(interval.contains_point(&0.0));     // Lower: included
/// assert!(!interval.contains_point(&1.0));    // Upper: excluded
///
/// // Test points near boundaries
/// assert!(interval.contains_point(&0.0001));  // Just above lower bound
/// assert!(interval.contains_point(&0.9999));  // Just below upper bound
/// assert!(!interval.contains_point(&-0.0001)); // Just below lower bound
/// assert!(!interval.contains_point(&1.0001));  // Just above upper bound
/// ```
///
/// ## Interval Containment and Set Operations
///
/// ### Containment Testing
/// ```rust
/// use grid1d::intervals::*;
///
/// let outer = IntervalLowerClosedUpperOpen::new(0.0, 3.0);  // [0, 3)
/// let inner = IntervalLowerClosedUpperOpen::new(0.5, 2.5);  // [0.5, 2.5)
/// let boundary = IntervalLowerClosedUpperOpen::new(0.0, 1.0); // [0, 1)
///
/// // Half-open intervals can contain smaller half-open intervals
/// assert!(outer.contains_interval(&inner));           // [0, 3) ⊇ [0.5, 2.5) ✓
/// assert!(outer.contains_interval(&boundary));        // [0, 3) ⊇ [0, 1) ✓
///
/// // Smaller intervals cannot contain larger ones
/// assert!(!inner.contains_interval(&outer));          // [0.5, 2.5) ⊉ [0, 3) ✗
/// assert!(!boundary.contains_interval(&outer));       // [0, 1) ⊉ [0, 3) ✗
///
/// // Cross-type containment
/// let closed = IntervalClosed::new(0.5, 2.0);           // [0.5, 2.0]
/// let open = IntervalOpen::new(0.5, 2.5);               // (0.5, 2.5)
/// let other_half = IntervalLowerOpenUpperClosed::new(0.5, 2.0); // (0.5, 2.0]
///
/// assert!(outer.contains_interval(&closed));          // [0, 3) ⊇ [0.5, 2.0] ✓
/// assert!(outer.contains_interval(&open));            // [0, 3) ⊇ (0.5, 2.5) ✓
/// assert!(outer.contains_interval(&other_half));      // [0, 3) ⊇ (0.5, 2.0] ✓
/// ```
///
/// ## Intersection Examples
///
/// ### Right Half-Open with Right Half-Open
/// ```rust
/// use grid1d::intervals::*;
///
/// let a = IntervalLowerClosedUpperOpen::new(0.0, 3.0);  // [0, 3)
/// let b = IntervalLowerClosedUpperOpen::new(1.0, 4.0);  // [1, 4)
///
/// if let Some(intersection) = a.intersection(&b) {
///     // Result: [1, 3) - intersection preserves right half-open boundaries
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
///             IntervalFinitePositiveLength::LowerClosedUpperOpen(half_open)
///         )) => {
///             assert_eq!(half_open.lower_bound_value(), &1.0);
///             assert_eq!(half_open.upper_bound_value(), &3.0);
///             assert_eq!(half_open.length().into_inner(), 2.0);
///         }
///         _ => unreachable!("Expected right half-open interval"),
///     }
/// }
/// ```
///
/// ### Right Half-Open with Other Types
/// ```rust
/// use grid1d::intervals::*;
///
/// let half_open = IntervalLowerClosedUpperOpen::new(0.0, 3.0);  // [0., 3.)
/// let closed = IntervalClosed::new(1.0, 4.0);                   // [1., 4.]
/// let open = IntervalOpen::new(2.0, 5.0);                       // (2., 5.)
/// let other_half = IntervalLowerOpenUpperClosed::new(1.5, 2.5); // (1.5, 2.5]
///
/// // Intersection of [0., 3.) with [1., 4.]
/// // Result: [1, 3)
/// let half_open_intersection_closed: IntervalLowerClosedUpperOpen<_>
///     = half_open.intersection(&closed).unwrap().try_into().unwrap();
/// assert_eq!(half_open_intersection_closed.lower_bound_value(), &1.0); // Closed from [1. ,4.]
/// assert_eq!(half_open_intersection_closed.upper_bound_value(), &3.0); // Open from [0., 3.)
///
/// // Intersection of [0., 3.) with (2., 5.)
/// // Result: (2, 3)
/// let half_open_intersection_open: IntervalOpen<_>
///     = half_open.intersection(&open).unwrap().try_into().unwrap();
/// assert_eq!(half_open_intersection_open.lower_bound_value(), &2.0); // Open from (2., 5.)
/// assert_eq!(half_open_intersection_open.upper_bound_value(), &3.0); // Open from [0., 3.]
///
/// // Intersection of [0., 3.) with (1.5, 2.5]
/// // Result: (1.5, 2.5]
/// let half_open_intersection_other_half: IntervalLowerOpenUpperClosed<_>
///     = half_open.intersection(&other_half).unwrap().try_into().unwrap();
/// assert_eq!(half_open_intersection_other_half.lower_bound_value(), &1.5); // Open from (1.5, 2.5]
/// assert_eq!(half_open_intersection_other_half.upper_bound_value(), &2.5); // Closed from (1.5, 2.5]
/// ```
///
/// ### Edge Cases: Touching Boundaries
/// ```rust
/// use grid1d::intervals::*;
///
/// let left = IntervalLowerOpenUpperClosed::new(0.0, 2.0);   // (0., 2.]
/// let right = IntervalLowerClosedUpperOpen::new(2.0, 4.0);  // [2., 4.)
///
/// /// // Intersection of (0., 2.] with [2., 4.]
/// // Result: [2.]
/// let singleton: IntervalSingleton<_>
///     = left.intersection(&right).unwrap().try_into().unwrap();
/// assert_eq!(singleton.value(), &2.0);
///
///
/// let left_open = IntervalOpen::new(0.0, 2.0);                  // (0., 2.)
/// let right_half = IntervalLowerClosedUpperOpen::new(2.0, 4.0); // [2., 4.)
///
/// // Open and right half-open meeting at boundary
/// assert!(left_open.intersection(&right_half).is_none()); // No overlap at excluded point 2
/// ```
///
/// ## Advanced Usage Patterns
///
///
/// ### Time Intervals and Scheduling
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Represents a time interval where start time is included but end time is excluded
/// #[derive(Debug)]
/// struct TimeSlot {
///     interval: IntervalLowerClosedUpperOpen<f64>  // Time in hours
/// }
///
/// impl TimeSlot {
///     fn new(start_hour: f64, end_hour: f64) -> Self {
///         Self {
///             interval: IntervalLowerClosedUpperOpen::new(start_hour, end_hour)
///         }
///     }
///     
///     /// Check if a specific time falls within this slot
///     fn contains_time(&self, time: f64) -> bool {
///         self.interval.contains_point(&time)
///     }
///     
///     /// Get duration in hours
///     fn duration(&self) -> f64 {
///         self.interval.length().into_inner()
///     }
///     
///     /// Check if this slot overlaps with another
///     fn overlaps_with(&self, other: &TimeSlot) -> bool {
///         self.interval.intersection(&other.interval).is_some()
///     }
///     
///     /// Get the overlap duration with another slot
///     fn overlap_duration(&self, other: &TimeSlot) -> f64 {
///         if let Some(intersection) = self.interval.intersection(&other.interval) {
///             match intersection {
///                 Interval::FiniteLength(IntervalFiniteLength::PositiveLength(positive)) => {
///                     positive.length().into_inner()
///                 }
///                 Interval::FiniteLength(IntervalFiniteLength::ZeroLength(_)) => 0.0,
///                 _ => 0.0,
///             }
///         } else {
///             0.0
///         }
///     }
/// }
///
/// let meeting1 = TimeSlot::new(9.0, 10.5);   // 9:00 to 10:30
/// let meeting2 = TimeSlot::new(10.0, 11.0);  // 10:00 to 11:00
/// let meeting3 = TimeSlot::new(10.5, 12.0);  // 10:30 to 12:00
///
/// assert!(meeting1.contains_time(9.0));      // Start time included
/// assert!(meeting1.contains_time(10.25));    // Time during meeting
/// assert!(!meeting1.contains_time(10.5));    // End time excluded
///
/// assert_eq!(meeting1.duration(), 1.5);      // 1.5 hours
/// assert_eq!(meeting2.duration(), 1.0);      // 1 hour
///
/// assert!(meeting1.overlaps_with(&meeting2)); // 9:00-10:30 overlaps 10:00-11:00
/// assert!(!meeting1.overlaps_with(&meeting3)); // 9:00-10:30 doesn't overlap 10:30-12:00 (boundary)
/// assert_eq!(meeting1.overlap_duration(&meeting2), 0.5); // 30 minutes overlap
/// ```
///
/// ### Interval Partitions and Subdivisions
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Create a partition of an interval into right half-open subintervals
/// fn create_partition(
///     start: f64,
///     end: f64,
///     num_parts: usize
/// ) -> Vec<IntervalLowerClosedUpperOpen<f64>> {
///     if num_parts == 0 {
///         return Vec::new();
///     }
///     
///     let step = (end - start) / num_parts as f64;
///     (0..num_parts)
///         .map(|i| {
///             let part_start = start + i as f64 * step;
///             let part_end = start + (i + 1) as f64 * step;
///             IntervalLowerClosedUpperOpen::new(part_start, part_end)
///         })
///         .collect()
/// }
///
/// /// Verify that a partition covers the original interval exactly
/// fn verify_partition_coverage(
///     original_start: f64,
///     original_end: f64,
///     partition: &[IntervalLowerClosedUpperOpen<f64>]
/// ) -> bool {
///     if partition.is_empty() {
///         return false;
///     }
///     
///     // Check that first interval starts at original start
///     if *partition[0].lower_bound_value() != original_start {
///         return false;
///     }
///     
///     // Check that last interval ends at original end
///     if *partition.last().unwrap().upper_bound_value() != original_end {
///         return false;
///     }
///     
///     // Check continuity: each interval starts where previous ends
///     for i in 1..partition.len() {
///         if *partition[i].lower_bound_value() != *partition[i-1].upper_bound_value() {
///             return false;
///         }
///     }
///     
///     // Check no overlaps (except at boundaries)
///     for i in 0..partition.len()-1 {
///         for j in i+1..partition.len() {
///             if let Some(intersection) = partition[i].intersection(&partition[j]) {
///                 // Only allow zero-length (singleton) intersections at boundaries
///                 match intersection {
///                     Interval::FiniteLength(IntervalFiniteLength::ZeroLength(_)) => continue,
///                     _ => return false,
///                 }
///             }
///         }
///     }
///     
///     true
/// }
///
/// let partition = create_partition(0.0, 10.0, 4);
/// assert_eq!(partition.len(), 4);
///
/// // Verify each part: [0,2.5), [2.5,5), [5,7.5), [7.5,10)
/// assert_eq!(partition[0].lower_bound_value(), &0.0);
/// assert_eq!(partition[0].upper_bound_value(), &2.5);
/// assert_eq!(partition[1].lower_bound_value(), &2.5);
/// assert_eq!(partition[1].upper_bound_value(), &5.0);
/// assert_eq!(partition[2].lower_bound_value(), &5.0);
/// assert_eq!(partition[2].upper_bound_value(), &7.5);
/// assert_eq!(partition[3].lower_bound_value(), &7.5);
/// assert_eq!(partition[3].upper_bound_value(), &10.0);
///
/// assert!(verify_partition_coverage(0.0, 10.0, &partition));
///
/// // Test point coverage - every point in [0,10) should be in exactly one partition interval
/// for test_point in [0.0, 1.0, 2.5, 4.0, 5.0, 7.5, 9.0] {
///     let containing_intervals: Vec<_> = partition.iter()
///         .enumerate()
///         .filter(|(_, interval)| interval.contains_point(&test_point))
///         .collect();
///     
///     assert_eq!(containing_intervals.len(), 1,
///                "Point {} should be in exactly one interval", test_point);
/// }
/// ```
///
/// ### Numerical Integration with Uniform Partitions
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Integrate a function over [a,b) using rectangular rule with right half-open subintervals
/// fn integrate_with_right_half_open_partition(
///     f: impl Fn(f64) -> f64,
///     domain: IntervalLowerClosedUpperOpen<f64>,
///     num_subdivisions: usize
/// ) -> f64 {
///     let start = *domain.lower_bound_value();
///     let end = *domain.upper_bound_value();
///     let step = domain.length().into_inner() / num_subdivisions as f64;
///     
///     let mut integral = 0.0;
///     
///     for i in 0..num_subdivisions {
///         let sub_start = start + i as f64 * step;
///         let sub_end = start + (i + 1) as f64 * step;
///         let subinterval = IntervalLowerClosedUpperOpen::new(sub_start, sub_end);
///         
///         // Use left endpoint for function evaluation (since it's included in [a,b))
///         let sample_point = *subinterval.lower_bound_value();
///         integral += f(sample_point) * subinterval.length().into_inner();
///     }
///     
///     integral
/// }
///
/// // Example: integrate f(x) = x² over [0, 3)
/// let domain = IntervalLowerClosedUpperOpen::new(0.0, 3.0);
/// let integral = integrate_with_right_half_open_partition(|x| x * x, domain, 1000);
///
/// // Should approximate ∫₀³ x² dx = [x³/3]₀³ = 9 (note: upper bound excluded doesn't affect continuous integrals significantly)
/// assert!((integral - 9.0).abs() < 0.1);
/// ```
///
/// ## Performance Characteristics
///
/// ### Memory Layout
/// - **Storage**: Two scalar values plus boundary type information
/// - **Size**: Identical to other bounded interval types (typically 16-24 bytes for f64)
/// - **Alignment**: Natural alignment of the scalar type
/// - **No heap allocation**: All data stored inline
///
/// ### Operation Complexity
/// | Operation | Time Complexity | Notes |
/// |-----------|-----------------|-------|
/// | **Point queries** | O(1) | Two comparisons: `a ≤ x < b` |
/// | **Containment tests** | O(1) | Boundary logic comparison |
/// | **Intersections** | O(1) | Pure computation, may allocate result |
/// | **Geometric calculations** | O(1) | Length, midpoint, symmetry |
/// | **Construction** | O(1) | Single validation check in debug mode |
///
/// ### Comparison with Other Interval Types
/// ```rust
/// use grid1d::intervals::*;
///
/// // All bounded interval types have identical computational complexity
/// let right_half_open = IntervalLowerClosedUpperOpen::new(0.0, 1.0);  // [0, 1)
/// let closed = IntervalClosed::new(0.0, 1.0);                         // [0, 1]
/// let open = IntervalOpen::new(0.0, 1.0);                             // (0, 1)
/// let left_half_open = IntervalLowerOpenUpperClosed::new(0.0, 1.0);   // (0, 1]
///
/// // Same geometric properties
/// assert_eq!(right_half_open.length(), closed.length());
/// assert_eq!(right_half_open.midpoint(), closed.midpoint());
/// assert_eq!(right_half_open.is_symmetric(), closed.is_symmetric());
///
/// // Different boundary behavior
/// assert!(right_half_open.contains_point(&0.0));     // Includes lower bound
/// assert!(!right_half_open.contains_point(&1.0));    // Excludes upper bound
/// assert!(closed.contains_point(&0.0));              // Includes both bounds
/// assert!(closed.contains_point(&1.0));              
/// assert!(!open.contains_point(&0.0));               // Excludes both bounds
/// assert!(!open.contains_point(&1.0));               
/// assert!(!left_half_open.contains_point(&0.0));     // Excludes lower bound
/// assert!(left_half_open.contains_point(&1.0));      // Includes upper bound
/// ```
///
/// ## Error Handling and Edge Cases
///
/// ### Construction Validation
/// ```rust,ignore
/// use grid1d::intervals::*;
///
/// // Valid construction
/// assert!(IntervalLowerClosedUpperOpen::try_new(0.0, 1.0).is_ok());
///
/// // Invalid construction
/// assert!(IntervalLowerClosedUpperOpen::try_new(1.0, 0.0).is_err());  // Inverted bounds
/// assert!(IntervalLowerClosedUpperOpen::try_new(1.0, 1.0).is_err());  // Zero length
/// assert!(IntervalLowerClosedUpperOpen::try_new(f64::NAN, 1.0).is_err()); // NaN bound
///
/// #[cfg(debug_assertions)]
/// {
///     // Panics in debug mode with invalid bounds
///     std::panic::catch_unwind(|| {
///         IntervalLowerClosedUpperOpen::new(1.0, 0.0);
///     }).expect_err("Should panic on invalid bounds");
/// }
/// ```
///
/// ### Boundary Precision Considerations
/// ```rust
/// use grid1d::intervals::*;
///
/// // Very small intervals near floating-point precision limits
/// let tiny_interval = IntervalLowerClosedUpperOpen::new(1.0, 1.0 + f64::EPSILON);
/// assert!(tiny_interval.length().into_inner() > 0.0);
/// assert!(tiny_interval.length().into_inner() < 1e-15);
///
/// // Points very close to boundaries
/// let interval = IntervalLowerClosedUpperOpen::new(0.0, 1.0);
/// assert!(interval.contains_point(&0.0));                    // Exact lower boundary
/// assert!(interval.contains_point(&f64::EPSILON));           // Just above lower bound
/// assert!(!interval.contains_point(&1.0));                   // Exact upper boundary
/// assert!(interval.contains_point(&(1.0 - f64::EPSILON)));   // Just below upper bound
/// ```
///
/// ### Common Pitfalls and Solutions
/// ```rust
/// use grid1d::intervals::*;
///
/// // WRONG: Assuming both boundaries are included
/// fn naive_boundary_check(interval: &IntervalLowerClosedUpperOpen<f64>) -> (bool, bool) {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     // This gives wrong results for the upper bound!
///     (
///         interval.contains_point(&lower),  // Correct: returns true
///         interval.contains_point(&upper)   // Wrong assumption: returns false
///     )
/// }
///
/// // CORRECT: Understanding the boundary semantics
/// fn correct_boundary_check(interval: &IntervalLowerClosedUpperOpen<f64>) -> (bool, bool) {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     (
///         interval.contains_point(&lower),  // Should be true (closed lower bound)
///         interval.contains_point(&upper)   // Should be false (open upper bound)
///     )
/// }
///
/// let interval = IntervalLowerClosedUpperOpen::new(0.0, 1.0);
/// assert_eq!(naive_boundary_check(&interval), (true, false));   // Correct by accident
/// assert_eq!(correct_boundary_check(&interval), (true, false)); // Correct by design
///
/// // WRONG: Creating ranges that don't partition correctly
/// fn naive_create_adjacent_ranges(start: f64, middle: f64, end: f64) ->
///     (IntervalLowerClosedUpperOpen<f64>, IntervalLowerClosedUpperOpen<f64>) {
///     // This creates a gap at the middle point!
///     (
///         IntervalLowerClosedUpperOpen::new(start, middle),    // [start, middle)
///         IntervalLowerClosedUpperOpen::new(middle + 0.001, end) // [middle+ε, end) - GAP!
///     )
/// }
///
/// // CORRECT: Ensuring proper adjacency
/// fn correct_create_adjacent_ranges(start: f64, middle: f64, end: f64) ->
///     (IntervalLowerClosedUpperOpen<f64>, IntervalLowerClosedUpperOpen<f64>) {
///     (
///         IntervalLowerClosedUpperOpen::new(start, middle),  // [start, middle)
///         IntervalLowerClosedUpperOpen::new(middle, end)     // [middle, end) - NO GAP
///     )
/// }
///
/// let (left, right) = correct_create_adjacent_ranges(0.0, 5.0, 10.0);
/// // Verify no gap: right interval starts exactly where left interval ends
/// assert_eq!(left.upper_bound_value(), right.lower_bound_value());
/// // Verify no overlap: intersection should be at most a single point
/// if let Some(intersection) = left.intersection(&right) {
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(_)) => {
///             // This is fine - single point intersection at boundary
///         }
///         _ => panic!("Should not have positive-length intersection"),
///     }
/// }
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`IntervalLowerClosedUpperOpen`]
/// - **Array indexing**: Representing index ranges like `[start, end)`
/// - **Time intervals**: Start time included, end time excluded
/// - **Partitions**: Creating non-overlapping subdivisions
/// - **Left-continuous functions**: Domains where functions are continuous from the left
/// - **Iteration ranges**: Programming constructs like `for i in start..end`
/// - **Histogram bins**: Where lower bound is included, upper bound excluded
///
/// ### When NOT to Use [`IntervalLowerClosedUpperOpen`]
/// - **Mathematical integration**: Usually want closed intervals `[a,b]` for definite integrals
/// - **Physical measurements**: Real-world ranges typically include both boundaries
/// - **Probability calculations**: Often need closed intervals for cumulative distributions
/// - **Optimization domains**: Constraint boundaries often need to be included
///
/// ### Recommended Patterns
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::{RealNative64StrictFiniteInDebug, RealScalar};
/// use try_create::TryNew;
///
/// // Define type aliases for your application
/// type Real = RealNative64StrictFiniteInDebug;
/// type RightHalfOpenInterval = IntervalLowerClosedUpperOpen<Real>;
///
/// // Pattern: Use for array-like indexing
/// fn create_index_range(start: usize, count: usize) -> RightHalfOpenInterval {
///     let end = start + count;
///     RightHalfOpenInterval::new(
///         Real::try_new(start as f64).unwrap(),
///         Real::try_new(end as f64).unwrap()
///     )
/// }
///
/// // Pattern: Verify boundary semantics explicitly
/// fn validate_range_boundaries(range: &RightHalfOpenInterval) -> bool {
///     let lower = range.lower_bound_value();
///     let upper = range.upper_bound_value();
///     
///     // Explicitly test the expected boundary behavior
///     range.contains_point(lower) && !range.contains_point(upper)
/// }
///
/// // Pattern: Create partitions with proper adjacency
/// fn create_uniform_partition(
///     total_range: RightHalfOpenInterval,
///     num_parts: usize
/// ) -> Vec<RightHalfOpenInterval> {
///     let start = total_range.lower_bound_value().clone();
///     let length = total_range.length().into_inner();
///     let part_length = length / Real::try_from_f64(num_parts as f64).unwrap();
///     
///     (0..num_parts)
///         .map(|i| {
///             let part_start = start.clone() + Real::try_from_f64(i as f64).unwrap() * part_length.clone();
///             let part_end = part_start.clone() + part_length.clone();
///             RightHalfOpenInterval::new(part_start, part_end)
///         })
///         .collect()
/// }
/// ```
///
/// ## Mathematical Correctness Guarantees
///
/// The [`IntervalLowerClosedUpperOpen`] struct maintains all essential mathematical properties:
///
/// - **Boundary semantics**: Lower bound always included, upper bound always excluded
/// - **Positive measure**: `length() = upper - lower > 0` for all valid instances
/// - **Set operations**: Containment and intersection follow standard interval arithmetic
/// - **Partition compatibility**: Adjacent right half-open intervals create proper partitions
/// - **Continuity properties**: Suitable for left-continuous function domains
/// - **Order preservation**: Respects the natural ordering of real numbers
/// - **Translation invariance**: Shifting preserves the half-open structure
///
/// These properties are preserved regardless of the scalar type used and are guaranteed
/// by the underlying implementations of the [`IntervalTrait`] and [`IntervalFinitePositiveLengthTrait`].
pub type IntervalLowerClosedUpperOpen<RealType> = IntervalBounded<RealType, Closed, Open>;
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
///  Container for all intervals with finite, positive, measurable length.
///
/// The [`IntervalFinitePositiveLength`] enum represents intervals where the upper bound
/// is strictly greater than the lower bound, resulting in a positive measure and
/// uncountably infinite cardinality. This enum encompasses all bounded intervals
/// with non-zero length, distinguished by their boundary inclusion behavior.
///
/// ## Mathematical Foundation
///
/// All variants represent intervals of the form `{x ∈ ℝ : a < x < b}` or similar,
/// where the specific boundary conditions depend on the variant. The fundamental
/// mathematical properties shared by all variants are:
///
/// - **Positive Measure**: `μ(I) = |b - a| > 0` for all variants
/// - **Non-empty Interior**: All intervals contain infinitely many interior points
/// - **Boundedness**: Both endpoints are finite real numbers
/// - **Connectedness**: All intervals are connected sets with no gaps
/// - **Uncountable Cardinality**: Each interval contains uncountably many points
///
/// ## Generic Over Any [`num_valid::RealScalar`] Type
///
/// **The [`IntervalFinitePositiveLength`] enum works seamlessly with any scalar type implementing [`num_valid::RealScalar`]:**
///
/// | Scalar Type | Performance | Validation | Best For |
/// |-------------|-------------|------------|----------|
/// | [`f64`] | ⚡⚡⚡ Maximum | ❌ None | Trusted input, maximum speed |
/// | [`RealNative64StrictFiniteInDebug`](num_valid::RealNative64StrictFiniteInDebug) | ⚡⚡⚡ **Same as f64** | ✅ **Debug only** | **Recommended for most uses** |
/// | [`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`)|
///
/// ## Enum Variants and Boundary Semantics
///
/// The four variants differ only in their boundary inclusion behavior:
///
/// ### [`IntervalFinitePositiveLength::Closed`] - `[a, b]`
/// **Both endpoints included**
/// - **Mathematical notation**: `[a, b] = {x ∈ ℝ : a ≤ x ≤ b}`
/// - **Boundary behavior**: Contains both `a` and `b`
/// - **Topology**: Compact set (closed and bounded)
/// - **Supremum/Infimum**: Both achieved at the endpoints
/// - **Use cases**: Probability ranges, integration bounds, physical constraints
///
/// ### [`IntervalFinitePositiveLength::Open`] - `(a, b)`
/// **Both endpoints excluded**
/// - **Mathematical notation**: `(a, b) = {x ∈ ℝ : a < x < b}`
/// - **Boundary behavior**: Excludes both `a` and `b`
/// - **Topology**: Open set, not compact
/// - **Supremum/Infimum**: Exist but not achieved
/// - **Use cases**: Strict inequalities, limit domains, convergence regions
///
/// ### [`IntervalFinitePositiveLength::LowerOpenUpperClosed`] - `(a, b]`
/// **Lower endpoint excluded, upper endpoint included**
/// - **Mathematical notation**: `(a, b] = {x ∈ ℝ : a < x ≤ b}`
/// - **Boundary behavior**: Excludes `a`, includes `b`
/// - **Topology**: Neither open nor closed (half-open)
/// - **Use cases**: Right-continuous functions, cumulative distributions
///
/// ### [`IntervalFinitePositiveLength::LowerClosedUpperOpen`] - `[a, b)`
/// **Lower endpoint included, upper endpoint excluded**
/// - **Mathematical notation**: `[a, b) = {x ∈ ℝ : a ≤ x < b}`
/// - **Boundary behavior**: Includes `a`, excludes `b`
/// - **Topology**: Neither open nor closed (half-open)
/// - **Use cases**: Left-continuous functions, array indexing, partitions
///
/// ## Construction Patterns
///
/// ### From Concrete Interval Types (Recommended)
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::TryNew;
///
/// // Use the recommended performance-optimal type
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Create specific interval types, then convert
/// let closed = IntervalClosed::new(
///     Real::try_new(0.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let positive: IntervalFinitePositiveLength<Real> = closed.into();
///
/// let open = IntervalOpen::new(
///     Real::try_new(-1.0).unwrap(),
///     Real::try_new(1.0).unwrap()
/// );
/// let positive: IntervalFinitePositiveLength<Real> = open.into();
/// ```
///
/// ### Direct Construction (Advanced)
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// // Build variants directly
/// let closed = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 1.0)
/// );
///
/// let open = IntervalFinitePositiveLength::Open(
///     IntervalOpen::new(0.0, 1.0)
/// );
///
/// let half_open_left = IntervalFinitePositiveLength::LowerOpenUpperClosed(
///     IntervalLowerOpenUpperClosed::new(0.0, 1.0)
/// );
///
/// let half_open_right = IntervalFinitePositiveLength::LowerClosedUpperOpen(
///     IntervalLowerClosedUpperOpen::new(0.0, 1.0)
/// );
/// ```
///
/// ## Core Operations
///
/// All [`IntervalFinitePositiveLength`] instances support the complete [`IntervalTrait`] interface
/// plus specialized operations for positive-length intervals:
///
/// ### Point Containment Testing
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let closed = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 1.0)
/// );
/// let open = IntervalFinitePositiveLength::Open(
///     IntervalOpen::new(0.0, 1.0)
/// );
///
/// // Closed interval includes endpoints
/// assert!(closed.contains_point(&0.0));    // Lower bound included
/// assert!(closed.contains_point(&0.5));    // Interior point
/// assert!(closed.contains_point(&1.0));    // Upper bound included
/// assert!(!closed.contains_point(&2.0));   // Outside interval
///
/// // Open interval excludes endpoints
/// assert!(!open.contains_point(&0.0));     // Lower bound excluded
/// assert!(open.contains_point(&0.5));      // Interior point
/// assert!(!open.contains_point(&1.0));     // Upper bound excluded
/// assert!(!open.contains_point(&2.0));     // Outside interval
/// ```
///
/// ### Geometric Properties
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let interval = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(-2.0, 3.0)
/// );
///
/// // All variants have identical geometric properties
/// assert_eq!(interval.length().into_inner(), 5.0);           // |b - a| = |3 - (-2)| = 5
/// assert_eq!(interval.midpoint(), 0.5);        // (a + b) / 2 = (-2 + 3) / 2 = 0.5
/// assert!(!interval.is_symmetric());            // a + b ≠ 0
///
/// // Bounds access
/// assert_eq!(interval.lower_bound_value(), &-2.0);
/// assert_eq!(interval.upper_bound_value(), &3.0);
///
/// // Consume to get bounds
/// let (lower, upper) = interval.into_bounds_pair();
/// assert_eq!((lower, upper), (-2.0, 3.0));
/// ```
///
/// ### Interval Containment Testing
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let outer_closed = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 2.0)
/// );
/// let inner_open = IntervalFinitePositiveLength::Open(
///     IntervalOpen::new(0.5, 1.5)
/// );
/// let boundary_half_open = IntervalFinitePositiveLength::LowerClosedUpperOpen(
///     IntervalLowerClosedUpperOpen::new(0.0, 1.0)
/// );
///
/// // Closed intervals can contain open intervals
/// assert!(outer_closed.contains_interval(&inner_open));           // [0,2] ⊇ (0.5,1.5) ✓
/// assert!(outer_closed.contains_interval(&boundary_half_open));   // [0,2] ⊇ [0,1) ✓
///
/// // Open intervals cannot contain intervals that include excluded boundaries
/// assert!(!inner_open.contains_interval(&outer_closed));          // (0.5,1.5) ⊉ [0,2] ✗
/// assert!(!inner_open.contains_interval(&boundary_half_open));    // (0.5,1.5) ⊉ [0,1) ✗
/// ```
///
/// ### Intersection Operations
/// ```rust
/// use grid1d::intervals::*;
///
/// let closed = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 2.0)
/// );
/// let open = IntervalFinitePositiveLength::Open(
///     IntervalOpen::new(1.0, 3.0)
/// );
///
/// if let Some(intersection) = closed.intersection(&open) {
///     // Result: (1.0, 2.0] - most restrictive boundary combination
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(positive)) => {
///             println!("Intersection length: {}", positive.length()); // 1.0
///             
///             // Can extract specific type if needed
///             if let IntervalFinitePositiveLength::LowerOpenUpperClosed(half_open) = positive {
///                 assert_eq!(half_open.lower_bound_value(), &1.0);
///                 assert_eq!(half_open.upper_bound_value(), &2.0);
///             }
///         }
///         _ => unreachable!("Positive-length intervals always have positive-length intersections"),
///     }
/// }
///
/// // Edge case: intersection at single point becomes singleton
/// let touching_left = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 1.0)
/// );
/// let touching_right = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(1.0, 2.0)
/// );
///
/// if let Some(intersection) = touching_left.intersection(&touching_right) {
///     // Result: singleton at point 1.0
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(singleton)) => {
///             assert_eq!(singleton.value(), &1.0);
///         }
///         _ => panic!("Expected singleton intersection"),
///     }
/// }
/// ```
///
/// ## Pattern Matching and Variant Analysis
///
/// ### Distinguishing Between Boundary Types
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// fn analyze_boundary_behavior(interval: &IntervalFinitePositiveLength<f64>) -> String {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     match interval {
///         IntervalFinitePositiveLength::Closed(_) => {
///             format!("[{}, {}] - both endpoints included", lower, upper)
///         }
///         IntervalFinitePositiveLength::Open(_) => {
///             format!("({}, {}) - both endpoints excluded", lower, upper)
///         }
///         IntervalFinitePositiveLength::LowerOpenUpperClosed(_) => {
///             format!("({}, {}] - lower open, upper closed", lower, upper)
///         }
///         IntervalFinitePositiveLength::LowerClosedUpperOpen(_) => {
///             format!("[{}, {}) - lower closed, upper open", lower, upper)
///         }
///     }
/// }
///
/// let closed = IntervalFinitePositiveLength::Closed(IntervalClosed::new(1.0, 2.0));
/// let open = IntervalFinitePositiveLength::Open(IntervalOpen::new(1.0, 2.0));
/// let half_open_left = IntervalFinitePositiveLength::LowerOpenUpperClosed(
///     IntervalLowerOpenUpperClosed::new(1.0, 2.0)
/// );
/// let half_open_right = IntervalFinitePositiveLength::LowerClosedUpperOpen(
///     IntervalLowerClosedUpperOpen::new(1.0, 2.0)
/// );
///
/// println!("{}", analyze_boundary_behavior(&closed));       // "[1, 2] - both endpoints included"
/// println!("{}", analyze_boundary_behavior(&open));         // "(1, 2) - both endpoints excluded"
/// println!("{}", analyze_boundary_behavior(&half_open_left)); // "(1, 2] - lower open, upper closed"
/// println!("{}", analyze_boundary_behavior(&half_open_right)); // "[1, 2) - lower closed, upper open"
/// ```
///
/// ### Boundary Endpoint Testing
/// ```rust
/// use grid1d::intervals::*;
///
/// fn test_boundary_inclusion(interval: &IntervalFinitePositiveLength<f64>) -> (bool, bool) {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     let includes_lower = interval.contains_point(&lower);
///     let includes_upper = interval.contains_point(&upper);
///     
///     (includes_lower, includes_upper)
/// }
///
/// let closed = IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 1.0));
/// let open = IntervalFinitePositiveLength::Open(IntervalOpen::new(0.0, 1.0));
/// let left_open = IntervalFinitePositiveLength::LowerOpenUpperClosed(
///     IntervalLowerOpenUpperClosed::new(0.0, 1.0)
/// );
/// let right_open = IntervalFinitePositiveLength::LowerClosedUpperOpen(
///     IntervalLowerClosedUpperOpen::new(0.0, 1.0)
/// );
///
/// assert_eq!(test_boundary_inclusion(&closed), (true, true));     // [0,1]: both included
/// assert_eq!(test_boundary_inclusion(&open), (false, false));     // (0,1): both excluded
/// assert_eq!(test_boundary_inclusion(&left_open), (false, true)); // (0,1]: lower excluded, upper included
/// assert_eq!(test_boundary_inclusion(&right_open), (true, false)); // [0,1): lower included, upper excluded
/// ```
///
/// ## Mathematical Applications
///
/// ### Integration Domains
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealScalar;
///
/// /// Compute definite integral using midpoint rule (simplified example)
/// fn integrate_midpoint<T: RealScalar + Clone>(
///     f: impl Fn(T) -> T,
///     domain: &IntervalFinitePositiveLength<T>,
///     num_subdivisions: usize
/// ) -> T {
///     let length = domain.length().into_inner();
///     let step = length.clone() / T::try_from_f64(num_subdivisions as f64).unwrap();
///     let mut sum = T::zero();
///     
///     for i in 0..num_subdivisions {
///         let offset = T::try_from_f64(i as f64).unwrap() * step.clone() + step.clone() / T::try_from_f64(2.0).unwrap();
///         let x = domain.lower_bound_value().clone() + offset;
///         sum = sum + f(x);
///     }
///     
///     sum * step
/// }
///
/// // Works with any boundary type - the integral is the same
/// let closed_domain = IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 1.0));
/// let open_domain = IntervalFinitePositiveLength::Open(IntervalOpen::new(0.0, 1.0));
///
/// let f = |x: f64| x * x; // f(x) = x²
///
/// let integral_closed = integrate_midpoint(&f, &closed_domain, 1000);
/// let integral_open = integrate_midpoint(&f, &open_domain, 1000);
///
/// // Both should approximate ∫₀¹ x² dx = 1/3 ≈ 0.333...
/// assert!((integral_closed - 1.0/3.0).abs() < 0.01);
/// assert!((integral_open - 1.0/3.0).abs() < 0.01);
/// ```
///
/// ### Probability Theory and Measure
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Represents a continuous uniform distribution
/// struct UniformDistribution {
///     support: IntervalFinitePositiveLength<f64>
/// }
///
/// impl UniformDistribution {
///     fn new(support: IntervalFinitePositiveLength<f64>) -> Self {
///         Self { support }
///     }
///     
///     /// Probability density function
///     fn pdf(&self, x: f64) -> f64 {
///         if self.support.contains_point(&x) {
///             1.0 / self.support.length().into_inner() // Uniform density
///         } else {
///             0.0 // Outside support
///         }
///     }
///     
///     /// Cumulative distribution function  
///     fn cdf(&self, x: f64) -> f64 {
///         let lower = *self.support.lower_bound_value();
///         let upper = *self.support.upper_bound_value();
///         
///         if x <= lower {
///             0.0
///         } else if x >= upper {
///             1.0
///         } else {
///             (x - lower) / self.support.length().into_inner()
///         }
///     }
///     
///     /// Sample from the distribution (simplified)
///     fn sample(&self) -> f64 {
///         let lower = *self.support.lower_bound_value();
///         let upper = *self.support.upper_bound_value();
///         
///         // In real implementation, use proper random number generator
///         lower + (upper - lower) * 0.5 // Deterministic "sample" for testing
///     }
/// }
///
/// // Works with any interval type
/// let uniform_closed = UniformDistribution::new(
///     IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 2.0))
/// );
/// let uniform_open = UniformDistribution::new(
///     IntervalFinitePositiveLength::Open(IntervalOpen::new(0.0, 2.0))
/// );
///
/// // Same density and behavior regardless of boundary type
/// assert_eq!(uniform_closed.pdf(1.0), 0.5);  // 1 / (2 - 0) = 0.5
/// assert_eq!(uniform_open.pdf(1.0), 0.5);    // Same density
///
/// assert_eq!(uniform_closed.cdf(1.0), 0.5);  // (1 - 0) / (2 - 0) = 0.5
/// assert_eq!(uniform_open.cdf(1.0), 0.5);    // Same CDF for interior points
/// ```
///
/// ### Optimization and Root Finding
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Simple bisection method for root finding
/// fn bisection_method(
///     f: impl Fn(f64) -> f64,
///     mut interval: IntervalFinitePositiveLength<f64>,
///     tolerance: f64
/// ) -> Option<f64> {
///     let mut iterations = 0;
///     const MAX_ITERATIONS: usize = 100;
///     
///     // Check that function changes sign over interval
///     let f_lower = f(*interval.lower_bound_value());
///     let f_upper = f(*interval.upper_bound_value());
///     
///     if f_lower * f_upper > 0.0 {
///         return None; // No sign change, no root guaranteed
///     }
///     
///     while interval.length().into_inner() > tolerance && iterations < MAX_ITERATIONS {
///         let midpoint = interval.midpoint();
///         let f_mid = f(midpoint);
///         
///         if f_mid.abs() < tolerance {
///             return Some(midpoint);
///         }
///         
///         let lower = *interval.lower_bound_value();
///         let upper = *interval.upper_bound_value();
///         
///         // Update interval to half containing the root
///         if f_lower * f_mid < 0.0 {
///             // Root in left half
///             interval = IntervalFinitePositiveLength::Closed(
///                 IntervalClosed::new(lower, midpoint)
///             );
///         } else {
///             // Root in right half  
///             interval = IntervalFinitePositiveLength::Closed(
///                 IntervalClosed::new(midpoint, upper)
///             );
///         }
///         
///         iterations += 1;
///     }
///     
///     Some(interval.midpoint())
/// }
///
/// // Find root of f(x) = x² - 2 (should be √2 ≈ 1.414)
/// let search_interval = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(1.0, 2.0)
/// );
///
/// let root = bisection_method(|x| x * x - 2.0, search_interval, 1e-6).unwrap();
/// assert!((root - 2.0_f64.sqrt()).abs() < 1e-6);
/// ```
///
/// ## Type Conversion and Specialization
///
/// ### Converting to Specific Interval Types
/// ```rust
/// use grid1d::intervals::*;
///
/// let positive = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 1.0)
/// );
///
/// // Convert to general interval enums
/// let finite: IntervalFiniteLength<f64> = positive.clone().into();
/// let general: Interval<f64> = positive.clone().into();
///
/// // Extract specific types when needed
/// match positive {
///     IntervalFinitePositiveLength::Closed(closed_interval) => {
///         println!("Working with closed interval: {:?}", closed_interval);
///         // Access closed-interval-specific behavior if needed
///         let (lower, upper) = closed_interval.into_bounds_pair();
///         assert_eq!((lower, upper), (0.0, 1.0));
///     }
///     IntervalFinitePositiveLength::Open(open_interval) => {
///         println!("Working with open interval: {:?}", open_interval);
///     }
///     IntervalFinitePositiveLength::LowerOpenUpperClosed(half_open) => {
///         println!("Working with left half-open interval: {:?}", half_open);
///     }
///     IntervalFinitePositiveLength::LowerClosedUpperOpen(half_open) => {
///         println!("Working with right half-open interval: {:?}", half_open);
///     }
/// }
/// ```
///
/// ### Extracting Bounds
/// ```rust
/// use grid1d::intervals::*;
///
/// fn extract_bounds_any_variant(
///     interval: IntervalFinitePositiveLength<f64>
/// ) -> (f64, f64) {
///     // Works uniformly across all variants
///     interval.into_bounds_pair()
/// }
///
/// let closed = IntervalFinitePositiveLength::Closed(IntervalClosed::new(1.0, 2.0));
/// let open = IntervalFinitePositiveLength::Open(IntervalOpen::new(3.0, 4.0));
/// let half_open = IntervalFinitePositiveLength::LowerOpenUpperClosed(
///     IntervalLowerOpenUpperClosed::new(5.0, 6.0)
/// );
///
/// assert_eq!(extract_bounds_any_variant(closed), (1.0, 2.0));
/// assert_eq!(extract_bounds_any_variant(open), (3.0, 4.0));
/// assert_eq!(extract_bounds_any_variant(half_open), (5.0, 6.0));
/// ```
///
/// ## Performance Characteristics
///
/// ### Memory Layout
/// - **Storage**: Two scalar values plus boundary type discriminant
/// - **Size**: Typically 16-24 bytes for f64, depending on alignment
/// - **Alignment**: Natural alignment of the scalar type
/// - **Enum overhead**: Small discriminant tag (1-8 bytes)
///
/// ### Operation Complexity
/// | Operation | Time Complexity | Space Complexity | Notes |
/// |-----------|----------------|------------------|-------|
/// | **Point queries** | O(1) | O(1) | 2 comparisons maximum |
/// | **Containment tests** | O(1) | O(1) | Boundary logic only |
/// | **Intersections** | O(1) | O(1) | Pure computation, may allocate result |
/// | **Geometric calculations** | O(1) | O(1) | Length, midpoint, etc. |
/// | **Pattern matching** | O(1) | O(1) | Compile-time optimized |
/// | **Bound extraction** | O(1) | O(1) | Zero-cost access |
///
/// ### Scalar Type Performance Impact
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::{RealNative64StrictFiniteInDebug, RealNative64StrictFinite, RealScalar};
///
/// // These have identical performance in release builds:
/// type FastPositiveInterval = IntervalFinitePositiveLength<f64>;
/// type OptimalPositiveInterval = IntervalFinitePositiveLength<RealNative64StrictFiniteInDebug>;
///
/// // This has small validation overhead:
/// type SafePositiveInterval = IntervalFinitePositiveLength<RealNative64StrictFinite>;
///
/// // Performance test function (generic over scalar type)
/// fn benchmark_positive_operations<T: RealScalar + Clone>(
///     interval1: &IntervalFinitePositiveLength<T>,
///     interval2: &IntervalFinitePositiveLength<T>
/// ) {
///     // All operations compile to identical assembly for f64 and RealNative64StrictFiniteInDebug
///     let _contains = interval1.contains_interval(interval2);
///     let _intersection = interval1.intersection(interval2);
///     let _length = interval1.length().into_inner();        // O(1) subtraction
///     let _midpoint = interval1.midpoint();   // O(1) arithmetic
///     let _symmetric = interval1.is_symmetric(); // O(1) comparison
///     
///     // Pattern matching is zero-cost
///     match interval1 {
///         IntervalFinitePositiveLength::Closed(_) => { /* ... */ }
///         IntervalFinitePositiveLength::Open(_) => { /* ... */ }
///         IntervalFinitePositiveLength::LowerOpenUpperClosed(_) => { /* ... */ }
///         IntervalFinitePositiveLength::LowerClosedUpperOpen(_) => { /* ... */ }
///     }
/// }
/// ```
///
/// ## Integration with Library Ecosystem
///
/// ### With Numerical Integration
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Adaptive quadrature integration
/// fn adaptive_integrate(
///     f: impl Fn(f64) -> f64 + Copy,
///     domain: &IntervalFinitePositiveLength<f64>,
///     tolerance: f64
/// ) -> f64 {
///     // The integration algorithm is identical regardless of boundary type
///     // because the measure is the same and boundary points have measure zero
///     
///     let length = domain.length().into_inner();
///     let midpoint = domain.midpoint();
///     
///     // Simpson's rule approximation (simplified)
///     let lower = *domain.lower_bound_value();
///     let upper = *domain.upper_bound_value();
///     
///     (length / 6.0) * (f(lower) + 4.0 * f(midpoint) + f(upper))
/// }
///
/// // All boundary types give the same integral (up to numerical precision)
/// let closed_domain = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(0.0, 1.0)
/// );
/// let open_domain = IntervalFinitePositiveLength::Open(
///     IntervalOpen::new(0.0, 1.0)
/// );
///
/// let f = |x: f64| x * x;
/// let integral_closed = adaptive_integrate(f, &closed_domain, 1e-6);
/// let integral_open = adaptive_integrate(f, &open_domain, 1e-6);
///
/// // Should be approximately equal (boundary points don't affect integral)
/// assert!((integral_closed - integral_open).abs() < 1e-10);
/// ```
///
/// ### With Mesh Generation
/// ```rust
/// use grid1d::intervals::*;
///
/// /// Generate uniform mesh points within any positive-length interval
/// fn generate_uniform_mesh(
///     domain: &IntervalFinitePositiveLength<f64>,
///     num_points: usize
/// ) -> Vec<f64> {
///     if num_points <= 1 {
///         return vec![domain.midpoint()];
///     }
///     
///     let lower = *domain.lower_bound_value();
///     let upper = *domain.upper_bound_value();
///     let step = (upper - lower) / (num_points - 1) as f64;
///     
///     (0..num_points)
///         .map(|i| lower + i as f64 * step)
///         .collect()
/// }
///
/// /// Generate mesh respecting boundary conditions
/// fn generate_boundary_aware_mesh(
///     domain: &IntervalFinitePositiveLength<f64>,
///     num_interior_points: usize
/// ) -> Vec<f64> {
///     let lower = *domain.lower_bound_value();
///     let upper = *domain.upper_bound_value();
///     
///     let mut points = Vec::new();
///     
///     // Include boundary points based on interval type
///     match domain {
///         IntervalFinitePositiveLength::Closed(_) => {
///             points.push(lower);
///         }
///         IntervalFinitePositiveLength::LowerClosedUpperOpen(_) => {
///             points.push(lower);
///         }
///         _ => {} // Open boundaries not included
///     }
///     
///     // Add interior points
///     for i in 1..=num_interior_points {
///         let t = i as f64 / (num_interior_points + 1) as f64;
///         points.push(lower + t * (upper - lower));
///     }
///     
///     // Include upper boundary if appropriate
///     match domain {
///         IntervalFinitePositiveLength::Closed(_) => {
///             points.push(upper);
///         }
///         IntervalFinitePositiveLength::LowerOpenUpperClosed(_) => {
///             points.push(upper);
///         }
///         _ => {} // Open boundaries not included
///     }
///     
///     points
/// }
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`IntervalFinitePositiveLength`]
/// - **Numerical analysis**: When you need positive-length domains but want to abstract boundary behavior
/// - **Generic algorithms**: Writing code that works with any finite positive interval
/// - **API design**: Accepting any bounded interval while excluding singletons
/// - **Mathematical modeling**: When the boundary type may vary but the domain must be non-degenerate
///
/// ### When to Use Specific Variants
/// - **Closed**: Integration bounds, probability ranges, physical constraints with inclusive limits
/// - **Open**: Limit domains, convergence regions, strict inequality constraints  
/// - **Half-open**: Partitions, cumulative functions, discrete-continuous interfaces
///
/// ### Recommended Patterns
/// ```rust
/// use grid1d::intervals::*;
/// use into_inner::IntoInner;
/// use num_valid::{RealNative64StrictFiniteInDebug, functions::Rounding};
/// use num::Zero;
///
/// // Define type aliases for your domain
/// type Real = RealNative64StrictFiniteInDebug;
/// type PositiveInterval = IntervalFinitePositiveLength<Real>;
/// type ClosedInterval = IntervalClosed<Real>;
///
/// // Pattern: Handle all variants uniformly when possible
/// fn compute_interval_statistics(interval: &PositiveInterval) -> (Real, Real, Real) {
///     let length = interval.length().into_inner();
///     let midpoint = interval.midpoint();
///     let lower = interval.lower_bound_value().clone();
///     
///     (length, midpoint, lower) // Same calculation for all variants
/// }
///
/// // Pattern: Specialize behavior only when boundary semantics matter
/// fn count_integer_endpoints(interval: &PositiveInterval) -> usize {
///     let lower = interval.lower_bound_value().clone();
///     let upper = interval.upper_bound_value().clone();
///     let lower_int = lower.clone().kernel_floor().into_inner() as i64;
///     let upper_int = upper.clone().kernel_ceil().into_inner() as i64;
///     
///     let mut count = 0;
///     
///     // Check lower boundary
///     if lower == lower_int as f64 && interval.contains_point(&lower) {
///         count += 1;
///     }
///     
///     // Check upper boundary  
///     if upper == upper_int as f64 && interval.contains_point(&upper) {
///         count += 1;
///     }
///     
///     count
/// }
///
/// // Pattern: Convert to most specific type when needed
/// fn require_closed_interval(interval: PositiveInterval) -> Result<ClosedInterval, String> {
///     match interval {
///         IntervalFinitePositiveLength::Closed(closed) => Ok(closed),
///         _ => Err("Algorithm requires closed interval for boundary conditions".to_string()),
///     }
/// }
/// ```
///
/// ## Error Handling and Edge Cases
///
/// ### Construction Validation
/// ```rust,ignore
/// use grid1d::intervals::*;
///
/// // All variants validate that lower < upper during construction
/// assert!(IntervalClosed::try_new(1.0, 0.0).is_err());  // Invalid bounds
/// assert!(IntervalOpen::try_new(1.0, 1.0).is_err());    // Zero length
/// assert!(IntervalLowerOpenUpperClosed::try_new(2.0, 1.0).is_err()); // Inverted
///
/// // In debug mode, invalid construction panics
/// #[cfg(debug_assertions)]
/// {
///     std::panic::catch_unwind(|| {
///         IntervalClosed::new(1.0, 0.0); // Panics in debug
///     }).expect_err("Should panic on invalid bounds");
/// }
/// ```
///
/// ### Numerical Precision Considerations
/// ```rust
/// use grid1d::intervals::*;
///
/// // CAUTION: Very small intervals near floating-point precision limits
/// let tiny_interval = IntervalFinitePositiveLength::Closed(
///     IntervalClosed::new(1.0, 1.0 + f64::EPSILON)
/// );
///
/// // Length is positive but very small
/// assert!(tiny_interval.length().into_inner() > 0.0);
/// assert!(tiny_interval.length().into_inner() < 1e-15);
///
/// // Midpoint calculation may have precision issues
/// let midpoint = tiny_interval.midpoint();
/// let expected_mid = 1.0 + f64::EPSILON / 2.0;
///
/// // Use appropriate tolerance for comparisons
/// assert!((midpoint - expected_mid).abs() < f64::EPSILON);
/// ```
///
/// ### Common Pitfalls and Solutions
/// ```rust
/// use grid1d::intervals::*;
///
/// // WRONG: Assuming all intervals include their endpoints
/// fn naive_endpoint_processing(interval: &IntervalFinitePositiveLength<f64>) -> (f64, f64) {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     // This ignores boundary semantics!
///     (lower, upper) // May include points not actually in the interval
/// }
///
/// // CORRECT: Check boundary inclusion when it matters
/// fn proper_endpoint_processing(interval: &IntervalFinitePositiveLength<f64>) -> Vec<f64> {
///     let lower = *interval.lower_bound_value();
///     let upper = *interval.upper_bound_value();
///     
///     let mut endpoints = Vec::new();
///     
///     if interval.contains_point(&lower) {
///         endpoints.push(lower);
///     }
///     if interval.contains_point(&upper) {
///         endpoints.push(upper);
///     }
///     
///     endpoints
/// }
///
/// let closed = IntervalFinitePositiveLength::Closed(IntervalClosed::new(0.0, 1.0));
/// let open = IntervalFinitePositiveLength::Open(IntervalOpen::new(0.0, 1.0));
///
/// assert_eq!(proper_endpoint_processing(&closed), vec![0.0, 1.0]); // Both endpoints
/// assert!(proper_endpoint_processing(&open).is_empty()); // No endpoints
/// ```
///
/// ## Mathematical Correctness Guarantees
///
/// The [`IntervalFinitePositiveLength`] enum maintains all essential mathematical properties:
///
/// - **Positive measure**: `μ(I) = |upper - lower| > 0` for all variants
/// - **Containment transitivity**: If A ⊆ B and B ⊆ C, then A ⊆ C
/// - **Intersection commutativity**: A ∩ B = B ∩ A
/// - **Intersection associativity**: (A ∩ B) ∩ C = A ∩ (B ∩ C)
/// - **Monotonicity**: If A ⊆ B, then μ(A) ≤ μ(B)
/// - **Boundary consistency**: Containment and intersection operations respect boundary semantics
/// - **Geometric invariance**: Length, midpoint, and symmetry properties are boundary-independent
///
/// These properties are preserved regardless of the scalar type used and are guaranteed
/// by the underlying implementations of the [`IntervalTrait`] and [`IntervalFinitePositiveLengthTrait`].
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
#[serde(bound(deserialize = "RealType: for<'a> Deserialize<'a>"))]
pub enum IntervalFinitePositiveLength<RealType: RealScalar> {
    /// A ***closed*** interval `[a,b]` where both endpoints are included.
    ///
    /// This variant represents the most restrictive form of inclusion, where both
    /// boundary points `a` and `b` are considered part of the interval. The interval
    /// contains all real numbers `x` such that `a ≤ x ≤ b`.
    ///
    /// **Mathematical properties:**
    /// - **Compactness**: Always compact (closed and bounded)
    /// - **Boundary inclusion**: Both `a` and `b` are contained
    /// - **Supremum/Infimum**: Both achieved at the endpoints
    /// - **Typical use cases**: Integration bounds, probability ranges, physical constraints
    Closed(IntervalClosed<RealType>),

    /// An ***open*** interval `(a,b)` where both endpoints are excluded.
    ///
    /// This variant represents the most permissive form of exclusion, where neither
    /// boundary point `a` nor `b` is considered part of the interval. The interval
    /// contains all real numbers `x` such that `a < x < b`.
    ///
    /// **Mathematical properties:**
    /// - **Openness**: Open set in the real line topology
    /// - **Boundary exclusion**: Neither `a` nor `b` are contained
    /// - **Supremum/Infimum**: Exist but are not achieved
    /// - **Typical use cases**: Limit domains, convergence regions, strict inequalities
    Open(IntervalOpen<RealType>),

    /// A ***left half-open*** interval `(a,b]` where the lower endpoint is excluded and the upper endpoint is included.
    ///
    /// This variant represents a mixed boundary condition where the lower bound `a`
    /// is not part of the interval but the upper bound `b` is included. The interval
    /// contains all real numbers `x` such that `a < x ≤ b`.
    ///
    /// **Mathematical properties:**
    /// - **Boundary behavior**: Excludes `a`, includes `b`
    /// - **Topology**: Neither open nor closed (half-open)
    /// - **Right-continuity**: Natural for right-continuous functions
    /// - **Typical use cases**: Cumulative distributions, partition elements
    LowerOpenUpperClosed(IntervalLowerOpenUpperClosed<RealType>),

    /// A ***right half-open*** interval `[a,b)` where the lower endpoint is included and the upper endpoint is excluded.
    ///
    /// This variant represents a mixed boundary condition where the lower bound `a`
    /// is part of the interval but the upper bound `b` is not included. The interval
    /// contains all real numbers `x` such that `a ≤ x < b`.
    ///
    /// **Mathematical properties:**
    /// - **Boundary behavior**: Includes `a`, excludes `b`
    /// - **Topology**: Neither open nor closed (half-open)
    /// - **Left-continuity**: Natural for left-continuous functions
    /// - **Typical use cases**: Array indexing ranges, interval partitions, discrete-continuous boundaries
    LowerClosedUpperOpen(IntervalLowerClosedUpperOpen<RealType>),
}

impl<RealType: RealScalar> IntervalFinitePositiveLength<RealType> {
    /// Constructs a finite-length interval with positive length from runtime-typed bounds.
    ///
    /// Similar to [`IntervalFiniteLength::try_from_runtime_bounds`], but excludes singleton intervals by
    /// construction. This function is used internally when the caller knows the interval
    /// must have positive length (e.g., in hull computation between bounded intervals).
    ///
    /// # Arguments
    ///
    /// * `lower_bound` - Runtime-typed lower bound (closed or open)
    /// * `upper_bound` - Runtime-typed upper bound (closed or open)
    ///
    /// # Bound Type Selection
    ///
    /// | Lower | Upper | Result Type | Notation |
    /// |-------|-------|-------------|----------|
    /// | Closed | Closed | [`IntervalClosed`] | `[a, b]` where `a < b` |
    /// | Closed | Open | [`IntervalLowerClosedUpperOpen`] | `[a, b)` |
    /// | Open | Closed | [`IntervalLowerOpenUpperClosed`] | `(a, b]` |
    /// | Open | Open | [`IntervalOpen`] | `(a, b)` |
    ///
    /// # Returns
    ///
    /// * `Ok(IntervalFinitePositiveLength)` - Successfully constructed interval with `lower < upper`
    /// * `Err(ErrorsIntervalConstruction)` - If `lower_bound ≥ upper_bound`
    ///
    /// # Difference from `IntervalFiniteLength::try_from_runtime_bounds`
    ///
    /// Unlike [`IntervalFiniteLength::try_from_runtime_bounds`], this function:
    /// - **Does not** check for singleton case (assumes `lower ≠ upper`)
    /// - **Always** creates positive-length intervals
    /// - Used when singleton case is excluded by construction
    ///
    /// # Examples
    ///
    /// ```ignore
    /// // Note: This is a private function used internally
    /// use grid1d::intervals::*;
    /// use grid1d::bounds::*;
    /// use try_create::New;
    ///
    /// // Half-open interval [0, 1)
    /// let lower: LowerBoundRuntime<f64> = LowerBoundClosed::new(0.0).into();
    /// let upper: UpperBoundRuntime<f64> = UpperBoundOpen::new(1.0).into();
    /// let interval = IntervalFinitePositiveLength::try_from_runtime_bounds(lower, upper).unwrap();
    /// ```
    ///
    /// # Use Cases
    ///
    /// - Hull computation between bounded intervals (always produces positive-length result)
    /// - Operations known to produce non-singleton intervals
    /// - Performance optimization when singleton check is unnecessary
    pub(crate) fn try_from_runtime_bounds(
        lower_bound: LowerBoundRuntime<RealType>,
        upper_bound: UpperBoundRuntime<RealType>,
    ) -> Result<Self, ErrorsIntervalConstruction<RealType>> {
        // Delegate to the appropriate concrete interval type's try_new() method
        // This ensures consistent validation and maintains a single source of truth
        match (lower_bound, upper_bound) {
            (LowerBoundRuntime::Closed(lower_bound), UpperBoundRuntime::Closed(upper_bound)) => {
                Ok(IntervalClosed::try_from_static_bounds(lower_bound, upper_bound)?.into())
            }
            (LowerBoundRuntime::Closed(lower_bound), UpperBoundRuntime::Open(upper_bound)) => Ok(
                IntervalLowerClosedUpperOpen::try_from_static_bounds(lower_bound, upper_bound)?
                    .into(),
            ),
            (LowerBoundRuntime::Open(lower_bound), UpperBoundRuntime::Closed(upper_bound)) => Ok(
                IntervalLowerOpenUpperClosed::try_from_static_bounds(lower_bound, upper_bound)?
                    .into(),
            ),
            (LowerBoundRuntime::Open(lower_bound), UpperBoundRuntime::Open(upper_bound)) => {
                Ok(IntervalOpen::try_from_static_bounds(lower_bound, upper_bound)?.into())
            }
        }
    }
}

impl<RealType: RealScalar> IntervalOperations for IntervalFinitePositiveLength<RealType> {}

impl<RealType: RealScalar> Contains for IntervalFinitePositiveLength<RealType> {
    /// Returns `true` if the 1D point/coordinate `x` is contained in the current [`IntervalFinitePositiveLength`].
    #[inline(always)]
    fn contains_point(&self, x: &RealType) -> bool {
        match self {
            IntervalFinitePositiveLength::Closed(interval) => interval.contains_point(x),
            IntervalFinitePositiveLength::Open(interval) => interval.contains_point(x),
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.contains_point(x)
            }
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.contains_point(x)
            }
        }
    }
}

impl<RealType: RealScalar> IntervalTrait for IntervalFinitePositiveLength<RealType> {}
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
/// A singleton interval `[a]` containing exactly one point.
///
/// The [`IntervalSingleton`] represents the most degenerate form of interval in mathematical
/// analysis - a set containing exactly one point. While mathematically an interval typically
/// implies a range of values, singleton intervals are essential for representing point masses,
/// discrete events, and boundary conditions in mathematical modeling.
///
/// ## Mathematical Properties
///
/// - **Notation**: `[a]` or `{a}` - a single point
/// - **Definition**: `{x ∈ ℝ : x = a}` - the set containing only the value `a`
/// - **Measure**: Always `0` (zero length in Lebesgue measure)
/// - **Cardinality**: Exactly `1` (finite set)
/// - **Topology**: Isolated point with empty interior
/// - **Connectedness**: Trivially connected (single point)
/// - **Compactness**: Always compact (finite sets are compact)
///
/// ## Use Cases
///
/// Singleton intervals are essential when:
/// - **Point masses**: Representing discrete probability masses at specific values
/// - **Boundary conditions**: Modeling fixed points in differential equations
/// - **Discrete events**: Representing events that occur at exact moments or values
/// - **Degenerate cases**: Handling limit cases where intervals collapse to points
/// - **Constraint satisfaction**: Representing exact equality constraints
/// - **Intersection results**: When two intervals meet at exactly one point
///
/// ## Comparison with Other Interval Types
///
/// | Property | Singleton `[a]` | Closed `[a,b]` | Open `(a,b)` |
/// |----------|-----------------|----------------|--------------|
/// | Contains `a` | ✓ | ✓ | ✗ |
/// | Contains other points | ✗ | ✓ (if `b > a`) | ✓ (if `b > a`) |
/// | Length | `0` | `b - a` | `b - a` |
/// | Measure | `0` | `b - a` | `b - a` |
/// | Cardinality | `1` | Uncountably infinite | Uncountably infinite |
/// | Compactness | Always | Always | Never |
/// | Interior | Empty set `∅` | `(a, b)` | `(a, b)` |
///
/// ## Construction
///
/// ### Safe Construction (Recommended)
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::{New, TryNew};
///
/// // Use the recommended performance-optimal type
/// type Real = RealNative64StrictFiniteInDebug;
///
/// // Create singleton with validated scalar
/// let singleton = IntervalSingleton::new(Real::try_new(42.0).unwrap());
/// assert_eq!(singleton.value(), &Real::try_new(42.0).unwrap());
/// ```
///
/// ### Raw Construction (Use with Caution)
/// ```rust
/// use grid1d::intervals::IntervalSingleton;
/// use try_create::New;
///
/// // Direct construction with raw f64 - only when you trust your data
/// let singleton = IntervalSingleton::new(3.14159);
/// assert_eq!(singleton.value(), &3.14159);
/// ```
///
/// ## Key Operations
///
/// ### Point Containment (Exact Match Only)
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let singleton = IntervalSingleton::new(42.0);
///
/// // Only the exact value is contained
/// assert!(singleton.contains_point(&42.0));      // Exact match
/// assert!(!singleton.contains_point(&42.1));     // Any other value
/// assert!(!singleton.contains_point(&41.9));     // Any other value
/// assert!(!singleton.contains_point(&0.0));      // Any other value
///
/// // Floating-point precision considerations
/// let epsilon = f64::EPSILON;
/// assert!(!singleton.contains_point(&(42.0 * (1.0 + epsilon)))); // Not exactly 42.0
/// ```
///
/// ### Value Access and Consumption
/// ```rust
/// use grid1d::intervals::*;
/// use into_inner::IntoInner;
/// use try_create::New;
///
/// let singleton = IntervalSingleton::new(123.456);
///
/// // Access the value by reference
/// assert_eq!(singleton.value(), &123.456);
///
/// // Consume the singleton to get the value
/// let value = singleton.into_inner();
/// assert_eq!(value, 123.456);
/// ```
///
/// ### Interval Containment Testing
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let singleton_a = IntervalSingleton::new(5.0);
/// let singleton_b = IntervalSingleton::new(5.0);
/// let singleton_c = IntervalSingleton::new(10.0);
/// let interval = IntervalClosed::new(0.0, 10.0);
///
/// // Singleton contains itself if values are equal
/// assert!(singleton_a.contains_interval(&singleton_b));  // Both contain 5.0
/// assert!(singleton_b.contains_interval(&singleton_a));  // Symmetric
///
/// // Different singletons don't contain each other
/// assert!(!singleton_a.contains_interval(&singleton_c)); // 5.0 ≠ 10.0
/// assert!(!singleton_c.contains_interval(&singleton_a)); // 10.0 ≠ 5.0
///
/// // Larger intervals can contain singletons
/// assert!(interval.contains_interval(&singleton_a));     // [0,10] contains {5}
/// assert!(interval.contains_interval(&singleton_c));     // [0,10] contains {10}
///
/// // Singletons cannot contain larger intervals
/// assert!(!singleton_a.contains_interval(&interval));    // {5} cannot contain [0,10]
/// ```
///
/// ## Intersection Examples
///
/// ### Singleton with Singleton
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let singleton_a = IntervalSingleton::new(7.0);
/// let singleton_b = IntervalSingleton::new(7.0);
/// let singleton_c = IntervalSingleton::new(8.0);
///
/// // Intersection of identical singletons
/// if let Some(intersection) = singleton_a.intersection(&singleton_b) {
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(result)) => {
///             assert_eq!(result.value(), &7.0);
///         }
///         _ => unreachable!("Expected singleton intersection"),
///     }
/// }
///
/// // Intersection of different singletons
/// assert!(singleton_a.intersection(&singleton_c).is_none()); // {7} ∩ {8} = ∅
/// ```
///
/// ### Singleton with Interval
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let singleton = IntervalSingleton::new(2.5);
/// let closed = IntervalClosed::new(2.0, 3.0);
/// let open = IntervalOpen::new(2.0, 3.0);
/// let outside = IntervalClosed::new(4.0, 5.0);
///
/// // Singleton intersects interval if contained
/// if let Some(intersection) = singleton.intersection(&closed) {
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(result)) => {
///             assert_eq!(result.value(), &2.5); // Returns the singleton
///         }
///         _ => unreachable!("Expected singleton intersection"),
///     }
/// }
///
/// // Same result for open interval (since 2.5 is in interior)
/// if let Some(intersection) = singleton.intersection(&open) {
///     match intersection {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(result)) => {
///             assert_eq!(result.value(), &2.5);
///         }
///         _ => unreachable!("Expected singleton intersection"),
///     }
/// }
///
/// // No intersection when singleton is outside interval
/// assert!(singleton.intersection(&outside).is_none()); // {2.5} ∩ [4,5] = ∅
/// ```
///
/// ### Edge Cases: Boundary Intersections
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let boundary_singleton = IntervalSingleton::new(1.0);
/// let closed = IntervalClosed::new(1.0, 2.0);        // [1, 2]
/// let open = IntervalOpen::new(1.0, 2.0);            // (1, 2)
/// let half_open_left = IntervalLowerOpenUpperClosed::new(1.0, 2.0); // (1, 2]
/// let half_open_right = IntervalLowerClosedUpperOpen::new(1.0, 2.0); // [1, 2)
///
/// // Intersection depends on boundary inclusion
/// assert!(boundary_singleton.intersection(&closed).is_some());      // {1} ∩ [1,2] = {1}
/// assert!(boundary_singleton.intersection(&open).is_none());        // {1} ∩ (1,2) = ∅
/// assert!(boundary_singleton.intersection(&half_open_left).is_none()); // {1} ∩ (1,2] = ∅
/// assert!(boundary_singleton.intersection(&half_open_right).is_some()); // {1} ∩ [1,2) = {1}
/// ```
///
/// ## Advanced Usage Patterns
///
/// ### Generic Programming with Different Scalar Types
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::{RealScalar, RealNative64StrictFinite};
/// use try_create::{New, TryNew};
///
/// // Generic function works with any RealScalar type
/// fn analyze_singleton<T: RealScalar + std::fmt::Display>(
///     singleton: &IntervalSingleton<T>
/// ) -> String {
///     format!("Singleton interval containing the single point: {}", singleton.value())
/// }
///
/// // Works with validated types
/// let safe_singleton = IntervalSingleton::new(
///     RealNative64StrictFinite::try_new(2.71828).unwrap()
/// );
/// println!("{}", analyze_singleton(&safe_singleton));
///
/// // Also works with raw f64
/// let raw_singleton = IntervalSingleton::new(2.71828);
/// println!("{}", analyze_singleton(&raw_singleton));
/// ```
///
/// ### Conversion to Other Interval Types
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let singleton = IntervalSingleton::new(1.5);
///
/// // Convert to category enums
/// let finite_zero: IntervalFiniteLength<f64> = singleton.clone().into();
/// let general: Interval<f64> = singleton.into();
///
/// // Pattern match on the general enum
/// match general {
///     Interval::FiniteLength(IntervalFiniteLength::ZeroLength(recovered_singleton)) => {
///         assert_eq!(recovered_singleton.value(), &1.5);
///     }
///     _ => unreachable!("Should be singleton interval"),
/// }
/// ```
///
/// ## Mathematical Applications
///
/// ### Probability Theory: Point Masses
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// /// Represents a discrete probability distribution with point masses
/// struct DiscreteDistribution {
///     point_masses: Vec<(IntervalSingleton<f64>, f64)>, // (point, probability)
/// }
///
/// impl DiscreteDistribution {
///     fn new() -> Self {
///         Self { point_masses: Vec::new() }
///     }
///     
///     /// Add a point mass at a specific value
///     fn add_point_mass(&mut self, value: f64, probability: f64) {
///         self.point_masses.push((IntervalSingleton::new(value), probability));
///     }
///     
///     /// Probability mass function
///     fn pmf(&self, x: f64) -> f64 {
///         for (point, prob) in &self.point_masses {
///             if point.contains_point(&x) {
///                 return *prob;
///             }
///         }
///         0.0 // No mass at this point
///     }
///     
///     /// Support: all points with non-zero probability
///     fn support(&self) -> Vec<f64> {
///         self.point_masses.iter()
///             .map(|(point, _)| *point.value())
///             .collect()
///     }
/// }
///
/// let mut discrete = DiscreteDistribution::new();
/// discrete.add_point_mass(1.0, 0.3);
/// discrete.add_point_mass(2.0, 0.5);
/// discrete.add_point_mass(3.0, 0.2);
///
/// assert_eq!(discrete.pmf(1.0), 0.3);
/// assert_eq!(discrete.pmf(1.5), 0.0);  // No mass between points
/// assert_eq!(discrete.pmf(2.0), 0.5);
/// assert_eq!(discrete.support(), vec![1.0, 2.0, 3.0]);
/// ```
///
/// ### Constraint Satisfaction: Exact Equality
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// /// Represents an exact equality constraint in optimization
/// struct EqualityConstraint {
///     feasible_set: IntervalSingleton<f64>
/// }
///
/// impl EqualityConstraint {
///     /// Create constraint: x = target_value
///     fn new(target_value: f64) -> Self {
///         Self {
///             feasible_set: IntervalSingleton::new(target_value)
///         }
///     }
///     
///     /// Check if a point satisfies the exact equality constraint
///     fn is_satisfied(&self, x: f64) -> bool {
///         self.feasible_set.contains_point(&x)
///     }
///     
///     /// Get the target value
///     fn target_value(&self) -> f64 {
///         *self.feasible_set.value()
///     }
///     
///     /// Constraint violation (distance from feasible point)
///     fn violation(&self, x: f64) -> f64 {
///         (x - self.target_value()).abs()
///     }
///     
///     /// Project a point onto the feasible set (returns target value)
///     fn project(&self, _x: f64) -> f64 {
///         self.target_value() // Always project to the single feasible point
///     }
/// }
///
/// let constraint = EqualityConstraint::new(5.0); // x = 5
///
/// assert!(constraint.is_satisfied(5.0));     // Exact match
/// assert!(!constraint.is_satisfied(5.001));  // Any deviation
/// assert!(!constraint.is_satisfied(4.999));  // Any deviation
///
/// assert_eq!(constraint.target_value(), 5.0);
/// assert_eq!(constraint.violation(5.0), 0.0);
/// assert_eq!(constraint.violation(6.0), 1.0);
/// assert_eq!(constraint.project(100.0), 5.0); // Always projects to target
/// ```
///
/// ### Dynamical Systems: Fixed Points
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// /// Represents a fixed point in a dynamical system
/// struct FixedPoint {
///     point: IntervalSingleton<f64>
/// }
///
/// impl FixedPoint {
///     fn new(value: f64) -> Self {
///         Self {
///             point: IntervalSingleton::new(value)
///         }
///     }
///     
///     /// Check if a trajectory passes through the fixed point
///     fn is_visited_by_trajectory(&self, trajectory: &[f64]) -> bool {
///         trajectory.iter().any(|&x| self.point.contains_point(&x))
///     }
///     
///     /// Distance from a point to the fixed point
///     fn distance_from(&self, x: f64) -> f64 {
///         (x - self.point.value()).abs()
///     }
///     
///     /// Check if a point is within epsilon of the fixed point
///     fn is_near(&self, x: f64, epsilon: f64) -> bool {
///         self.distance_from(x) < epsilon
///     }
/// }
///
/// // Example: fixed point of f(x) = x/2 + 1 is x = 2
/// let fixed_point = FixedPoint::new(2.0);
///
/// let trajectory = vec![4.0, 3.0, 2.5, 2.25, 2.125]; // Converging to 2.0
/// assert!(!fixed_point.is_visited_by_trajectory(&trajectory)); // Doesn't hit exactly
///
/// let exact_trajectory = vec![4.0, 3.0, 2.0, 2.0, 2.0]; // Hits fixed point
/// assert!(fixed_point.is_visited_by_trajectory(&exact_trajectory));
///
/// assert!((fixed_point.distance_from(2.1) - 0.1).abs() < 1.0e-15);
/// assert!(fixed_point.is_near(1.999, 0.01));
/// assert!(!fixed_point.is_near(1.99, 0.001));
/// ```
///
/// ### Signal Processing: Impulse Functions
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// /// Represents a Dirac delta function approximation (impulse at specific point)
/// struct DiracDelta {
///     impulse_location: IntervalSingleton<f64>,
///     amplitude: f64,
/// }
///
/// impl DiracDelta {
///     fn new(location: f64, amplitude: f64) -> Self {
///         Self {
///             impulse_location: IntervalSingleton::new(location),
///             amplitude,
///         }
///     }
///     
///     /// Sample the delta function (infinite at impulse location, zero elsewhere)
///     fn sample(&self, x: f64) -> f64 {
///         if self.impulse_location.contains_point(&x) {
///             self.amplitude // In practice, this represents infinite amplitude
///         } else {
///             0.0
///         }
///     }
///     
///     /// Get the location of the impulse
///     fn location(&self) -> f64 {
///         *self.impulse_location.value()
///     }
///     
///     /// Integrate over any interval containing the impulse
///     fn integrate_over(&self, interval: &impl IntervalTrait<RealType = f64>) -> f64 {
///         if interval.contains_interval(&self.impulse_location) {
///             self.amplitude // Integral of delta function is its amplitude
///         } else {
///             0.0 // No impulse in this interval
///         }
///     }
/// }
///
/// let impulse = DiracDelta::new(0.0, 1.0); // Unit impulse at t=0
///
/// assert_eq!(impulse.sample(0.0), 1.0);     // At impulse location
/// assert_eq!(impulse.sample(0.001), 0.0);   // Anywhere else
/// assert_eq!(impulse.location(), 0.0);
///
/// let containing_interval = IntervalClosed::new(-1.0, 1.0);
/// let outside_interval = IntervalClosed::new(1.0, 2.0);
///
/// assert_eq!(impulse.integrate_over(&containing_interval), 1.0); // Contains impulse
/// assert_eq!(impulse.integrate_over(&outside_interval), 0.0);    // Doesn't contain impulse
/// ```
///
/// ## Performance Characteristics
///
/// ### Memory Layout
/// - **Storage**: Single scalar value (8 bytes for f64)
/// - **Alignment**: Natural alignment of the scalar type
/// - **No heap allocation**: All data stored inline
/// - **Minimal overhead**: Most compact interval representation
///
/// ### Operation Complexity
/// | Operation | Time Complexity | Notes |
/// |-----------|-----------------|-------|
/// | **Point queries** | O(1) | Single equality comparison |
/// | **Containment tests** | O(1) | Equality check only |
/// | **Intersections** | O(1) | Equality or emptiness check |
/// | **Value access** | O(1) | Direct field access |
/// | **Construction** | O(1) | Single validation check in debug mode |
///
/// ### Comparison with Other Interval Types
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// // Singleton has the smallest memory footprint
/// let singleton = IntervalSingleton::new(42.0);
/// let open = IntervalOpen::new(41.9, 42.1); // Approximation
///
/// // Closed interval with equal lower and upper bounds (degenerate interval) cannot be built
///
/// // All represent roughly the same mathematical concept but:
/// // - Singleton: 8 bytes, exact point representation
/// // - Open: 16+ bytes, approximate neighborhood
///
/// let n = 1000000;
/// let one_div_n = 1.0e-6;
///
/// // Operations on singleton are fastest
/// let start = std::time::Instant::now();
/// for i in 0..n {
///     let x = i as f64 * one_div_n;
///     let _ = singleton.contains_point(&x); // Single equality check
/// }
/// let singleton_time = start.elapsed();
///
/// let start = std::time::Instant::now();
/// for i in 0..n {
///     let x = i as f64 * one_div_n;
///     let _ = open.contains_point(&x); // Two comparisons + equality check
/// }
/// let open_time = start.elapsed();
///
/// // singleton_time should be faster than open_time
/// ```
///
/// ## Error Handling and Edge Cases
///
/// ### Construction Validation
/// ```rust,ignore
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// // Valid construction
/// let singleton = IntervalSingleton::new(123.456);
/// assert_eq!(singleton.value(), &123.456);
///
/// #[cfg(debug_assertions)]
/// {
///     // Debug mode validation catches invalid values
///     std::panic::catch_unwind(|| {
///         IntervalSingleton::new(f64::NAN); // Panics in debug
///     }).expect_err("Should panic on NaN");
///     
///     std::panic::catch_unwind(|| {
///         IntervalSingleton::new(f64::INFINITY); // Panics in debug
///     }).expect_err("Should panic on infinity");
/// }
/// ```
///
/// ### Floating-Point Precision Considerations
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// let singleton = IntervalSingleton::new(0.1);
///
/// // Exact comparison required for containment
/// let computed_value = 0.1_f64; // This might have precision issues
/// let sum_value = 0.05_f64 + 0.05_f64; // This might be different
///
/// // These might not be equal due to floating-point representation
/// if singleton.contains_point(&computed_value) {
///     println!("Exact match");
/// } else {
///     println!("Precision difference: {} vs {}", singleton.value(), computed_value);
/// }
///
/// // For approximate comparisons, use interval neighborhoods instead
/// let epsilon = 1e-10;
/// let neighborhood = IntervalClosed::new(0.1 - epsilon, 0.1 + epsilon);
/// assert!(neighborhood.contains_point(&computed_value)); // More robust
/// ```
///
/// ### Common Pitfalls and Solutions
/// ```rust
/// use grid1d::intervals::*;
/// use try_create::New;
///
/// // WRONG: Expecting singleton to behave like an interval with length
/// fn naive_midpoint(singleton: &IntervalSingleton<f64>) -> f64 {
///     // This doesn't make sense - singleton has no "middle"
///     *singleton.value() // The value IS the midpoint
/// }
///
/// // CORRECT: Understanding that singleton represents a single point
/// fn get_singleton_value(singleton: &IntervalSingleton<f64>) -> f64 {
///     *singleton.value() // The single point value
/// }
///
/// // WRONG: Trying to sample "within" a singleton
/// fn naive_sample_from_singleton(singleton: &IntervalSingleton<f64>) -> Vec<f64> {
///     // Can't sample multiple distinct points from a singleton!
///     vec![*singleton.value(); 10] // All samples are the same
/// }
///
/// // CORRECT: Understanding that singleton contains only one point
/// fn sample_singleton_correctly(singleton: &IntervalSingleton<f64>, n: usize) -> Vec<f64> {
///     vec![*singleton.value(); n] // All samples must be the same value
/// }
///
/// let singleton = IntervalSingleton::new(7.0);
/// assert_eq!(get_singleton_value(&singleton), 7.0);
/// assert_eq!(sample_singleton_correctly(&singleton, 5), vec![7.0, 7.0, 7.0, 7.0, 7.0]);
/// ```
///
/// ## Best Practices and Recommendations
///
/// ### When to Use [`IntervalSingleton`]
/// - **Point masses**: Discrete probability distributions
/// - **Exact constraints**: Equality constraints in optimization
/// - **Fixed points**: Dynamical systems analysis
/// - **Boundary conditions**: Specific value requirements
/// - **Degenerate limits**: When intervals collapse to points
/// - **Discrete events**: Events at specific values/times
///
/// ### When NOT to Use [`IntervalSingleton`]
/// - **Ranges of values**: Use proper intervals like [`IntervalClosed`]
/// - **Approximate equality**: Use small intervals instead
/// - **Continuous domains**: Use positive-length intervals
/// - **Tolerances**: Use intervals with appropriate width
///
/// ### Recommended Patterns
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealNative64StrictFiniteInDebug;
/// use try_create::{New, TryNew};
/// use into_inner::IntoInner;
///
/// // Define type aliases for your application
/// type Real = RealNative64StrictFiniteInDebug;
/// type Singleton = IntervalSingleton<Real>;
///
/// // Pattern: Use for exact point representations
/// fn create_point_constraint(target: Real) -> Singleton {
///     Singleton::new(target)
/// }
///
/// // Pattern: Convert to interval when you need neighborhoods
/// fn expand_to_neighborhood(singleton: Singleton, epsilon: Real) -> IntervalClosed<Real> {
///     let value = singleton.into_inner();
///     IntervalClosed::new(value.clone() - &epsilon, value + epsilon)
/// }
///
/// // Pattern: Check exact equality, not approximate
/// fn is_exactly_at_point(singleton: &Singleton, candidate: &Real) -> bool {
///     singleton.contains_point(candidate) // Exact equality required
/// }
///
/// // Pattern: Handle singleton vs interval cases explicitly
/// fn handle_interval_or_singleton(input: Interval<Real>) -> String {
///     match input {
///         Interval::FiniteLength(IntervalFiniteLength::ZeroLength(singleton)) => {
///             format!("Exact point: {}", singleton.value())
///         }
///         Interval::FiniteLength(IntervalFiniteLength::PositiveLength(interval)) => {
///             format!("Range: [{}, {}]", interval.lower_bound_value(), interval.upper_bound_value())
///         }
///         Interval::InfiniteLength(_) => {
///             "Unbounded interval".to_string()
///         }
///     }
/// }
/// ```
///
/// ## Mathematical Correctness Guarantees
///
/// The [`IntervalSingleton`] struct maintains all essential mathematical properties:
///
/// - **Set equality**: Two singletons are equal if and only if they contain the same value
/// - **Containment**: A singleton contains another interval if and only if that interval is the same singleton
/// - **Intersection**: The intersection of a singleton with any interval is either the singleton (if contained) or empty
/// - **Measure**: The measure is always exactly zero
/// - **Compactness**: Always compact (finite sets are compact)
/// - **Closure**: The singleton is its own closure
/// - **Interior**: The interior is always empty
/// - **Boundary**: The singleton is its own boundary
///
/// These properties are preserved regardless of the scalar type used and are guaranteed
/// by the underlying implementations of the [`IntervalTrait`].
#[derive(Debug, Clone, PartialEq, Eq, IntoInner, Serialize, Deserialize)]
#[serde(bound(deserialize = "RealType: for<'a> Deserialize<'a>"))]
pub struct IntervalSingleton<RealType: RealScalar>(RealType);

impl<RealType: RealScalar> IntervalSingleton<RealType> {
    /// Get the value of the singleton.
    #[inline(always)]
    pub fn value(&self) -> &RealType {
        &self.0
    }
}

impl<RealType: RealScalar> New for IntervalSingleton<RealType> {
    fn new(value: RealType) -> Self {
        debug_assert!(value.is_finite(), "The input value {value} is not finite!");
        Self(value)
    }
}

impl<RealType: RealScalar> IntervalOperations for IntervalSingleton<RealType> {}

impl<RealType: RealScalar> Contains for IntervalSingleton<RealType> {
    /// Returns `true` if the 1D point/coordinate `x` is equal to the value of the current **singleton**.
    #[inline(always)]
    fn contains_point(&self, x: &RealType) -> bool {
        x == &self.0
    }
}

impl<RealType: RealScalar> IntervalTrait for IntervalSingleton<RealType> {}
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
impl<RealType: RealScalar> IntervalBoundsRuntime for IntervalSingleton<RealType> {
    type RealType = RealType;

    #[inline(always)]
    fn lower_bound_runtime(&self) -> Option<LowerBoundRuntime<RealType>> {
        Some(IntervalBoundRuntime::Closed(LowerBoundClosed::new(
            self.value().clone(),
        )))
    }

    #[inline(always)]
    fn upper_bound_runtime(&self) -> Option<UpperBoundRuntime<RealType>> {
        Some(IntervalBoundRuntime::Closed(UpperBoundClosed::new(
            self.value().clone(),
        )))
    }
}

impl<RealType: RealScalar> IntervalBoundsRuntime for IntervalFinitePositiveLength<RealType> {
    type RealType = RealType;

    #[inline(always)]
    fn lower_bound_runtime(&self) -> Option<LowerBoundRuntime<RealType>> {
        match self {
            IntervalFinitePositiveLength::Open(interval) => interval.lower_bound_runtime(),
            IntervalFinitePositiveLength::Closed(interval) => interval.lower_bound_runtime(),
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.lower_bound_runtime()
            }
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.lower_bound_runtime()
            }
        }
    }

    #[inline(always)]
    fn upper_bound_runtime(&self) -> Option<UpperBoundRuntime<RealType>> {
        match self {
            IntervalFinitePositiveLength::Open(interval) => interval.upper_bound_runtime(),
            IntervalFinitePositiveLength::Closed(interval) => interval.upper_bound_runtime(),
            IntervalFinitePositiveLength::LowerClosedUpperOpen(interval) => {
                interval.upper_bound_runtime()
            }
            IntervalFinitePositiveLength::LowerOpenUpperClosed(interval) => {
                interval.upper_bound_runtime()
            }
        }
    }
}
//------------------------------------------------------------------------------------------------

/// Enum representing the different types of sub-intervals that can appear in a partition of a finite interval with positive length.
///
/// The [`SubIntervalInPartition`] enum encodes the **position-dependent boundary semantics** that arise
/// when partitioning an interval into adjacent, non-overlapping sub-intervals. The type of each sub-interval
/// depends on its position (first, middle, or last) and the boundary semantics of the original domain being partitioned.
///
/// ## Mathematical Foundation
///
/// When partitioning an interval domain `D` with partition points `{p₀, p₁, ..., pₙ}`, the resulting
/// sub-intervals must satisfy two critical requirements:
/// 1. **Complete coverage**: `⋃ᵢ Iᵢ = D` (union of all sub-intervals equals the original domain)
/// 2. **Non-overlapping**: `Iᵢ ∩ Iⱼ = ∅` for `i ≠ j` (disjoint interiors, at most boundary point overlap)
///
/// To achieve this while preserving the domain's boundary semantics, different sub-intervals require
/// different boundary types based on their position in the partition.
///
/// ## Core Design Philosophy
///
/// ### Why Four Variants?
/// The library uses a **partition strategy** that ensures:
/// - **Uniqueness**: Every point in the domain belongs to exactly one sub-interval
/// - **Boundary preservation**: Domain boundary semantics are maintained
/// - **Consistency**: Internal partition points never create ambiguity
///
/// This is achieved by following a specific rule: **partition points belong to the left sub-interval**
/// (except at the domain boundaries where the domain's own semantics apply).
///
/// ### The Partition Strategy
///
/// Given partition points `{p₀, p₁, ..., pₙ}` and domain `D`:
///
/// | Position | Sub-Interval Type | Mathematical Notation | Reason |
/// |----------|-------------------|----------------------|---------|
/// | **Only** (n=1) | Matches domain exactly | Depends on `D` | Single interval = domain itself |
/// | **First** (i=0) | Depends on domain | `[p₀, p₁)` for `[a,b]`/`[a,b)`; `(p₀, p₁]` for `(a,b)`/`(a,b]` | Inherits domain's boundary type |
/// | **Middle** (0<i<n-1) | Determined by `MiddleIntervalInPartition` | `[pᵢ, pᵢ₊₁)` for `[a,b]`/`[a,b)`; `(pᵢ, pᵢ₊₁]` for `(a,b)`/`(a,b]` | Avoids ambiguity at partition points |
/// | **Last** (i=n-1) | Depends on domain | `[pₙ₋₁, pₙ]` for `[a,b]`; `(pₙ₋₁, pₙ)` for `(a,b)`; `[pₙ₋₁, pₙ)` for `[a,b)`; `(pₙ₋₁, pₙ]` for `(a,b]` | Inherits domain's boundary type |
///
/// ### Why This Strategy?
/// Consider a closed domain `[0, 2]` partitioned at `{0, 1, 2}`:
///
/// **Strategy 1: All closed `[pᵢ, pᵢ₊₁]`** ❌
/// - Sub-intervals: `[0, 1]`, `[1, 2]`
/// - Problem: Point `1` belongs to **both** intervals (ambiguity!)
/// - Violates uniqueness requirement
///
/// **Strategy 2: All open `(pᵢ, pᵢ₊₁)`** ❌
/// - Sub-intervals: `(0, 1)`, `(1, 2)`
/// - Problem: Partition points `0, 1, 2` belong to **no** interval (gaps!)
/// - Violates complete coverage
///
/// **Strategy 3: Our approach (for `[a, b]` domain)** ✅
/// - Sub-intervals: `[0, 1)`, `[1, 2]`
/// - Result: Point `0 ∈ [0, 1)`, point `1 ∈ [1, 2]`, point `2 ∈ [1, 2]`
/// - Every point belongs to exactly one interval!
///
/// ## Variants
///
/// ### [`SubIntervalInPartition::Single`]
///
/// **Used when**: The partition consists of a **single interval** (n=1, meaning only 2 partition points).
///
/// **Semantics**: The sub-interval is **identical** to the original domain.
///
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// // Single interval partition of closed domain [0, 1]
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0]);
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
///
/// // The only interval matches the domain exactly
/// let only_interval = grid.interval(&IntervalId::new(0));
/// // Type: SubIntervalInPartition::Single(IntervalClosed<f64>)
/// // Result: [0, 1] - same as domain
///
/// // Single interval partition of open domain (0, 1)
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0]);
/// let grid = Grid1D::<IntervalOpen<f64>>::try_from_sorted(coords).unwrap();
///
/// // The only interval matches the domain exactly
/// let only_interval = grid.interval(&IntervalId::new(0));
/// // Type: SubIntervalInPartition::Single(IntervalOpen<f64>)
/// // Result: (0, 1) - same as domain
/// ```
///
/// **Mathematical property**: `Only(D) = D` for any domain `D`.
///
/// ### [`SubIntervalInPartition::First`]
///
/// **Used when**: This is the **leftmost** sub-interval in a multi-interval partition.
///
/// **Semantics**: The **lower bound** inherits the domain's lower boundary semantics.
/// The **inner (right) bound** type follows the domain's partition strategy.
///
/// **Type construction** — depends on the domain:
///
/// | Domain | First Interval Type | Notation |
/// |--------|---------------------|----------|
/// | `[a, b]` | `IntervalLowerClosedUpperOpen` | `[p₀, p₁)` |
/// | `(a, b)` | `IntervalLowerOpenUpperClosed` | `(p₀, p₁]` |
/// | `[a, b)` | `IntervalLowerClosedUpperOpen` | `[p₀, p₁)` |
/// | `(a, b]` | `IntervalLowerOpenUpperClosed` | `(p₀, p₁]` |
///
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// // Closed domain [0, 2] partitioned at {0, 1, 2}
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0]);
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
///
/// let first_interval = grid.interval(&IntervalId::new(0));
/// // Type: SubIntervalInPartition::First(IntervalClosed<f64>)
/// // Result: [0, 1) - closed lower bound - open upper bound
/// assert!(first_interval.contains_point(&0.0)); // ✓ includes lower bound
/// assert!(!first_interval.contains_point(&1.0)); // ✗ excludes upper partition point
///
/// // Open domain (0, 2) partitioned at {0, 1, 2}
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0]);
/// let grid = Grid1D::<IntervalOpen<f64>>::try_from_sorted(coords).unwrap();
///
/// let first_interval = grid.interval(&IntervalId::new(0));
/// // Type: SubIntervalInPartition::First(IntervalLowerOpenUpperClosed<f64>)
/// // Result: (0, 1] - open lower bound inherited from domain
/// assert!(!first_interval.contains_point(&0.0)); // ✗ excludes lower bound
/// assert!(first_interval.contains_point(&1.0));  // ✓ includes upper partition point
/// ```
///
/// **Inner bound convention**: For `[a, b]` and `[a, b)` domains the inner bound is open so
/// that `p₁` belongs to the **next** sub-interval. For `(a, b)` and `(a, b]` domains the inner
/// bound is closed so that `p₁` belongs to the First interval itself. Either way every partition
/// point belongs to exactly one sub-interval.
///
/// ### [`SubIntervalInPartition::Middle`]
///
/// **Used when**: This is an **interior** sub-interval (neither first nor last) in the partition.
///
/// **Semantics**: The concrete type is [`PartitionDomain1D::MiddleIntervalInPartition`](Grid1DIntervalBuilder::MiddleIntervalInPartition),
/// determined by the domain's [`Grid1DIntervalBuilder`] implementation.
///
/// | Domain | Middle Interval Type | Notation |
/// |--------|----------------------|----------|
/// | `[a, b]` | `IntervalLowerClosedUpperOpen` | `[pₖ, pₖ₊₁)` |
/// | `(a, b)` | `IntervalLowerOpenUpperClosed` | `(pₖ, pₖ₊₁]` |
/// | `[a, b)` | `IntervalLowerClosedUpperOpen` | `[pₖ, pₖ₊₁)` |
/// | `(a, b]` | `IntervalLowerOpenUpperClosed` | `(pₖ, pₖ₊₁]` |
///
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// // Domain [0, 3] partitioned at {0, 1, 2, 3}
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0, 3.0]);
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
///
/// // Middle intervals: [1, 2) and [2, 3) (last is not middle!)
/// let middle_interval = grid.interval(&IntervalId::new(1));
/// // Type: SubIntervalInPartition::Middle(IntervalLowerClosedUpperOpen<f64>)
/// // Result: [1, 2)
/// assert!(middle_interval.contains_point(&1.0)); // ✓ includes left partition point
/// assert!(!middle_interval.contains_point(&2.0));  // ✗ excludes right partition point
///
/// // This prevents point 1.0 from belonging to both intervals [0, 1) and [1, 2)
/// let first_interval = grid.interval(&IntervalId::new(0));
/// assert!(!first_interval.contains_point(&1.0));   // Point 1 belongs to middle interval
/// ```
///
/// **Partition convention**: The middle interval's boundary type determines which sub-interval
/// each interior partition point belongs to:
/// - `[pₖ, pₖ₊₁)` (used by `[a,b]`, `[a,b)`): `pₖ` belongs to the **right** sub-interval
/// - `(pₖ, pₖ₊₁]` (used by `(a,b)`, `(a,b]`): `pₖ` belongs to the **left** sub-interval
///
/// Either convention guarantees that every partition point belongs to exactly one sub-interval.
/// Custom `Grid1DIntervalBuilder` implementations may use a different middle interval type,
/// provided the same uniqueness guarantee holds.
///
/// ### [`SubIntervalInPartition::Last`]
///
/// **Used when**: This is the **rightmost** sub-interval in a multi-interval partition.
///
/// **Semantics**: The **upper bound** inherits the domain's upper boundary semantics.
/// The **inner (left) bound** type follows the domain's partition strategy.
///
/// **Type construction** — depends on the domain:
///
/// | Domain | Last Interval Type | Notation |
/// |--------|-------------------|----------|
/// | `[a, b]` | `IntervalClosed` | `[pₙ₋₁, pₙ]` |
/// | `(a, b)` | `IntervalOpen` | `(pₙ₋₁, pₙ)` |
/// | `[a, b)` | `IntervalLowerClosedUpperOpen` | `[pₙ₋₁, pₙ)` |
/// | `(a, b]` | `IntervalLowerOpenUpperClosed` | `(pₙ₋₁, pₙ]` |
///
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// // Closed domain [0, 2] partitioned at {0, 1, 2}
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0]);
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
///
/// let last_interval = grid.interval(&IntervalId::new(1));
/// // Type: SubIntervalInPartition::Last(IntervalClosed<f64>)
/// // Result: [1, 2] - closed upper bound inherited from domain
/// assert!(last_interval.contains_point(&1.0)); // ✓ includes left partition point
/// assert!(last_interval.contains_point(&2.0));  // ✓ includes upper bound from domain
///
/// // Semi-open domain [0, 2) partitioned at {0, 1, 2}
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0]);
/// let grid = Grid1D::<IntervalLowerClosedUpperOpen<f64>>::try_from_sorted(coords).unwrap();
///
/// let last_interval = grid.interval(&IntervalId::new(1));
/// // Type: SubIntervalInPartition::Last(IntervalLowerClosedUpperOpen<f64>)
/// // Result: [1, 2) - open upper bound
/// assert!(last_interval.contains_point(&1.0)); // ✓ includes lower bound
/// assert!(!last_interval.contains_point(&2.0)); // ✗ excludes upper bound from domain
/// ```
///
/// **Inner bound convention**: For `[a, b]` and `[a, b)` domains the inner bound is closed so
/// that `pₙ₋₁` belongs to the Last interval itself. For `(a, b)` and `(a, b]` domains the inner
/// bound is open so that `pₙ₋₁` belongs to the previous sub-interval. Either way the
/// partition's uniqueness property is maintained.
///
/// ## Comprehensive Domain Examples
///
/// ### Closed Domain `[a, b]`
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0, 3.0]);
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(coords).unwrap();
///
/// // Partition result:
/// // - First:  [0, 1) (closed-open)
/// // - Middle: [1, 2) (closed-open)
/// // - Last:   [2, 3] (closed-closed)
///
/// let interval_0 = grid.interval(&IntervalId::new(0));
/// let interval_1 = grid.interval(&IntervalId::new(1));
/// let interval_2 = grid.interval(&IntervalId::new(2));
///
/// // Point 0 belongs only to first interval
/// use grid1d::intervals::{Contains, IntervalTrait};
/// assert!(interval_0.contains_point(&0.0));
/// assert!(!interval_1.contains_point(&0.0));
/// assert!(!interval_2.contains_point(&0.0));
///
/// // Point 1 belongs only to second interval (left boundary)
/// assert!(!interval_0.contains_point(&1.0));
/// assert!(interval_1.contains_point(&1.0));
/// assert!(!interval_2.contains_point(&1.0));
///
/// // Point 2 belongs only to last interval (left boundary)
/// assert!(!interval_0.contains_point(&2.0));
/// assert!(!interval_1.contains_point(&2.0));
/// assert!(interval_2.contains_point(&2.0));
///
/// // Point 3 belongs only to last interval (right boundary)
/// assert!(!interval_0.contains_point(&3.0));
/// assert!(!interval_1.contains_point(&3.0));
/// assert!(interval_2.contains_point(&3.0));
/// ```
///
/// ### Open Domain `(a, b)`
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0, 3.0]);
/// let grid = Grid1D::<IntervalOpen<f64>>::try_from_sorted(coords).unwrap();
///
/// // Partition result:
/// // - First:  (0, 1] (open-closed, inherits open lower from domain)
/// // - Middle: (1, 2] (open-closed)
/// // - Last:   (2, 3) (open-open, inherits open upper from domain)
///
/// let interval_0 = grid.interval(&IntervalId::new(0));
/// let interval_2 = grid.interval(&IntervalId::new(2));
///
/// // Point 0 belongs to no interval (domain excludes it)
/// assert!(!interval_0.contains_point(&0.0));
///
/// // Point 3 belongs to no interval (domain excludes it)
/// assert!(!interval_2.contains_point(&3.0));
/// ```
///
/// ### Semi-Open Domain `[a, b)`
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// let coords = SortedSet::from_unsorted(vec![0.0, 1.0, 2.0, 3.0]);
/// let grid = Grid1D::<IntervalLowerClosedUpperOpen<f64>>::try_from_sorted(coords).unwrap();
///
/// // Partition result (all intervals share the domain's type [x, y)):
/// // - First:  [0, 1) (closed-open, inherits closed lower from domain)
/// // - Middle: [1, 2) (closed-open, partition point belongs to right interval)
/// // - Last:   [2, 3) (closed-open, inherits open upper from domain)
///
/// let interval_0 = grid.interval(&IntervalId::new(0));
/// let interval_2 = grid.interval(&IntervalId::new(2));
///
/// // Point 0 belongs to first interval (domain includes it)
/// assert!(interval_0.contains_point(&0.0));
///
/// // Point 3 belongs to no interval (domain excludes it)
/// assert!(!interval_2.contains_point(&3.0));
/// ```
///
/// ## Generic Type Parameter
///
/// ```rust,ignore
/// pub enum SubIntervalInPartition<PartitionDomain1D: Grid1DIntervalBuilder> {
///     Only(PartitionDomain1D),
///     First(PartitionDomain1D::FirstIntervalInPartition),
///     Middle(PartitionDomain1D::MiddleIntervalInPartition),
///     Last(PartitionDomain1D::LastIntervalInPartition),
/// }
/// ```
///
/// **`PartitionDomain1D`**: The type of the original domain being partitioned.
/// - Must implement [`Grid1DIntervalBuilder`] trait
/// - Determines the scalar type (`RealType`) used for coordinates
/// - Defines boundary semantics for first and last sub-intervals
///
/// ## Mathematical Guarantees
///
/// The [`SubIntervalInPartition`] enum, combined with the partition strategy, guarantees:
///
/// ### 1. Complete Coverage
/// ```text
/// ∀x ∈ D, ∃i : x ∈ Iᵢ
/// ```
/// Every point in the domain belongs to at least one sub-interval.
///
/// ### 2. Uniqueness (No Overlap)
/// ```text
/// ∀x ∈ D, |{i : x ∈ Iᵢ}| = 1
/// ```
/// Every point in the domain belongs to **exactly one** sub-interval.
///
/// ### 3. Exact Reconstruction
/// ```text
/// ⋃ᵢ Iᵢ = D
/// ```
/// The union of all sub-intervals exactly equals the original domain.
///
/// ### 4. Partition Point Assignment
/// ```text
/// ∀i ∈ (0, n-1) : pᵢ ∈ Iᵢ₋₁  (partition points belong to the left interval)
/// ```
/// Internal partition points always belong to the **left** sub-interval.
///
/// ### 5. Boundary Preservation
/// ```text
/// lower_bound(I₀) = lower_bound(D)  (semantics preserved)
/// upper_bound(Iₙ₋₁) = upper_bound(D)  (semantics preserved)
/// ```
/// Domain boundary semantics are preserved in the first and last sub-intervals.
///
/// ## Trait Implementations
///
/// The enum implements all standard interval operations through its variants:
///
/// - **[`IntervalTrait`]**: Full interval arithmetic and containment testing
/// - **[`IntervalFinitePositiveLengthTrait`]**: Length, midpoint, bounds access
/// - **[`GetLowerBoundValue`], [`GetUpperBoundValue`]**: Boundary value access
/// - **[`Contains`]**: Point containment testing
/// - **[`IntervalBoundsRuntime`]**: Runtime boundary type inspection
///
/// All operations delegate to the appropriate variant's implementation.
///
/// ## Performance Characteristics
///
/// | Operation | Time Complexity | Notes |
/// |-----------|-----------------|-------|
/// | **Variant matching** | O(1) | Single branch prediction |
/// | **Point containment** | O(1) | Delegates to concrete interval type |
/// | **Bounds access** | O(1) | Direct field access |
/// | **Length calculation** | O(1) | Cached or computed once |
/// | **Intersection** | O(1) | Delegates to interval operations |
///
/// ## Integration with Grid System
///
/// The [`SubIntervalInPartition`] enum is the **return type** for the
/// [`HasIntervals::interval`](crate::HasIntervals::interval) method, which is implemented by:
/// - [`Grid1DUniform`](crate::Grid1DUniform): Uniform spacing grids
/// - [`Grid1DNonUniform`](crate::Grid1DNonUniform): Non-uniform spacing grids
/// - [`Grid1D`](crate::Grid1D): Unified grid interface
/// - [`Grid1DUnion`](crate::Grid1DUnion): Union of multiple grids
///
/// ## Best Practices
///
/// ### Pattern Matching
/// ```rust
/// use grid1d::{*, intervals::*, scalars::*};
/// use sorted_vec::partial::SortedSet;
///
/// let grid = Grid1D::<IntervalClosed<f64>>::try_from_sorted(
///     SortedSet::from_unsorted(vec![0.0, 1.0, 2.0])
/// ).unwrap();
///
/// let interval = grid.interval(&IntervalId::new(0));
///
/// match interval {
///     SubIntervalInPartition::Single(domain) => {
///         println!("Single interval partition matching domain");
///     }
///     SubIntervalInPartition::First(first) => {
///         println!("First interval with inherited lower bound");
///     }
///     SubIntervalInPartition::Middle(middle) => {
///         println!("Middle interval (type depends on domain: `(a,b]` or `[a,b)` for built-in types)");
///     }
///     SubIntervalInPartition::Last(last) => {
///         println!("Last interval with inherited upper bound");
///     }
/// }
/// ```
///
/// ### Generic Operations
/// ```rust
/// use grid1d::{Grid1DIntervalBuilder, Grid1D, Grid1DTrait, intervals::*, scalars::IntervalId};
/// use sorted_vec::partial::SortedSet;
/// use num_valid::RealScalar;
///
/// fn analyze_sub_interval<T, D>(interval: &SubIntervalInPartition<D>) -> String
/// where
///     T: RealScalar,
///     D: Grid1DIntervalBuilder<RealType = T>,
///     SubIntervalInPartition<D>: IntervalFinitePositiveLengthTrait<RealType = T>,
/// {
///     format!(
///         "Interval: [{}, {}], length: {}",
///         interval.lower_bound_value(),
///         interval.upper_bound_value(),
///         interval.length()
///     )
/// }
/// ```
///
/// ## See Also
///
/// - [`Grid1DIntervalBuilder`]: Trait for constructing sub-intervals based on domain type
/// - [`Grid1DTrait`](crate::Grid1DTrait): Core trait that uses this enum for interval access
/// - [`IntervalTrait`]: Common interval operations
/// - [`Grid1D`](crate::Grid1D), [`Grid1DUniform`](crate::Grid1DUniform), [`Grid1DNonUniform`](crate::Grid1DNonUniform): Grid types that produce these sub-intervals
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
#[serde(bound(deserialize = "PartitionDomain1D: for<'a> Deserialize<'a>"))]
pub enum SubIntervalInPartition<PartitionDomain1D: Grid1DIntervalBuilder> {
    /// A partition with a single interval covering the entire domain.
    Single(PartitionDomain1D),

    /// The first interval in a multi-interval partition.
    First(PartitionDomain1D::FirstIntervalInPartition),

    /// A middle interval in a multi-interval partition.
    /// The concrete type is [`Grid1DIntervalBuilder::MiddleIntervalInPartition`];
    /// for `[a, b]` and `[a, b)` domains this is `[a, b)`, for `(a, b)` and `(a, b]` it is `(a, b]`.
    Middle(PartitionDomain1D::MiddleIntervalInPartition),

    /// The last interval in a multi-interval partition.
    Last(PartitionDomain1D::LastIntervalInPartition),
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> IntervalBoundsRuntime
    for SubIntervalInPartition<PartitionDomain1D>
{
    type RealType = PartitionDomain1D::RealType;

    #[inline]
    fn lower_bound_runtime(&self) -> Option<LowerBoundRuntime<PartitionDomain1D::RealType>> {
        match self {
            SubIntervalInPartition::Single(interval) => interval.lower_bound_runtime(),
            SubIntervalInPartition::First(interval) => interval.lower_bound_runtime(),
            SubIntervalInPartition::Middle(interval) => interval.lower_bound_runtime(),
            SubIntervalInPartition::Last(interval) => interval.lower_bound_runtime(),
        }
    }

    #[inline]
    fn upper_bound_runtime(&self) -> Option<UpperBoundRuntime<PartitionDomain1D::RealType>> {
        match self {
            SubIntervalInPartition::Single(interval) => interval.upper_bound_runtime(),
            SubIntervalInPartition::First(interval) => interval.upper_bound_runtime(),
            SubIntervalInPartition::Middle(interval) => interval.upper_bound_runtime(),
            SubIntervalInPartition::Last(interval) => interval.upper_bound_runtime(),
        }
    }
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> IntervalOperations
    for SubIntervalInPartition<PartitionDomain1D>
{
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> Contains
    for SubIntervalInPartition<PartitionDomain1D>
{
    #[inline(always)]
    fn contains_point(&self, x: &PartitionDomain1D::RealType) -> bool {
        match self {
            SubIntervalInPartition::Single(interval) => interval.contains_point(x),
            SubIntervalInPartition::First(interval) => interval.contains_point(x),
            SubIntervalInPartition::Middle(interval) => interval.contains_point(x),
            SubIntervalInPartition::Last(interval) => interval.contains_point(x),
        }
    }
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> IntervalTrait
    for SubIntervalInPartition<PartitionDomain1D>
{
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> GetLowerBoundValue
    for SubIntervalInPartition<PartitionDomain1D>
{
    type LowerBoundValue = PartitionDomain1D::RealType;

    #[inline(always)]
    fn lower_bound_value(&self) -> &Self::LowerBoundValue {
        match self {
            SubIntervalInPartition::Single(interval) => interval.lower_bound_value(),
            SubIntervalInPartition::First(interval) => interval.lower_bound_value(),
            SubIntervalInPartition::Middle(interval) => interval.lower_bound_value(),
            SubIntervalInPartition::Last(interval) => interval.lower_bound_value(),
        }
    }

    #[inline(always)]
    fn is_lower_bound_closed(&self) -> bool {
        match self {
            SubIntervalInPartition::Single(interval) => interval.is_lower_bound_closed(),
            SubIntervalInPartition::First(interval) => interval.is_lower_bound_closed(),
            SubIntervalInPartition::Middle(interval) => interval.is_lower_bound_closed(),
            SubIntervalInPartition::Last(interval) => interval.is_lower_bound_closed(),
        }
    }
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> GetUpperBoundValue
    for SubIntervalInPartition<PartitionDomain1D>
{
    type UpperBoundValue = PartitionDomain1D::RealType;

    #[inline(always)]
    fn upper_bound_value(&self) -> &Self::UpperBoundValue {
        match self {
            SubIntervalInPartition::Single(interval) => interval.upper_bound_value(),
            SubIntervalInPartition::First(interval) => interval.upper_bound_value(),
            SubIntervalInPartition::Middle(interval) => interval.upper_bound_value(),
            SubIntervalInPartition::Last(interval) => interval.upper_bound_value(),
        }
    }

    /// Check if the upper bound of this sub-interval is closed (inclusive).
    ///
    /// Returns `true` if the upper boundary point is included in the interval,
    /// `false` if it is excluded.
    ///
    /// # Partition Semantics
    ///
    /// The upper bound closure depends on the variant and the partition domain type:
    ///
    /// - **`SubIntervalInPartition::Single`**: Inherits the upper bound type from the partition domain
    ///   - Closed domain `[a, b]` → closed upper bound
    ///   - Open domain `(a, b)` → open upper bound
    ///   - Half-open `[a, b)` → open upper bound
    ///
    /// - **`SubIntervalInPartition::First`**: For partitions with multiple intervals, first intervals typically
    ///   have closed upper bounds to ensure gap-free coverage. The exact type depends
    ///   on the partition domain's implementation of [`Grid1DIntervalBuilder`].
    ///
    /// - **`SubIntervalInPartition::Middle`**: Determined by [`Grid1DIntervalBuilder::MiddleIntervalInPartition`].
    ///   - For `[a, b]` and `[a, b)` domains this is `[a, b)` (open upper bound).
    ///   - For `(a, b)` and `(a, b]` domains this is `(a, b]` (closed upper bound).
    ///
    /// - **`SubIntervalInPartition::Last`**: Inherits the upper bound type from the partition domain
    ///   - Closed domain `[a, b]` → closed upper bound
    ///   - Open domain `(a, b)` → open upper bound
    ///   - Half-open `[a, b)` → open upper bound
    ///   - Half-open `(a, b]` → closed upper bound
    ///
    /// # Examples
    ///
    /// ```
    /// use grid1d::intervals::*;
    ///
    /// // Only variant with closed domain [0, 10]
    /// let only_closed: SubIntervalInPartition<IntervalClosed<f64>> =
    ///     SubIntervalInPartition::Single(IntervalClosed::new(0.0, 10.0));
    /// assert!(only_closed.is_upper_bound_closed());  // Closed upper bound
    ///
    /// // Only variant with open domain (0, 10)
    /// let only_open: SubIntervalInPartition<IntervalOpen<f64>> =
    ///     SubIntervalInPartition::Single(IntervalOpen::new(0.0, 10.0));
    /// assert!(!only_open.is_upper_bound_closed());   // Open upper bound
    ///
    /// // Middle variant: for IntervalClosed domain, MiddleIntervalInPartition = IntervalLowerClosedUpperOpen
    /// let middle: SubIntervalInPartition<IntervalClosed<f64>> =
    ///     SubIntervalInPartition::Middle(IntervalLowerClosedUpperOpen::new(3.0, 7.0));
    /// assert!(!middle.is_upper_bound_closed());
    ///
    /// // Last variant with half-open domain [0, 10)
    /// let last_half_open: SubIntervalInPartition<IntervalLowerClosedUpperOpen<f64>> =
    ///     SubIntervalInPartition::Last(IntervalLowerClosedUpperOpen::new(7.0, 10.0));
    /// assert!(!last_half_open.is_upper_bound_closed()); // Open upper bound from domain
    /// ```
    ///
    /// # See Also
    ///
    /// - [`is_upper_bound_open`](GetUpperBoundValue::is_upper_bound_open) - Logical complement
    /// - [`is_lower_bound_closed`](GetLowerBoundValue::is_lower_bound_closed) - Lower bound closure
    /// - [`upper_bound_value`](GetUpperBoundValue::upper_bound_value) - Get the upper bound value
    #[inline(always)]
    fn is_upper_bound_closed(&self) -> bool {
        match self {
            SubIntervalInPartition::Single(interval) => interval.is_upper_bound_closed(),
            SubIntervalInPartition::First(interval) => interval.is_upper_bound_closed(),
            SubIntervalInPartition::Middle(interval) => interval.is_upper_bound_closed(),
            SubIntervalInPartition::Last(interval) => interval.is_upper_bound_closed(),
        }
    }

    /// Check if the upper bound of this sub-interval is open (exclusive).
    ///
    /// Returns `true` if the upper boundary point is excluded from the interval,
    /// `false` if it is included. This is the logical complement of
    /// [`is_upper_bound_closed`](GetUpperBoundValue::is_upper_bound_closed).
    ///
    /// # Partition Semantics
    ///
    /// The upper bound openness depends on the variant and partition domain type:
    ///
    /// - **`SubIntervalInPartition::Single`**: Inherits from the partition domain
    ///   - Closed domain `[a, b]` → closed upper (returns `false`)
    ///   - Open domain `(a, b)` → open upper (returns `true`)
    ///   - Half-open `[a, b)` → open upper (returns `true`)
    ///   - Half-open `(a, b]` → closed upper (returns `false`)
    ///
    /// - **`SubIntervalInPartition::First`**: Typically has closed upper bounds for gap-free partitions,
    ///   so usually returns `false`. The exact behavior depends on the
    ///   [`Grid1DIntervalBuilder`] implementation.
    ///
    /// - **`SubIntervalInPartition::Middle`**: Returns `true` for `[a, b]` and `[a, b)` domains (upper bound open).
    ///   Returns `false` for `(a, b)` and `(a, b]` domains (upper bound closed). Custom
    ///   `SubIntervalInPartition::Middle` types may behave differently.
    ///
    /// - **`SubIntervalInPartition::Last`**: Inherits from the partition domain
    ///   - Closed domain → closed upper (returns `false`)
    ///   - Open domain → open upper (returns `true`)
    ///   - Half-open `[a, b)` → open upper (returns `true`)
    ///   - Half-open `(a, b]` → closed upper (returns `false
    ///
    /// # Examples
    ///
    /// ```
    /// use grid1d::intervals::*;
    ///
    /// // Only variant with half-open domain [0, 10)
    /// let only_half_open: SubIntervalInPartition<IntervalLowerClosedUpperOpen<f64>> =
    ///     SubIntervalInPartition::Single(IntervalLowerClosedUpperOpen::new(0.0, 10.0));
    /// assert!(only_half_open.is_upper_bound_open());  // Open upper bound
    ///
    /// // Middle variant: for IntervalClosed domain, MiddleIntervalInPartition = IntervalLowerClosedUpperOpen
    /// let middle: SubIntervalInPartition<IntervalClosed<f64>> =
    ///     SubIntervalInPartition::Middle(IntervalLowerClosedUpperOpen::new(3.0, 7.0));
    /// assert!(middle.is_upper_bound_open());
    ///
    /// // Verify complement relationship
    /// assert_eq!(
    ///     middle.is_upper_bound_open(),
    ///     !middle.is_upper_bound_closed()
    /// );
    /// ```
    ///
    /// # Implementation Note
    ///
    /// This method has a default implementation in [`GetUpperBoundValue`] as
    /// `!self.is_upper_bound_closed()`. It delegates to the underlying interval's
    /// implementation through pattern matching.
    ///
    /// # See Also
    ///
    /// - [`is_upper_bound_closed`](GetUpperBoundValue::is_upper_bound_closed) - Logical complement
    /// - [`is_lower_bound_open`](GetLowerBoundValue::is_lower_bound_open) - Lower bound openness
    #[inline(always)]
    fn is_upper_bound_open(&self) -> bool {
        !self.is_upper_bound_closed()
    }
}

impl<PartitionDomain1D: Grid1DIntervalBuilder> IntervalFinitePositiveLengthTrait
    for SubIntervalInPartition<PartitionDomain1D>
{
    /// Return the inner bounds of the interval as a pair in which the first entry of the pair
    /// is the lower bound and the second entry of the pair is the upper bound.
    ///
    /// # Note
    /// This function consumes `Self` (no cloning of the bounds is needed).
    #[inline(always)]
    fn into_bounds_pair(self) -> (Self::RealType, Self::RealType) {
        match self {
            SubIntervalInPartition::Single(interval) => interval.into_bounds_pair(),
            SubIntervalInPartition::First(interval) => interval.into_bounds_pair(),
            SubIntervalInPartition::Middle(interval) => interval.into_bounds_pair(),
            SubIntervalInPartition::Last(interval) => interval.into_bounds_pair(),
        }
    }

    fn translate(&self, amount: Self::RealType) -> Self {
        match self {
            SubIntervalInPartition::Single(interval) => {
                SubIntervalInPartition::Single(interval.translate(amount))
            }
            SubIntervalInPartition::First(interval) => {
                SubIntervalInPartition::First(interval.translate(amount))
            }
            SubIntervalInPartition::Middle(interval) => {
                SubIntervalInPartition::Middle(interval.translate(amount))
            }
            SubIntervalInPartition::Last(interval) => {
                SubIntervalInPartition::Last(interval.translate(amount))
            }
        }
    }

    fn scale(&self, factor: &PositiveRealScalar<Self::RealType>) -> Self {
        match self {
            SubIntervalInPartition::Single(interval) => {
                SubIntervalInPartition::Single(interval.scale(factor))
            }
            SubIntervalInPartition::First(interval) => {
                SubIntervalInPartition::First(interval.scale(factor))
            }
            SubIntervalInPartition::Middle(interval) => {
                SubIntervalInPartition::Middle(interval.scale(factor))
            }
            SubIntervalInPartition::Last(interval) => {
                SubIntervalInPartition::Last(interval.scale(factor))
            }
        }
    }
    fn expand(&self, amount: &PositiveRealScalar<Self::RealType>) -> Self {
        match self {
            SubIntervalInPartition::Single(interval) => {
                SubIntervalInPartition::Single(interval.expand(amount))
            }
            SubIntervalInPartition::First(interval) => {
                SubIntervalInPartition::First(interval.expand(amount))
            }
            SubIntervalInPartition::Middle(interval) => {
                SubIntervalInPartition::Middle(interval.expand(amount))
            }
            SubIntervalInPartition::Last(interval) => {
                SubIntervalInPartition::Last(interval.expand(amount))
            }
        }
    }
}

// IntervalFinitePositiveLength hull implementation
impl<RealType: RealScalar, I2: IntervalFinitePositiveLengthTrait<RealType = RealType>>
    IntervalHull<I2> for IntervalFinitePositiveLength<RealType>
{
    type Output = IntervalFinitePositiveLength<RealType>;

    #[inline]
    fn hull(&self, other: &I2) -> Self::Output {
        compute_hull_bounded_intervals(self, other)
    }
}