grid1d 0.4.0

//! Core traits for grid1d.
//!
//! This module contains all the fundamental trait definitions used throughout the `grid1d` library.
//! These traits define the core abstractions for domains, coordinates, interval partitions,
//! and point transformations.
//!
//! # Trait Overview
//!
//! | Trait | Purpose | Key Methods |
//! |-------|---------|-------------|
//! | [`HasDomain1D`] | Domain access | `domain()` |
//! | [`HasCoords1D`] | Coordinate access | `coords()`, `num_points()` |
//! | [`BuildIntervalInPartition`] | Sub-interval construction rules | `build_first_interval()`, `build_last_interval()` |
//! | [`IntervalPartition`] | Main partition interface | `interval()`, `find_interval_id_of_point()`, `refine()` |
//! | [`TransformedPoints1D`] | Point transformations | `points_param_interval()`, `points_phys_interval()` |
//!
//! # Trait Hierarchy
//!
//! ```text
//! IntervalPartition
//!     ├── HasCoords1D  (coordinate access)
//!     └── HasDomain1D  (domain access)
//!             └── Domain1D: BuildIntervalInPartition
//!
//! BuildIntervalInPartition
//!     └── IntervalFinitePositiveLengthTrait
//!
//! TransformedPoints1D (independent)
//! ```
//!
//! # Usage
//!
//! Most users will interact with these traits through the concrete types like
//! [`Grid1D`](crate::Grid1D), [`Grid1DUniform`](crate::Grid1DUniform), and
//! [`Grid1DNonUniform`](crate::Grid1DNonUniform). However, understanding these traits
//! is essential for writing generic code that works with any grid type.
//!
//! ## Example: Generic Function Over Any Partition
//!
//! ```rust
//! use grid1d::{IntervalPartition, HasCoords1D, HasDomain1D, scalars::IntervalId};
//! use grid1d::intervals::IntervalTrait;
//!
//! fn analyze_partition<P>(partition: &P) -> (usize, usize)
//! where
//!     P: IntervalPartition,
//! {
//!     let num_points = *partition.num_points().as_ref();
//!     let num_intervals = *partition.num_intervals().as_ref();
//!     (num_points, num_intervals)
//! }
//! ```

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

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

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

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

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

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

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

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

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

