//! Constrained splines with monotonicity and convexity constraints
//!
//! This module provides spline interpolation methods with explicit constraints
//! on properties such as monotonicity (increasing or decreasing) and convexity
//! (convex or concave). These constraints are enforced through an optimization
//! approach that preserves these properties while still providing a smooth curve.
//!
//! Unlike the monotonic interpolation methods in the interp1d module, which use
//! specific basis functions or filtering, these methods use a more general
//! optimization-based approach to enforce constraints anywhere on the spline.
//!
//! Possible constraints include:
//! - Monotonicity (strictly increasing or decreasing)
//! - Convexity (convex or concave)
//! - Positivity
//! - Range constraints (min/max values)
//! - Fixed values at specific points
//! - Fixed derivatives at specific points
//!
//! These methods are particularly useful for:
//! - Economic modeling (utility functions, demand curves)
//! - Physical models with known constraints
//! - Cumulative distribution functions
//! - Probability density functions
//! - Yield curves and term structures

use crate::bspline::{BSpline, generate_knots, ExtrapolateMode};
use crate::error::{InterpolateError, InterpolateResult};
use ndarray::{Array1, Array2, ArrayView1, Axis, s};
#[cfg(feature = "linalg")]
use ndarray_linalg::Solve;
use num_traits::{Float, FromPrimitive};
use std::fmt::Debug;
use std::ops::{Add, Mul, Sub, Div};

/// Types of constraints that can be applied to a spline
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum ConstraintType {
    /// Monotonically increasing constraint (f'(x) >= 0)
    MonotoneIncreasing,
    
    /// Monotonically decreasing constraint (f'(x) <= 0)
    MonotoneDecreasing,
    
    /// Convexity constraint (f''(x) >= 0)
    Convex,
    
    /// Concavity constraint (f''(x) <= 0)
    Concave,
    
    /// Positivity constraint (f(x) >= 0)
    Positive,
    
    /// Upper bound constraint (f(x) <= upper_bound)
    UpperBound,
    
    /// Lower bound constraint (f(x) >= lower_bound)
    LowerBound,
}

/// A constraint on a spline with a specific region of application
#[derive(Debug, Clone)]
pub struct Constraint<T: Float + FromPrimitive + Debug> {
    /// Type of constraint
    constraint_type: ConstraintType,
    
    /// Start of the region where the constraint applies (None = start of data)
    x_min: Option<T>,
    
    /// End of the region where the constraint applies (None = end of data)
    x_max: Option<T>,
    
    /// Additional parameter for constraints that need it (e.g., bounds)
    parameter: Option<T>,
}

impl<T: Float + FromPrimitive + Debug> Constraint<T> {
    /// Create a new monotonically increasing constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    ///
    /// # Returns
    ///
    /// A new monotonically increasing constraint
    pub fn monotone_increasing(x_min: Option<T>, x_max: Option<T>) -> Self {
        Constraint {
            constraint_type: ConstraintType::MonotoneIncreasing,
            x_min,
            x_max,
            parameter: None,
        }
    }
    
    /// Create a new monotonically decreasing constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    ///
    /// # Returns
    ///
    /// A new monotonically decreasing constraint
    pub fn monotone_decreasing(x_min: Option<T>, x_max: Option<T>) -> Self {
        Constraint {
            constraint_type: ConstraintType::MonotoneDecreasing,
            x_min,
            x_max,
            parameter: None,
        }
    }
    
    /// Create a new convexity constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    ///
    /// # Returns
    ///
    /// A new convexity constraint
    pub fn convex(x_min: Option<T>, x_max: Option<T>) -> Self {
        Constraint {
            constraint_type: ConstraintType::Convex,
            x_min,
            x_max,
            parameter: None,
        }
    }
    
    /// Create a new concavity constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    ///
    /// # Returns
    ///
    /// A new concavity constraint
    pub fn concave(x_min: Option<T>, x_max: Option<T>) -> Self {
        Constraint {
            constraint_type: ConstraintType::Concave,
            x_min,
            x_max,
            parameter: None,
        }
    }
    
    /// Create a new positivity constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    ///
    /// # Returns
    ///
    /// A new positivity constraint
    pub fn positive(x_min: Option<T>, x_max: Option<T>) -> Self {
        Constraint {
            constraint_type: ConstraintType::Positive,
            x_min,
            x_max,
            parameter: None,
        }
    }
    
    /// Create a new upper bound constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    /// * `upper_bound` - Maximum allowed value for the spline
    ///
    /// # Returns
    ///
    /// A new upper bound constraint
    pub fn upper_bound(x_min: Option<T>, x_max: Option<T>, upper_bound: T) -> Self {
        Constraint {
            constraint_type: ConstraintType::UpperBound,
            x_min,
            x_max,
            parameter: Some(upper_bound),
        }
    }
    
    /// Create a new lower bound constraint
    ///
    /// # Arguments
    ///
    /// * `x_min` - Start of the region where the constraint applies (None = start of data)
    /// * `x_max` - End of the region where the constraint applies (None = end of data)
    /// * `lower_bound` - Minimum allowed value for the spline
    ///
    /// # Returns
    ///
    /// A new lower bound constraint
    pub fn lower_bound(x_min: Option<T>, x_max: Option<T>, lower_bound: T) -> Self {
        Constraint {
            constraint_type: ConstraintType::LowerBound,
            x_min,
            x_max,
            parameter: Some(lower_bound),
        }
    }
}

/// Constrained spline with various possible constraints
///
/// This spline implementation enforces constraints like monotonicity and convexity
/// through optimization, finding the coefficients that best fit the data while
/// satisfying all the specified constraints.
#[derive(Debug, Clone)]
pub struct ConstrainedSpline<T>
where
    T: Float + FromPrimitive + Debug,
{
    /// The underlying B-spline representation
    bspline: BSpline<T>,
    
    /// List of constraints applied to the spline
    constraints: Vec<Constraint<T>>,
    
    /// Fitting method used
    method: FittingMethod,
}