// ============================================================================
// BuildIntervalInPartition
// ============================================================================

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

    /// The interval type used for middle sub-intervals in a partition.
    ///
    /// A "middle" interval is any sub-interval that is neither the first nor the last
    /// in a multi-interval partition. The concrete type encodes the boundary semantics
    /// appropriate for interior intervals of the given domain type.
    ///
    /// For all built-in domain types this defaults to [`IntervalLowerOpenUpperClosed`],
    /// which ensures that every interior partition point belongs to exactly one interval.
    /// Custom implementations may override this to a different interval type provided the
    /// partition-coverage guarantee is maintained.
    type MiddleIntervalInPartition: IntervalFinitePositiveLengthTrait<RealType = Self::RealType>
        + Into<Interval<Self::RealType>>
        + Into<IntervalFinitePositiveLength<Self::RealType>>;

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

    /// Builds a unique interval representing the entire domain.
    ///
    /// For partitions with a single interval, this returns a clone of self.
    fn build_unique_interval(&self) -> Self {
        self.clone()
    }

    /// Constructs the first sub-interval in a partition.
    ///
    /// This method creates the initial sub-interval of a partition, taking into account
    /// the boundary semantics of the domain type.
    ///
    /// # Parameters
    ///
    /// - `lower_bound`: The lower bound of the sub-interval (typically `p_0`)
    /// - `upper_bound`: The upper bound of the sub-interval (typically `p_1`)
    ///
    /// # Returns
    ///
    /// The first sub-interval with appropriate boundary semantics.
    fn build_first_interval(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Self::FirstIntervalInPartition;

    /// Constructs a middle sub-interval in a partition.
    ///
    /// The concrete return type is [`MiddleIntervalInPartition`](Self::MiddleIntervalInPartition),
    /// which encodes the appropriate boundary semantics for interior intervals of this domain.
    /// For all built-in domain types this produces a `(lower, upper]` interval, ensuring
    /// that every interior partition point belongs to exactly one sub-interval.
    ///
    /// # Parameters
    ///
    /// - `lower_bound`: The lower bound of the sub-interval
    /// - `upper_bound`: The upper bound of the sub-interval
    fn build_middle_interval(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Self::MiddleIntervalInPartition;

    /// Constructs the last sub-interval in a partition.
    ///
    /// This method creates the final sub-interval of a partition, taking into account
    /// the boundary semantics of the domain type.
    ///
    /// # Parameters
    ///
    /// - `lower_bound`: The lower bound of the sub-interval (typically `p_{n-1}`)
    /// - `upper_bound`: The upper bound of the sub-interval (typically `p_n`)
    ///
    /// # Returns
    ///
    /// The last sub-interval with appropriate boundary semantics.
    fn build_last_interval(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Self::LastIntervalInPartition;

    /// Whether interior partition boundary points `pₖ` (for `0 < k < n`) belong to the
    /// **right** sub-interval (the one starting at `pₖ`).
    ///
    /// - `true`: sub-intervals are left-closed `[pₖ, pₖ₊₁)` — used by `[a,b]` and `[a,b)` domains.
    /// - `false`: sub-intervals are right-closed `(pₖ, pₖ₊₁]` — used by `(a,b)` and `(a,b]` domains.
    ///
    /// This constant drives the binary-search strategy in
    /// [`IntervalPartition::try_find_interval_id_of_point`].
    const LOWER_BOUND_OWNS_BOUNDARY: bool;
}

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

    #[inline(always)]
    fn build_first_interval(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Self::FirstIntervalInPartition {
        Self::FirstIntervalInPartition::new(lower_bound, upper_bound)
    }

    #[inline(always)]
    fn build_middle_interval(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Self::MiddleIntervalInPartition {
        Self::MiddleIntervalInPartition::new(lower_bound, upper_bound)
    }

    #[inline(always)]
    fn build_last_interval(
        lower_bound: Self::RealType,
        upper_bound: Self::RealType,
    ) -> Self::LastIntervalInPartition {
        Self::LastIntervalInPartition::new(lower_bound, upper_bound)
    }
}

// ============================================================================
// IntervalPartition
// ============================================================================

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

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

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

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

    /// Returns the number of intervals defining the partition.
    ///
    /// This value is equal to `self.num_points() - 1`.
    fn num_intervals(&self) -> NumIntervals {
        NumIntervals::try_new(self.num_points().as_ref() - 1).unwrap()
    }

    /// Returns the `i`-th interval as [`SubIntervalInPartition`].
    ///
    /// # Panics
    ///
    /// In debug mode, panics if the interval ID is out of bounds.
    fn interval(&self, i: &IntervalId) -> SubIntervalInPartition<Self::Domain1D> {
        more_asserts::debug_assert_lt!(i.as_ref(), self.num_intervals().as_ref());

        let i = *i.as_ref();
        let num_intervals = *self.num_intervals().as_ref();

        // Case with only one interval: it's the domain itself
        if num_intervals == 1 {
            let sub_interval = Self::Domain1D::build_unique_interval(self.domain());
            return SubIntervalInPartition::Only(sub_interval);
        }

        let coords = self.coords().as_ref();

        let lower_bound = coords[i].clone();
        let upper_bound = coords[i + 1].clone();

        if i == 0 {
            // First interval
            let sub_interval = Self::Domain1D::build_first_interval(lower_bound, upper_bound);
            SubIntervalInPartition::First(sub_interval)
        } else if i == num_intervals - 1 {
            // Last interval
            let sub_interval = Self::Domain1D::build_last_interval(lower_bound, upper_bound);
            SubIntervalInPartition::Last(sub_interval)
        } else {
            // Middle interval
            SubIntervalInPartition::Middle(Self::Domain1D::build_middle_interval(
                lower_bound,
                upper_bound,
            ))
        }
    }

    /// Returns an iterator over all intervals in the partition.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     IntervalPartition, Grid1D,
    ///     intervals::*,
    ///     scalars::NumIntervals,
    /// };
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 2.0),
    ///     NumIntervals::try_new(2).unwrap()
    /// );
    ///
    /// for (id, interval) in grid.iter_intervals() {
    ///     println!("Interval {}: {:?}", id.as_ref(), interval);
    /// }
    /// ```
    fn iter_intervals(
        &self,
    ) -> impl Iterator<Item = (IntervalId, SubIntervalInPartition<Self::Domain1D>)> + '_ {
        (0..*self.num_intervals().as_ref()).map(move |i| {
            let id = IntervalId::new(i);
            let interval = self.interval(&id);
            (id, interval)
        })
    }

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

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

    /// Returns the ID of the interval that contains the point `x`.
    ///
    /// The containment rule for boundary points depends on the domain type:
    /// - For left-closed sub-intervals (`[a,b]` and `[a,b)` domains), interior boundary point
    ///   `pₖ` belongs to the **right** interval (the one that starts at `pₖ`).
    /// - For right-closed sub-intervals (`(a,b)` and `(a,b]` domains), `pₖ` belongs to the
    ///   **left** interval (the one that ends at `pₖ`).
    ///
    /// # Panics
    ///
    /// In debug builds, panics if `x` is not contained within the domain.
    /// Use [`try_find_interval_id_of_point`](IntervalPartition::try_find_interval_id_of_point)
    /// for a safe alternative.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Grid1DUniform, IntervalPartition,
    ///     intervals::*,
    ///     scalars::{IntervalId, NumIntervals},
    /// };
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1DUniform::new(
    ///     IntervalClosed::new(0.0, 2.0),
    ///     NumIntervals::try_new(2).unwrap()
    /// );
    ///
    /// assert_eq!(grid.find_interval_id_of_point(&0.5), IntervalId::new(0));
    /// assert_eq!(grid.find_interval_id_of_point(&1.0), IntervalId::new(1)); // Boundary → right (domain is [a,b])
    /// assert_eq!(grid.find_interval_id_of_point(&1.5), IntervalId::new(1));
    /// ```
    fn find_interval_id_of_point(&self, x: &Self::Point1DType) -> IntervalId {
        debug_assert!(
            self.domain().contains_point(x),
            "The value {} is not in the interval {:?} spanned by the grid1d!",
            x,
            self.domain()
        );

        self.try_find_interval_id_of_point(x)
            .expect("Precondition violated: point is outside the domain.")
    }

    /// Finds the ID of the interval containing `x`, returning `None` if outside the domain.
    ///
    /// This is a safe alternative to [`find_interval_id_of_point`](IntervalPartition::find_interval_id_of_point).
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{
    ///     Grid1DUniform, IntervalPartition,
    ///     intervals::*,
    ///     scalars::{IntervalId, NumIntervals},
    /// };
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1DUniform::new(
    ///     IntervalClosed::new(0.0, 2.0),
    ///     NumIntervals::try_new(2).unwrap()
    /// );
    ///
    /// assert_eq!(grid.try_find_interval_id_of_point(&1.0), Some(IntervalId::new(1))); // Boundary → right for [a,b]
    /// assert_eq!(grid.try_find_interval_id_of_point(&3.0), None); // Outside domain
    /// ```
    fn try_find_interval_id_of_point(&self, x: &Self::Point1DType) -> Option<IntervalId>;

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

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

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

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

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

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

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

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

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

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

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

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

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

        intersecting_intervals
    }

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

    /// Performs uniform refinement where all intervals are subdivided equally.
    ///
    /// # Parameters
    ///
    /// - `num_extra_points_each_interval`: Number of additional points to insert in each interval
    ///
    /// # Returns
    ///
    /// A [`crate::Grid1DUniformRefinement`] containing the refined grid
    /// and bidirectional mappings.
    ///
    /// # Example
    ///
    /// ```rust
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 1.0),
    ///     NumIntervals::try_new(4).unwrap()
    /// );
    ///
    /// // Add 1 point per interval (doubles resolution)
    /// let refined = grid.refine_uniform(&PositiveNumPoints1D::try_new(1).unwrap());
    /// assert_eq!(refined.refined_grid().num_intervals().as_ref(), &8);
    /// ```
    fn refine_uniform(
        self,
        num_extra_points_each_interval: &PositiveNumPoints1D,
    ) -> Grid1DUniformRefinement<Self>;

    /// Performs selective refinement of specified intervals.
    ///
    /// Only the intervals in the map are refined; others remain unchanged.
    /// The [`BTreeMap`] ensures unique interval IDs and ordered iteration.
    ///
    /// # Parameters
    ///
    /// - `intervals_to_refine`: Map from interval IDs to number of extra points
    ///
    /// # Panics
    ///
    /// In debug mode, panics if any interval ID is out of bounds.
    ///
    /// # Example
    ///
    /// ```rust
    /// use std::collections::BTreeMap;
    /// use grid1d::{*, intervals::*, scalars::*};
    /// use try_create::TryNew;
    ///
    /// let grid = Grid1D::uniform(
    ///     IntervalClosed::new(0.0, 1.0),
    ///     NumIntervals::try_new(4).unwrap()
    /// );
    ///
    /// let mut plan = BTreeMap::new();
    /// plan.insert(IntervalId::new(1), PositiveNumPoints1D::try_new(2).unwrap());
    /// plan.insert(IntervalId::new(3), PositiveNumPoints1D::try_new(1).unwrap());
    ///
    /// let refined = grid.refine(&plan);
    /// assert!(refined.was_refined(&IntervalId::new(1)));
    /// assert!(!refined.was_refined(&IntervalId::new(0)));
    /// ```
    fn refine(
        self,
        intervals_to_refine: &BTreeMap<IntervalId, PositiveNumPoints1D>,
    ) -> Grid1DNonUniformRefinement<Self>;
}

// ============================================================================
// TransformedPoints1D
// ============================================================================

/// A trait for point transformations between parametric and physical domains.
///
/// This trait provides access to coordinates in both a reference (parametric) interval
/// and a physical (target) interval. It is useful for finite element methods, isoparametric
/// mappings, and coordinate transformations.
///
/// # Type Parameters
///
/// - `ParamStorage`: Storage type for parametric coordinates (owned, borrowed, Arc, Rc)
/// - `PhysStorage`: Storage type for physical coordinates
///
/// # Associated Types
///
/// - `ParamDomainType`: The interval type for the parametric domain (e.g., `[-1, 1]`)
/// - `PhysDomainType`: The interval type for the physical domain
///
/// # Example
///
/// ```rust
/// use grid1d::{Coords1D, Coords1DInDomain, intervals::IntervalClosed};
/// use grid1d::TransformedPoints1D;
/// use sorted_vec::partial::SortedSet;
///
/// // A struct implementing TransformedPoints1D could map
/// // reference coordinates [-1, 1] to physical coordinates [a, b]
/// ```
pub trait TransformedPoints1D<ParamStorage, PhysStorage>
where
    ParamStorage:
        std::borrow::Borrow<Coords1D<<Self::ParamDomainType as IntervalBoundsRuntime>::RealType>>,
    PhysStorage:
        std::borrow::Borrow<Coords1D<<Self::PhysDomainType as IntervalBoundsRuntime>::RealType>>,
{
    /// The interval type for the parametric domain.
    ///
    /// This represents the reference interval from which points are transformed,
    /// typically a standard interval like `[-1, 1]` or `[0, 1]`.
    type ParamDomainType: IntervalTrait;

    /// The interval type for the physical domain.
    ///
    /// This represents the target physical interval to which parametric points are mapped.
    /// Must use the same real number type as the parametric domain.
    type PhysDomainType: IntervalTrait<
        RealType = <Self::ParamDomainType as IntervalBoundsRuntime>::RealType,
    >;

    /// Returns the coordinates in the parametric interval.
    ///
    /// These values are sorted in ascending order.
    fn points_param_interval(&self) -> &Coords1DInDomain<Self::ParamDomainType, ParamStorage>;

    /// Returns the coordinates in the physical interval.
    ///
    /// These values are sorted in ascending order.
    fn points_phys_interval(&self) -> &Coords1DInDomain<Self::PhysDomainType, PhysStorage>;
}

// ============================================================================
// Type Aliases for TransformedPoints1D
// ============================================================================

/// Type alias for [`TransformedPoints1D`] with owned coordinate storage for both domains.
///
/// This is the most common usage pattern where both parametric and physical coordinate
/// storage is owned by the implementing type.
pub type TransformedPoints1DOwned<ParamDomain, PhysDomain> = dyn TransformedPoints1D<
        Coords1D<<ParamDomain as IntervalBoundsRuntime>::RealType>,
        Coords1D<<PhysDomain as IntervalBoundsRuntime>::RealType>,
        ParamDomainType = ParamDomain,
        PhysDomainType = PhysDomain,
    >;

/// Type alias for [`TransformedPoints1D`] with borrowed coordinate storage for both domains.
///
/// Useful for zero-cost access patterns where coordinates are borrowed from external sources.
pub type TransformedPoints1DBorrowed<'a, ParamDomain, PhysDomain> = dyn TransformedPoints1D<
        &'a Coords1D<<ParamDomain as IntervalBoundsRuntime>::RealType>,
        &'a Coords1D<<PhysDomain as IntervalBoundsRuntime>::RealType>,
        ParamDomainType = ParamDomain,
        PhysDomainType = PhysDomain,
    > + 'a;

/// Type alias for [`TransformedPoints1D`] with Arc coordinate storage for both domains.
///
/// Enables thread-safe sharing of coordinate data across multiple threads.
pub type TransformedPoints1DArc<ParamDomain, PhysDomain> = dyn TransformedPoints1D<
        Arc<Coords1D<<ParamDomain as IntervalBoundsRuntime>::RealType>>,
        Arc<Coords1D<<PhysDomain as IntervalBoundsRuntime>::RealType>>,
        ParamDomainType = ParamDomain,
        PhysDomainType = PhysDomain,
    >;

/// Type alias for [`TransformedPoints1D`] with Rc coordinate storage for both domains.
///
/// Enables efficient single-threaded sharing of coordinate data with reference counting.
pub type TransformedPoints1DRc<ParamDomain, PhysDomain> = dyn TransformedPoints1D<
        std::rc::Rc<Coords1D<<ParamDomain as IntervalBoundsRuntime>::RealType>>,
        std::rc::Rc<Coords1D<<PhysDomain as IntervalBoundsRuntime>::RealType>>,
        ParamDomainType = ParamDomain,
        PhysDomainType = PhysDomain,
    >;

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

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

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

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

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

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

    #[test]
    fn test_build_interval_in_partition_closed() {
        let first = IntervalClosed::build_first_interval(0.0, 1.0);
        let middle = IntervalClosed::build_middle_interval(1.0, 2.0);
        let last = IntervalClosed::build_last_interval(2.0, 3.0);

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

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

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

    #[test]
    fn test_build_interval_in_partition_open() {
        let first = IntervalOpen::build_first_interval(0.0, 1.0);
        let last = IntervalOpen::build_last_interval(2.0, 3.0);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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