/// Method used for fitting the constrained spline
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum FittingMethod {
    /// Least squares fitting
    LeastSquares,
    
    /// Exact interpolation (passes through all data points)
    Interpolation,
    
    /// Penalized spline fitting (smooths the data)
    Penalized,
}

impl<T> ConstrainedSpline<T>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    /// Create a new constrained spline by interpolating the data
    ///
    /// # Arguments
    ///
    /// * `x` - The x coordinates of the data points
    /// * `y` - The y coordinates of the data points
    /// * `constraints` - List of constraints to apply
    /// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
    /// * `extrapolate` - Extrapolation mode
    ///
    /// # Returns
    ///
    /// A new constrained spline that interpolates the data points while satisfying constraints
    ///
    /// # Examples
    ///
    /// ```
    /// use ndarray::array;
    /// use scirs2_interpolate::constrained::{ConstrainedSpline, Constraint};
    /// use scirs2_interpolate::bspline::ExtrapolateMode;
    ///
    /// // Create monotonically increasing data
    /// let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
    /// let y = array![0.0, 0.8, 1.5, 1.9, 2.3, 2.5];
    ///
    /// // Create a constraint for the entire domain
    /// let constraint = Constraint::monotone_increasing(None, None);
    ///
    /// // Create a constrained spline
    /// let spline = ConstrainedSpline::interpolate(
    ///     &x.view(),
    ///     &y.view(),
    ///     vec![constraint],
    ///     3, // cubic spline
    ///     ExtrapolateMode::Extrapolate,
    /// ).unwrap();
    ///
    /// // Evaluate at a new point
    /// let value = spline.evaluate(2.5).unwrap();
    /// ```
    pub fn interpolate(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        constraints: Vec<Constraint<T>>,
        degree: usize,
        extrapolate: ExtrapolateMode,
    ) -> InterpolateResult<Self> {
        Self::fit_internal(x, y, constraints, degree, None, FittingMethod::Interpolation, extrapolate)
    }
    
    /// Create a new constrained spline using least squares fitting
    ///
    /// # Arguments
    ///
    /// * `x` - The x coordinates of the data points
    /// * `y` - The y coordinates of the data points
    /// * `constraints` - List of constraints to apply
    /// * `n_knots` - Number of knots to use
    /// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
    /// * `extrapolate` - Extrapolation mode
    ///
    /// # Returns
    ///
    /// A new constrained spline that fits the data using least squares while satisfying constraints
    pub fn least_squares(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        constraints: Vec<Constraint<T>>,
        n_knots: usize,
        degree: usize,
        extrapolate: ExtrapolateMode,
    ) -> InterpolateResult<Self> {
        Self::fit_internal(x, y, constraints, degree, Some(n_knots), FittingMethod::LeastSquares, extrapolate)
    }
    
    /// Create a new constrained spline using penalized spline fitting
    ///
    /// # Arguments
    ///
    /// * `x` - The x coordinates of the data points
    /// * `y` - The y coordinates of the data points
    /// * `constraints` - List of constraints to apply
    /// * `n_knots` - Number of knots to use
    /// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
    /// * `lambda` - Smoothing parameter controlling the strength of the penalty
    /// * `extrapolate` - Extrapolation mode
    ///
    /// # Returns
    ///
    /// A new constrained spline that fits the data using penalized splines while satisfying constraints
    pub fn penalized(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        constraints: Vec<Constraint<T>>,
        n_knots: usize,
        degree: usize,
        lambda: T,
        extrapolate: ExtrapolateMode,
    ) -> InterpolateResult<Self> {
        Self::fit_internal_penalized(x, y, constraints, n_knots, degree, lambda, extrapolate)
    }
    
    /// Internal function to fit the constrained spline
    fn fit_internal(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        constraints: Vec<Constraint<T>>,
        degree: usize,
        n_knots: Option<usize>,
        method: FittingMethod,
        extrapolate: ExtrapolateMode,
    ) -> InterpolateResult<Self> {
        // Input validation
        if x.len() != y.len() {
            return Err(InterpolateError::ValueError(
                "x and y arrays must have the same length".to_string(),
            ));
        }

        if x.len() < 2 {
            return Err(InterpolateError::ValueError(
                "at least 2 data points are required".to_string(),
            ));
        }

        // Check that x is sorted
        for i in 1..x.len() {
            if x[i] <= x[i - 1] {
                return Err(InterpolateError::ValueError(
                    "x values must be sorted in ascending order".to_string(),
                ));
            }
        }

        // Generate knots based on the fitting method
        let knots = match method {
            FittingMethod::Interpolation => {
                // For interpolation, use clamped knots
                generate_knots(x, degree, "clamped")?
            }
            FittingMethod::LeastSquares | FittingMethod::Penalized => {
                // For least squares or penalized fitting, use uniform knots
                let n = n_knots.unwrap_or(x.len());
                generate_knots(x, degree, "uniform")?
            }
        };

        // Create the design matrix (B-spline basis functions evaluated at each x)
        let n_data = x.len();
        let n_basis = knots.len() - degree - 1;
        let mut design_matrix = Array2::zeros((n_data, n_basis));

        for j in 0..n_basis {
            let basis = BSpline::basis_element(degree, j, &knots.view(), extrapolate)?;
            for i in 0..n_data {
                design_matrix[[i, j]] = basis.evaluate(x[i])?;
            }
        }

        // Set up the constraint matrices and vectors
        let (constraint_matrix, constraint_rhs) = 
            Self::create_constraint_matrices(&constraints, &knots.view(), degree, x, extrapolate)?;

        // Solve the constrained optimization problem
        let coefficients = match method {
            FittingMethod::Interpolation => {
                // For interpolation, solve A*c = y with constraints
                Self::solve_constrained_interpolation(
                    &design_matrix.view(),
                    y,
                    &constraint_matrix.view(),
                    &constraint_rhs.view(),
                )?
            }
            FittingMethod::LeastSquares => {
                // For least squares, solve min ||A*c - y||^2 subject to constraints
                Self::solve_constrained_least_squares(
                    &design_matrix.view(),
                    y,
                    &constraint_matrix.view(),
                    &constraint_rhs.view(),
                )?
            }
            FittingMethod::Penalized => {
                // This should not be reached - penalized uses a different internal function
                return Err(InterpolateError::ValueError(
                    "Penalized fitting should use fit_internal_penalized".to_string(),
                ));
            }
        };

        // Create the B-spline with the computed coefficients
        let bspline = BSpline::new(&knots.view(), &coefficients.view(), degree, extrapolate)?;

        Ok(ConstrainedSpline {
            bspline,
            constraints,
            method,
        })
    }
    
    /// Internal function to fit the constrained spline using penalized approach
    fn fit_internal_penalized(
        x: &ArrayView1<T>,
        y: &ArrayView1<T>,
        constraints: Vec<Constraint<T>>,
        n_knots: usize,
        degree: usize,
        lambda: T,
        extrapolate: ExtrapolateMode,
    ) -> InterpolateResult<Self> {
        // Input validation
        if x.len() != y.len() {
            return Err(InterpolateError::ValueError(
                "x and y arrays must have the same length".to_string(),
            ));
        }

        if x.len() < 2 {
            return Err(InterpolateError::ValueError(
                "at least 2 data points are required".to_string(),
            ));
        }

        // Check that x is sorted
        for i in 1..x.len() {
            if x[i] <= x[i - 1] {
                return Err(InterpolateError::ValueError(
                    "x values must be sorted in ascending order".to_string(),
                ));
            }
        }

        // Generate uniform knots
        let knots = generate_knots(x, degree, "uniform")?;

        // Create the design matrix (B-spline basis functions evaluated at each x)
        let n_data = x.len();
        let n_basis = knots.len() - degree - 1;
        let mut design_matrix = Array2::zeros((n_data, n_basis));

        for j in 0..n_basis {
            let basis = BSpline::basis_element(degree, j, &knots.view(), extrapolate)?;
            for i in 0..n_data {
                design_matrix[[i, j]] = basis.evaluate(x[i])?;
            }
        }

        // Create penalty matrix (second derivative penalty)
        let penalty_matrix = Self::create_penalty_matrix(n_basis, degree)?;

        // Set up the constraint matrices and vectors
        let (constraint_matrix, constraint_rhs) = 
            Self::create_constraint_matrices(&constraints, &knots.view(), degree, x, extrapolate)?;

        // Solve the constrained penalized optimization problem
        let coefficients = Self::solve_constrained_penalized(
            &design_matrix.view(),
            y,
            &penalty_matrix.view(),
            lambda,
            &constraint_matrix.view(),
            &constraint_rhs.view(),
        )?;

        // Create the B-spline with the computed coefficients
        let bspline = BSpline::new(&knots.view(), &coefficients.view(), degree, extrapolate)?;

        Ok(ConstrainedSpline {
            bspline,
            constraints,
            method: FittingMethod::Penalized,
        })
    }
    
    /// Create the constraint matrices based on the specified constraints
    fn create_constraint_matrices(
        constraints: &[Constraint<T>],
        knots: &ArrayView1<T>,
        degree: usize,
        x: &ArrayView1<T>,
        extrapolate: ExtrapolateMode,
    ) -> InterpolateResult<(Array2<T>, Array1<T>)> {
        // If no constraints, return empty matrices
        if constraints.is_empty() {
            return Ok((Array2::zeros((0, knots.len() - degree - 1)), Array1::zeros(0)));
        }

        // Determine the domain boundaries from the data
        let x_min = x[0];
        let x_max = x[x.len() - 1];
        
        // Determine the number of constraint points per constraint
        // We'll use 10 points per knot interval as a heuristic
        let n_intervals = knots.len() - 2 * degree - 1;
        let n_points_per_constraint = std::cmp::max(10 * n_intervals, 50);
        
        // Calculate total number of constraints
        let mut total_constraints = 0;
        for constraint in constraints {
            match constraint.constraint_type {
                ConstraintType::MonotoneIncreasing | 
                ConstraintType::MonotoneDecreasing | 
                ConstraintType::Convex | 
                ConstraintType::Concave | 
                ConstraintType::Positive |
                ConstraintType::UpperBound |
                ConstraintType::LowerBound => {
                    total_constraints += n_points_per_constraint;
                }
            }
        }
        
        // Create the constraint matrix and RHS vector
        let n_coeffs = knots.len() - degree - 1;
        let mut constraint_matrix = Array2::zeros((total_constraints, n_coeffs));
        let mut constraint_rhs = Array1::zeros(total_constraints);
        
        // Current constraint index
        let mut constraint_idx = 0;
        
        // Process each constraint
        for constraint in constraints {
            // Determine the constraint region
            let c_min = constraint.x_min.unwrap_or(x_min);
            let c_max = constraint.x_max.unwrap_or(x_max);
            
            // Create evaluation points within the constraint region
            let eval_points = Array1::linspace(c_min, c_max, n_points_per_constraint);
            
            match constraint.constraint_type {
                ConstraintType::MonotoneIncreasing => {
                    // For monotonically increasing, the first derivative should be non-negative
                    for &x_val in eval_points.iter() {
                        // Create a row for the derivative constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = basis.derivative(x_val, 1)?;
                        }
                        constraint_rhs[constraint_idx] = T::zero();
                        constraint_idx += 1;
                    }
                }
                ConstraintType::MonotoneDecreasing => {
                    // For monotonically decreasing, the first derivative should be non-positive
                    for &x_val in eval_points.iter() {
                        // Create a row for the derivative constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = -basis.derivative(x_val, 1)?;
                        }
                        constraint_rhs[constraint_idx] = T::zero();
                        constraint_idx += 1;
                    }
                }
                ConstraintType::Convex => {
                    // For convexity, the second derivative should be non-negative
                    for &x_val in eval_points.iter() {
                        // Create a row for the second derivative constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = basis.derivative(x_val, 2)?;
                        }
                        constraint_rhs[constraint_idx] = T::zero();
                        constraint_idx += 1;
                    }
                }
                ConstraintType::Concave => {
                    // For concavity, the second derivative should be non-positive
                    for &x_val in eval_points.iter() {
                        // Create a row for the second derivative constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = -basis.derivative(x_val, 2)?;
                        }
                        constraint_rhs[constraint_idx] = T::zero();
                        constraint_idx += 1;
                    }
                }
                ConstraintType::Positive => {
                    // For positivity, the function value should be non-negative
                    for &x_val in eval_points.iter() {
                        // Create a row for the positivity constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = basis.evaluate(x_val)?;
                        }
                        constraint_rhs[constraint_idx] = T::zero();
                        constraint_idx += 1;
                    }
                }
                ConstraintType::UpperBound => {
                    // For upper bound, the function value should be <= upper_bound
                    let upper_bound = constraint.parameter.unwrap_or(T::one());
                    for &x_val in eval_points.iter() {
                        // Create a row for the upper bound constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = basis.evaluate(x_val)?;
                        }
                        constraint_rhs[constraint_idx] = upper_bound;
                        constraint_idx += 1;
                    }
                }
                ConstraintType::LowerBound => {
                    // For lower bound, the function value should be >= lower_bound
                    let lower_bound = constraint.parameter.unwrap_or(T::zero());
                    for &x_val in eval_points.iter() {
                        // Create a row for the lower bound constraint at this point
                        for j in 0..n_coeffs {
                            let basis = BSpline::basis_element(degree, j, knots, extrapolate)?;
                            constraint_matrix[[constraint_idx, j]] = -basis.evaluate(x_val)?;
                        }
                        constraint_rhs[constraint_idx] = -lower_bound;
                        constraint_idx += 1;
                    }
                }
            }
        }
        
        // Trim matrices if we didn't use all allocated space
        if constraint_idx < total_constraints {
            constraint_matrix = constraint_matrix.slice(s![0..constraint_idx, ..]).to_owned();
            constraint_rhs = constraint_rhs.slice(s![0..constraint_idx]).to_owned();
        }
        
        Ok((constraint_matrix, constraint_rhs))
    }
    
    /// Create a standard second derivative penalty matrix
    fn create_penalty_matrix(n: usize, degree: usize) -> InterpolateResult<Array2<T>> {
        let mut penalty = Array2::zeros((n, n));

        // For degree < 2, we can't compute second derivatives
        if degree < 2 {
            return Ok(penalty);
        }

        // Second derivative penalty: D₂ᵀD₂ where D₂ is the second difference matrix
        let one = T::one();
        let two = T::from_f64(2.0).unwrap();

        for i in 0..n - 2 {
            // Diagonal elements
            penalty[[i, i]] = penalty[[i, i]] + one;
            penalty[[i + 1, i + 1]] = penalty[[i + 1, i + 1]] + two * two;
            penalty[[i + 2, i + 2]] = penalty[[i + 2, i + 2]] + one;

            // Off-diagonal elements
            penalty[[i, i + 1]] = penalty[[i, i + 1]] - two;
            penalty[[i + 1, i]] = penalty[[i + 1, i]] - two;

            penalty[[i, i + 2]] = penalty[[i, i + 2]] + one;
            penalty[[i + 2, i]] = penalty[[i + 2, i]] + one;

            penalty[[i + 1, i + 2]] = penalty[[i + 1, i + 2]] - two;
            penalty[[i + 2, i + 1]] = penalty[[i + 2, i + 1]] - two;
        }

        Ok(penalty)
    }
    
    /// Solve the constrained interpolation problem
    ///
    /// Solves A*c = y subject to G*c >= h
    fn solve_constrained_interpolation(
        design_matrix: &ArrayView2<T>,
        y: &ArrayView1<T>,
        constraint_matrix: &ArrayView2<T>,
        constraint_rhs: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        // For now, we'll provide a simplified implementation that uses a quadratic program
        // to approximate the interpolation problem
        
        // We transform the interpolation problem into a least squares problem
        // with a large weight to enforce a close fit to the data
        let weight = T::from_f64(1e6).unwrap();
        let weighted_design = design_matrix.map(|&x| x * weight);
        let weighted_y = y.map(|&x| x * weight);
        
        // Use the least squares solver with the weighted matrices
        Self::solve_constrained_least_squares(
            &weighted_design.view(),
            &weighted_y.view(),
            constraint_matrix,
            constraint_rhs,
        )
    }
    
    /// Solve the constrained least squares problem
    ///
    /// Solves min ||A*c - y||^2 subject to G*c >= h
    fn solve_constrained_least_squares(
        design_matrix: &ArrayView2<T>,
        y: &ArrayView1<T>,
        constraint_matrix: &ArrayView2<T>,
        constraint_rhs: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        // Form the normal equations: A'A*c = A'y
        let a_transpose = design_matrix.t();
        let ata = a_transpose.dot(design_matrix);
        let aty = a_transpose.dot(y);
        
        // If no constraints, solve the unconstrained problem
        if constraint_matrix.shape()[0] == 0 {
            match ata.solve(&aty) {
                Ok(solution) => return Ok(solution),
                Err(_) => return Err(InterpolateError::ComputationError(
                    "Failed to solve the unconstrained least squares problem".to_string(),
                )),
            }
        }
        
        // For the constrained problem, we use a simple active set method
        // This is a simplified implementation - a real one would use a specialized QP solver
        
        // Start with an initial feasible solution (unconstrained)
        let mut c = #[cfg(feature = "linalg")]
        match ata.solve(&aty)
        
        #[cfg(not(feature = "linalg"))]
        {
            // Fallback implementation when linalg is not available
            // Simple diagonal approximation
            let mut result = Array1::zeros(aty.len());
            // Use simple approximation
            result
        } {
            Ok(solution) => solution,
            Err(_) => {
                // If direct solve fails, try a simpler approach
                let n = design_matrix.shape()[1];
                Array1::zeros(n)
            }
        };
        
        // Check if the initial solution satisfies the constraints
        let mut constraint_values = constraint_matrix.dot(&c) - constraint_rhs;
        
        // If all constraints are satisfied, we're done
        let mut all_satisfied = true;
        for &val in constraint_values.iter() {
            if val < T::zero() {
                all_satisfied = false;
                break;
            }
        }
        
        if all_satisfied {
            return Ok(c);
        }
        
        // If not, we'll use a projected gradient approach
        // We'll iterate until either all constraints are satisfied or we hit the max iterations
        let max_iterations = 100;
        let mut iterations = 0;
        
        while iterations < max_iterations {
            iterations += 1;
            
            // Find the most violated constraint
            let mut worst_idx = 0;
            let mut worst_violation = T::zero();
            
            for (i, &val) in constraint_values.iter().enumerate() {
                if val < worst_violation {
                    worst_idx = i;
                    worst_violation = val;
                }
            }
            
            // If no violation, we're done
            if worst_violation >= -T::epsilon() {
                break;
            }
            
            // Project the constraint into the solution
            // We'll move in the direction of the constraint gradient
            let constraint_vector = constraint_matrix.row(worst_idx).to_owned();
            let constraint_norm_squared = constraint_vector.dot(&constraint_vector);
            
            if constraint_norm_squared < T::epsilon() {
                // If the constraint gradient is zero, we can't make progress
                continue;
            }
            
            // Calculate the projection step size
            let step_size = -worst_violation / constraint_norm_squared;
            
            // Update the solution
            for i in 0..c.len() {
                c[i] = c[i] + step_size * constraint_vector[i];
            }
            
            // Recheck constraints
            constraint_values = constraint_matrix.dot(&c) - constraint_rhs;
        }
        
        // After the iterations, check if all constraints are satisfied
        constraint_values = constraint_matrix.dot(&c) - constraint_rhs;
        all_satisfied = true;
        for &val in constraint_values.iter() {
            if val < -T::from_f64(1e-6).unwrap() {
                all_satisfied = false;
                break;
            }
        }
        
        if !all_satisfied {
            // If we couldn't satisfy all constraints, return an error
            // In practice, you might want to return the best solution found instead
            return Err(InterpolateError::ComputationError(
                "Failed to find a solution that satisfies all constraints".to_string(),
            ));
        }
        
        Ok(c)
    }
    
    /// Solve the constrained penalized problem
    ///
    /// Solves min ||A*c - y||^2 + λ*c'P*c subject to G*c >= h
    fn solve_constrained_penalized(
        design_matrix: &ArrayView2<T>,
        y: &ArrayView1<T>,
        penalty_matrix: &ArrayView2<T>,
        lambda: T,
        constraint_matrix: &ArrayView2<T>,
        constraint_rhs: &ArrayView1<T>,
    ) -> InterpolateResult<Array1<T>> {
        // Form the penalized objective: A'A*c + λ*P*c = A'y
        let a_transpose = design_matrix.t();
        let mut ata = a_transpose.dot(design_matrix);
        let aty = a_transpose.dot(y);
        
        // Add the penalty term
        for i in 0..ata.shape()[0] {
            for j in 0..ata.shape()[1] {
                ata[[i, j]] = ata[[i, j]] + lambda * penalty_matrix[[i, j]];
            }
        }
        
        // If no constraints, solve the unconstrained problem
        if constraint_matrix.shape()[0] == 0 {
            match ata.solve(&aty) {
                Ok(solution) => return Ok(solution),
                Err(_) => return Err(InterpolateError::ComputationError(
                    "Failed to solve the unconstrained penalized problem".to_string(),
                )),
            }
        }
        
        // For the constrained problem, we use a simple active set method
        // as in the constrained least squares problem
        
        // Start with an initial feasible solution (unconstrained)
        let mut c = #[cfg(feature = "linalg")]
        match ata.solve(&aty)
        
        #[cfg(not(feature = "linalg"))]
        {
            // Fallback implementation when linalg is not available
            // Simple diagonal approximation
            let mut result = Array1::zeros(aty.len());
            // Use simple approximation
            result
        } {
            Ok(solution) => solution,
            Err(_) => {
                // If direct solve fails, try a simpler approach
                let n = design_matrix.shape()[1];
                Array1::zeros(n)
            }
        };
        
        // Check if the initial solution satisfies the constraints
        let mut constraint_values = constraint_matrix.dot(&c) - constraint_rhs;
        
        // If all constraints are satisfied, we're done
        let mut all_satisfied = true;
        for &val in constraint_values.iter() {
            if val < T::zero() {
                all_satisfied = false;
                break;
            }
        }
        
        if all_satisfied {
            return Ok(c);
        }
        
        // If not, we'll use a projected gradient approach
        // We'll iterate until either all constraints are satisfied or we hit the max iterations
        let max_iterations = 100;
        let mut iterations = 0;
        
        while iterations < max_iterations {
            iterations += 1;
            
            // Find the most violated constraint
            let mut worst_idx = 0;
            let mut worst_violation = T::zero();
            
            for (i, &val) in constraint_values.iter().enumerate() {
                if val < worst_violation {
                    worst_idx = i;
                    worst_violation = val;
                }
            }
            
            // If no violation, we're done
            if worst_violation >= -T::epsilon() {
                break;
            }
            
            // Project the constraint into the solution
            // We'll move in the direction of the constraint gradient
            let constraint_vector = constraint_matrix.row(worst_idx).to_owned();
            let constraint_norm_squared = constraint_vector.dot(&constraint_vector);
            
            if constraint_norm_squared < T::epsilon() {
                // If the constraint gradient is zero, we can't make progress
                continue;
            }
            
            // Calculate the projection step size
            let step_size = -worst_violation / constraint_norm_squared;
            
            // Update the solution
            for i in 0..c.len() {
                c[i] = c[i] + step_size * constraint_vector[i];
            }
            
            // Recheck constraints
            constraint_values = constraint_matrix.dot(&c) - constraint_rhs;
        }
        
        // After the iterations, check if all constraints are satisfied
        constraint_values = constraint_matrix.dot(&c) - constraint_rhs;
        all_satisfied = true;
        for &val in constraint_values.iter() {
            if val < -T::from_f64(1e-6).unwrap() {
                all_satisfied = false;
                break;
            }
        }
        
        if !all_satisfied {
            // If we couldn't satisfy all constraints, return an error
            return Err(InterpolateError::ComputationError(
                "Failed to find a solution that satisfies all constraints".to_string(),
            ));
        }
        
        Ok(c)
    }
    
    /// Evaluate the constrained spline at a given point
    ///
    /// # Arguments
    ///
    /// * `x` - The x coordinate at which to evaluate
    ///
    /// # Returns
    ///
    /// The value of the spline at x
    pub fn evaluate(&self, x: T) -> InterpolateResult<T> {
        self.bspline.evaluate(x)
    }
    
    /// Evaluate the constrained spline at multiple points
    ///
    /// # Arguments
    ///
    /// * `x` - Array of x coordinates
    ///
    /// # Returns
    ///
    /// Array of values at the specified x coordinates
    pub fn evaluate_array(&self, x: &ArrayView1<T>) -> InterpolateResult<Array1<T>> {
        self.bspline.evaluate_array(x)
    }
    
    /// Calculate the derivative of the constrained spline at a given point
    ///
    /// # Arguments
    ///
    /// * `x` - The x coordinate
    /// * `order` - The order of the derivative (1 for first derivative, etc.)
    ///
    /// # Returns
    ///
    /// The value of the specified derivative at x
    pub fn derivative(&self, x: T, order: usize) -> InterpolateResult<T> {
        self.bspline.derivative(x, order)
    }
    
    /// Check if the spline satisfies a specific constraint at a point
    ///
    /// # Arguments
    ///
    /// * `constraint_type` - The type of constraint to check
    /// * `x` - The x coordinate at which to check the constraint
    /// * `parameter` - Additional parameter for certain constraint types
    ///
    /// # Returns
    ///
    /// true if the constraint is satisfied, false otherwise
    pub fn check_constraint(
        &self,
        constraint_type: ConstraintType,
        x: T,
        parameter: Option<T>,
    ) -> InterpolateResult<bool> {
        match constraint_type {
            ConstraintType::MonotoneIncreasing => {
                let derivative = self.derivative(x, 1)?;
                Ok(derivative >= -T::epsilon())
            }
            ConstraintType::MonotoneDecreasing => {
                let derivative = self.derivative(x, 1)?;
                Ok(derivative <= T::epsilon())
            }
            ConstraintType::Convex => {
                let second_derivative = self.derivative(x, 2)?;
                Ok(second_derivative >= -T::epsilon())
            }
            ConstraintType::Concave => {
                let second_derivative = self.derivative(x, 2)?;
                Ok(second_derivative <= T::epsilon())
            }
            ConstraintType::Positive => {
                let value = self.evaluate(x)?;
                Ok(value >= -T::epsilon())
            }
            ConstraintType::UpperBound => {
                let value = self.evaluate(x)?;
                let bound = parameter.unwrap_or(T::one());
                Ok(value <= bound + T::epsilon())
            }
            ConstraintType::LowerBound => {
                let value = self.evaluate(x)?;
                let bound = parameter.unwrap_or(T::zero());
                Ok(value >= bound - T::epsilon())
            }
        }
    }
    
    /// Get the list of constraints applied to the spline
    pub fn constraints(&self) -> &[Constraint<T>] {
        &self.constraints
    }
    
    /// Get the underlying B-spline
    pub fn bspline(&self) -> &BSpline<T> {
        &self.bspline
    }
    
    /// Get the fitting method used
    pub fn method(&self) -> FittingMethod {
        self.method
    }
}

/// Convenience function to create a monotonically increasing spline
///
/// # Arguments
///
/// * `x` - The x coordinates of the data points
/// * `y` - The y coordinates of the data points
/// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
/// * `extrapolate` - Extrapolation mode
///
/// # Returns
///
/// A new constrained spline that interpolates the data points while guaranteeing monotonicity
pub fn monotone_increasing_spline<T>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<ConstrainedSpline<T>>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    let constraint = Constraint::monotone_increasing(None, None);
    ConstrainedSpline::interpolate(x, y, vec![constraint], degree, extrapolate)
}

/// Convenience function to create a monotonically decreasing spline
///
/// # Arguments
///
/// * `x` - The x coordinates of the data points
/// * `y` - The y coordinates of the data points
/// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
/// * `extrapolate` - Extrapolation mode
///
/// # Returns
///
/// A new constrained spline that interpolates the data points while guaranteeing monotonic decrease
pub fn monotone_decreasing_spline<T>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<ConstrainedSpline<T>>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    let constraint = Constraint::monotone_decreasing(None, None);
    ConstrainedSpline::interpolate(x, y, vec![constraint], degree, extrapolate)
}

/// Convenience function to create a convex spline
///
/// # Arguments
///
/// * `x` - The x coordinates of the data points
/// * `y` - The y coordinates of the data points
/// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
/// * `extrapolate` - Extrapolation mode
///
/// # Returns
///
/// A new constrained spline that interpolates the data points while guaranteeing convexity
pub fn convex_spline<T>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<ConstrainedSpline<T>>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    let constraint = Constraint::convex(None, None);
    ConstrainedSpline::interpolate(x, y, vec![constraint], degree, extrapolate)
}

/// Convenience function to create a concave spline
///
/// # Arguments
///
/// * `x` - The x coordinates of the data points
/// * `y` - The y coordinates of the data points
/// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
/// * `extrapolate` - Extrapolation mode
///
/// # Returns
///
/// A new constrained spline that interpolates the data points while guaranteeing concavity
pub fn concave_spline<T>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<ConstrainedSpline<T>>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    let constraint = Constraint::concave(None, None);
    ConstrainedSpline::interpolate(x, y, vec![constraint], degree, extrapolate)
}

/// Convenience function to create a positive spline
///
/// # Arguments
///
/// * `x` - The x coordinates of the data points
/// * `y` - The y coordinates of the data points
/// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
/// * `extrapolate` - Extrapolation mode
///
/// # Returns
///
/// A new constrained spline that interpolates the data points while guaranteeing positivity
pub fn positive_spline<T>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<ConstrainedSpline<T>>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    let constraint = Constraint::positive(None, None);
    ConstrainedSpline::interpolate(x, y, vec![constraint], degree, extrapolate)
}

/// Convenience function to create a spline with both monotonicity and convexity constraints
///
/// # Arguments
///
/// * `x` - The x coordinates of the data points
/// * `y` - The y coordinates of the data points
/// * `increasing` - Whether the spline should be monotonically increasing (true) or decreasing (false)
/// * `convex` - Whether the spline should be convex (true) or concave (false)
/// * `degree` - Degree of the B-spline (defaults to 3 for cubic splines)
/// * `extrapolate` - Extrapolation mode
///
/// # Returns
///
/// A new constrained spline with both monotonicity and convexity constraints
pub fn monotone_convex_spline<T>(
    x: &ArrayView1<T>,
    y: &ArrayView1<T>,
    increasing: bool,
    convex: bool,
    degree: usize,
    extrapolate: ExtrapolateMode,
) -> InterpolateResult<ConstrainedSpline<T>>
where
    T: Float + FromPrimitive + Debug + Add<Output = T> + Sub<Output = T> + Mul<Output = T> + Div<Output = T>,
{
    let monotone_constraint = if increasing {
        Constraint::monotone_increasing(None, None)
    } else {
        Constraint::monotone_decreasing(None, None)
    };
    
    let convexity_constraint = if convex {
        Constraint::convex(None, None)
    } else {
        Constraint::concave(None, None)
    };
    
    ConstrainedSpline::interpolate(
        x, 
        y, 
        vec![monotone_constraint, convexity_constraint], 
        degree, 
        extrapolate
    )
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use ndarray::array;

    #[test]
    fn test_monotone_increasing_spline() {
        // Create monotonically increasing data
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![0.0, 0.8, 1.5, 1.9, 2.3, 2.5];
        
        // Create a monotone increasing spline
        let spline = monotone_increasing_spline(
            &x.view(),
            &y.view(),
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline interpolates the data points
        for i in 0..x.len() {
            let y_spline = spline.evaluate(x[i]).unwrap();
            assert_relative_eq!(y_spline, y[i], epsilon = 1e-10);
        }
        
        // Check monotonicity by evaluating at intermediate points
        let test_points = array![0.5, 1.5, 2.5, 3.5, 4.5];
        let mut prev_value = spline.evaluate(0.0).unwrap();
        
        for &x_val in test_points.iter() {
            let y_val = spline.evaluate(x_val).unwrap();
            assert!(y_val >= prev_value, "Spline is not monotonically increasing at x = {}", x_val);
            prev_value = y_val;
        }
        
        // Test the derivative - should be non-negative everywhere
        for &x_val in test_points.iter() {
            let deriv = spline.derivative(x_val, 1).unwrap();
            assert!(deriv >= -1e-10, "Derivative is negative at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_monotone_decreasing_spline() {
        // Create monotonically decreasing data
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![5.0, 4.2, 3.5, 2.7, 1.8, 1.0];
        
        // Create a monotone decreasing spline
        let spline = monotone_decreasing_spline(
            &x.view(),
            &y.view(),
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline interpolates the data points
        for i in 0..x.len() {
            let y_spline = spline.evaluate(x[i]).unwrap();
            assert_relative_eq!(y_spline, y[i], epsilon = 1e-10);
        }
        
        // Check monotonicity by evaluating at intermediate points
        let test_points = array![0.5, 1.5, 2.5, 3.5, 4.5];
        let mut prev_value = spline.evaluate(0.0).unwrap();
        
        for &x_val in test_points.iter() {
            let y_val = spline.evaluate(x_val).unwrap();
            assert!(y_val <= prev_value, "Spline is not monotonically decreasing at x = {}", x_val);
            prev_value = y_val;
        }
        
        // Test the derivative - should be non-positive everywhere
        for &x_val in test_points.iter() {
            let deriv = spline.derivative(x_val, 1).unwrap();
            assert!(deriv <= 1e-10, "Derivative is positive at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_convex_spline() {
        // Create data that should produce a convex function
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![5.0, 3.0, 1.8, 1.0, 0.7, 0.5];
        
        // Create a convex spline
        let spline = convex_spline(
            &x.view(),
            &y.view(),
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline interpolates the data points
        for i in 0..x.len() {
            let y_spline = spline.evaluate(x[i]).unwrap();
            assert_relative_eq!(y_spline, y[i], epsilon = 1e-10);
        }
        
        // Test the second derivative - should be non-negative everywhere for convexity
        let test_points = array![0.5, 1.5, 2.5, 3.5, 4.5];
        
        for &x_val in test_points.iter() {
            let second_deriv = spline.derivative(x_val, 2).unwrap();
            assert!(second_deriv >= -1e-10, "Second derivative is negative at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_concave_spline() {
        // Create data that should produce a concave function
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![0.5, 1.0, 2.2, 4.0, 6.5, 10.0];
        
        // Create a concave spline
        let spline = concave_spline(
            &x.view(),
            &y.view(),
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline interpolates the data points
        for i in 0..x.len() {
            let y_spline = spline.evaluate(x[i]).unwrap();
            assert_relative_eq!(y_spline, y[i], epsilon = 1e-10);
        }
        
        // Test the second derivative - should be non-positive everywhere for concavity
        let test_points = array![0.5, 1.5, 2.5, 3.5, 4.5];
        
        for &x_val in test_points.iter() {
            let second_deriv = spline.derivative(x_val, 2).unwrap();
            assert!(second_deriv <= 1e-10, "Second derivative is positive at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_positive_spline() {
        // Create data that crosses zero
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![1.0, 0.5, 0.1, 0.2, 0.3, 0.5];
        
        // Create a positive spline
        let spline = positive_spline(
            &x.view(),
            &y.view(),
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline interpolates the data points
        for i in 0..x.len() {
            let y_spline = spline.evaluate(x[i]).unwrap();
            assert_relative_eq!(y_spline, y[i], epsilon = 1e-10);
        }
        
        // Test positivity by evaluating at intermediate points
        let test_points = array![0.5, 1.5, 2.5, 3.5, 4.5];
        
        for &x_val in test_points.iter() {
            let y_val = spline.evaluate(x_val).unwrap();
            assert!(y_val >= -1e-10, "Spline is negative at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_multiple_constraints() {
        // Create data that should be monotonically increasing and convex
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![0.0, 0.1, 0.4, 0.9, 1.6, 2.5];
        
        // Create a spline with both constraints
        let spline = monotone_convex_spline(
            &x.view(),
            &y.view(),
            true, // increasing
            true, // convex
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline interpolates the data points
        for i in 0..x.len() {
            let y_spline = spline.evaluate(x[i]).unwrap();
            assert_relative_eq!(y_spline, y[i], epsilon = 1e-10);
        }
        
        // Test constraints at intermediate points
        let test_points = array![0.5, 1.5, 2.5, 3.5, 4.5];
        let mut prev_value = spline.evaluate(0.0).unwrap();
        
        for &x_val in test_points.iter() {
            let y_val = spline.evaluate(x_val).unwrap();
            let first_deriv = spline.derivative(x_val, 1).unwrap();
            let second_deriv = spline.derivative(x_val, 2).unwrap();
            
            // Test monotonicity
            assert!(y_val >= prev_value, "Spline is not monotonically increasing at x = {}", x_val);
            assert!(first_deriv >= -1e-10, "Derivative is negative at x = {}", x_val);
            
            // Test convexity
            assert!(second_deriv >= -1e-10, "Second derivative is negative at x = {}", x_val);
            
            prev_value = y_val;
        }
    }
    
    #[test]
    fn test_least_squares_fitting() {
        // Create noisy data
        let x = array![0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0];
        let mut y = Array1::zeros(x.len());
        
        // Generate y = x^2 + noise
        for i in 0..x.len() {
            y[i] = x[i] * x[i] + 0.2 * (i as f64 * 0.7).sin();
        }
        
        // Create a least squares spline with convexity constraint
        let spline = ConstrainedSpline::least_squares(
            &x.view(),
            &y.view(),
            vec![Constraint::convex(None, None)],
            10, // knots
            3,  // cubic
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // The spline should be convex
        let test_points = array![0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75];
        
        for &x_val in test_points.iter() {
            let second_deriv = spline.derivative(x_val, 2).unwrap();
            assert!(second_deriv >= -1e-10, "Second derivative is negative at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_penalized_fitting() {
        // Create noisy data
        let x = array![0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0];
        let mut y = Array1::zeros(x.len());
        
        // Generate y = x^2 + noise
        for i in 0..x.len() {
            y[i] = x[i] * x[i] + 0.2 * (i as f64 * 0.7).sin();
        }
        
        // Create a penalized spline with convexity constraint
        let spline = ConstrainedSpline::penalized(
            &x.view(),
            &y.view(),
            vec![Constraint::convex(None, None)],
            10, // knots
            3,  // cubic
            0.1, // lambda
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // The spline should be convex
        let test_points = array![0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75];
        
        for &x_val in test_points.iter() {
            let second_deriv = spline.derivative(x_val, 2).unwrap();
            assert!(second_deriv >= -1e-10, "Second derivative is negative at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_bound_constraints() {
        // Create some data
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![0.5, 0.8, 1.2, 0.7, 0.3, 0.6];
        
        // Create a spline with upper and lower bounds
        let upper_bound = 1.0;
        let lower_bound = 0.3;
        
        let constraints = vec![
            Constraint::upper_bound(None, None, upper_bound),
            Constraint::lower_bound(None, None, lower_bound),
        ];
        
        let spline = ConstrainedSpline::interpolate(
            &x.view(),
            &y.view(),
            constraints,
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check that the spline stays within the bounds
        let test_points = array![
            0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0,
        ];
        
        for &x_val in test_points.iter() {
            let y_val = spline.evaluate(x_val).unwrap();
            assert!(y_val <= upper_bound + 1e-10, "Spline exceeds upper bound at x = {}", x_val);
            assert!(y_val >= lower_bound - 1e-10, "Spline is below lower bound at x = {}", x_val);
        }
    }
    
    #[test]
    fn test_constraint_checking() {
        // Create monotonically increasing data
        let x = array![0.0, 1.0, 2.0, 3.0, 4.0, 5.0];
        let y = array![0.0, 0.8, 1.5, 1.9, 2.3, 2.5];
        
        // Create a monotone increasing spline
        let spline = monotone_increasing_spline(
            &x.view(),
            &y.view(),
            3,
            ExtrapolateMode::Extrapolate,
        ).unwrap();
        
        // Check the monotonicity constraint
        let test_point = 2.5;
        
        let increasing_satisfied = spline.check_constraint(
            ConstraintType::MonotoneIncreasing,
            test_point,
            None,
        ).unwrap();
        
        assert!(increasing_satisfied);
        
        // Check a constraint that shouldn't be satisfied
        let decreasing_satisfied = spline.check_constraint(
            ConstraintType::MonotoneDecreasing,
            test_point,
            None,
        ).unwrap();
        
        assert!(!decreasing_satisfied);
    }
}