oxiz_math/

simplex_parametric.rs

//! Parametric Simplex for Sensitivity Analysis.
#![allow(dead_code, missing_docs)] // Under development
//!
//! Extends simplex to solve parametric linear programs where objective
//! coefficients or RHS values depend on a parameter λ.
//!
//! ## Applications
//!
//! - Sensitivity analysis: how does optimal value change with parameter?
//! - Optimal control: parametric objectives
//! - Theory propagation: bounds propagation as parameter varies
//!
//! ## Algorithm
//!
//! - Solve for critical parameter values where basis changes
//! - Construct piecewise-linear optimal value function
//! - Report optimal solutions for each parameter interval
//!
//! ## References
//!
//! - Gass & Saaty: "The Computational Algorithm for the Parametric Objective Function"
//! - Z3's `math/lp/parametric_simplex.cpp`

// TODO: Re-enable after SimplexSolver API is available
// use crate::simplex::SimplexSolver;
#[allow(unused_imports)]
use crate::prelude::*;
use num_rational::BigRational;
use num_traits::{One, Zero};

/// Variable identifier.
pub type VarId = usize;

/// Parameter type (which coefficient/RHS is parametric).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ParametricType {
    /// Objective coefficient is parametric: c_j(λ) = c_j + λ * d_j
    Objective,
    /// RHS is parametric: b_i(λ) = b_i + λ * e_i
    Rhs,
}

/// A breakpoint in the parametric solution.
#[derive(Debug, Clone)]
pub struct Breakpoint {
    /// Parameter value at breakpoint.
    pub lambda: BigRational,
    /// Optimal value at this breakpoint.
    pub optimal_value: BigRational,
    /// Basis at this breakpoint.
    pub basis: Vec<VarId>,
}

/// Parametric interval with constant basis.
#[derive(Debug, Clone)]
pub struct ParametricInterval {
    /// Lower bound on λ (inclusive).
    pub lambda_min: BigRational,
    /// Upper bound on λ (exclusive).
    pub lambda_max: BigRational,
    /// Optimal value function: z(λ) = a + b*λ
    pub value_slope: BigRational,
    pub value_intercept: BigRational,
    /// Basis for this interval.
    pub basis: Vec<VarId>,
}

/// Configuration for parametric simplex.
#[derive(Debug, Clone)]
pub struct ParametricSimplexConfig {
    /// Parametric type.
    pub param_type: ParametricType,
    /// Maximum number of breakpoints to compute.
    pub max_breakpoints: usize,
    /// Lambda range to explore.
    pub lambda_min: BigRational,
    pub lambda_max: BigRational,
}

impl Default for ParametricSimplexConfig {
    fn default() -> Self {
        Self {
            param_type: ParametricType::Objective,
            max_breakpoints: 1000,
            lambda_min: BigRational::zero(),
            lambda_max: BigRational::from_integer(100.into()),
        }
    }
}

/// Statistics for parametric simplex.
#[derive(Debug, Clone, Default)]
pub struct ParametricSimplexStats {
    /// Number of breakpoints computed.
    pub breakpoints: u64,
    /// Number of intervals.
    pub intervals: u64,
    /// Simplex iterations.
    pub iterations: u64,
}

/// Result of parametric simplex.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum ParametricSimplexResult {
    /// Solution computed successfully.
    Success,
    /// Problem is infeasible for all λ.
    Infeasible,
    /// Problem is unbounded for some λ.
    Unbounded,
    /// Reached breakpoint limit.
    BreakpointLimit,
}

/// Parametric simplex solver.
#[derive(Debug)]
pub struct ParametricSimplexSolver {
    /// Configuration.
    config: ParametricSimplexConfig,
    /// TODO: Base simplex solver (disabled until SimplexSolver API available).
    // simplex: SimplexSolver,
    /// Parametric coefficients (d_j for objective, e_i for RHS).
    parametric_coeffs: FxHashMap<usize, BigRational>,
    /// Breakpoints found.
    breakpoints: Vec<Breakpoint>,
    /// Intervals with constant basis.
    intervals: Vec<ParametricInterval>,
    /// Statistics.
    stats: ParametricSimplexStats,
}

impl ParametricSimplexSolver {
    /// Create a new parametric simplex solver.
    ///
    /// TODO: Disabled until SimplexSolver API is available.
    pub fn new(config: ParametricSimplexConfig) -> Self {
        Self {
            config,
            // simplex,
            parametric_coeffs: FxHashMap::default(),
            breakpoints: Vec::new(),
            intervals: Vec::new(),
            stats: ParametricSimplexStats::default(),
        }
    }

    /// Create with default configuration.
    pub fn default_config() -> Self {
        Self::new(ParametricSimplexConfig::default())
    }

    /// Set parametric coefficient for a variable/constraint.
    pub fn set_parametric_coeff(&mut self, id: usize, coeff: BigRational) {
        self.parametric_coeffs.insert(id, coeff);
    }

    /// Solve the parametric LP.
    ///
    /// Computes breakpoints and intervals where basis changes.
    pub fn solve(&mut self) -> ParametricSimplexResult {
        // Start at λ = lambda_min
        let _current_lambda = self.config.lambda_min.clone();

        // TODO: Re-enable after SimplexSolver API is available
        // Update simplex with current parameter value
        // self.update_simplex_at_lambda(&current_lambda);

        // Solve initial problem
        // if !self.simplex.solve() {
        //     return ParametricSimplexResult::Infeasible;
        // }

        // Placeholder: return Infeasible until SimplexSolver is available
        ParametricSimplexResult::Infeasible

        // TODO: Re-enable after SimplexSolver API is available
        // Get initial basis
        // let mut current_basis = self.get_current_basis();

        // Compute initial breakpoint
        // let initial_breakpoint = Breakpoint {
        //     lambda: current_lambda.clone(),
        //     optimal_value: self.simplex.objective_value(),
        //     basis: current_basis.clone(),
        // };
        // self.breakpoints.push(initial_breakpoint);
        // self.stats.breakpoints += 1;

        /*
        // Iterate to find breakpoints
        while current_lambda < self.config.lambda_max {
            if self.breakpoints.len() >= self.config.max_breakpoints {
                return ParametricSimplexResult::BreakpointLimit;
            }

            // Compute next breakpoint
            let next_lambda = self.compute_next_breakpoint(&current_lambda);

            if next_lambda >= self.config.lambda_max {
                // Create final interval
                let interval = self.create_interval(
                    current_lambda.clone(),
                    self.config.lambda_max.clone(),
                    &current_basis,
                );
                self.intervals.push(interval);
                self.stats.intervals += 1;
                break;
            }

            // Create interval [current_lambda, next_lambda)
            let interval = self.create_interval(
                current_lambda.clone(),
                next_lambda.clone(),
                &current_basis,
            );
            self.intervals.push(interval);
            self.stats.intervals += 1;

            // TODO: Re-enable after SimplexSolver API is available
            // Update to next lambda
            // current_lambda = next_lambda.clone();
            // self.update_simplex_at_lambda(&current_lambda);

            // Re-optimize with new parameter value
            // if !self.simplex.solve() {
            //     return ParametricSimplexResult::Infeasible;
            // }

            // Update basis
            // current_basis = self.get_current_basis();

            // TODO: Re-enable after SimplexSolver API is available
            // Record breakpoint
            // let breakpoint = Breakpoint {
            //     lambda: current_lambda.clone(),
            //     optimal_value: self.simplex.objective_value(),
            //     basis: current_basis.clone(),
            // };
            // self.breakpoints.push(breakpoint);
            // self.stats.breakpoints += 1;
        }

        ParametricSimplexResult::Success
        */
    }

    /// Update simplex with parameter value λ.
    ///
    /// TODO: Re-enable after SimplexSolver API is available.
    fn update_simplex_at_lambda(&mut self, _lambda: &BigRational) {
        // match self.config.param_type {
        //     ParametricType::Objective => {
        //         // Update objective: c_j(λ) = c_j + λ * d_j
        //         for (&var, d_j) in &self.parametric_coeffs {
        //             let base_coeff = self.simplex.get_objective_coeff(var);
        //             let new_coeff = base_coeff + lambda * d_j;
        //             self.simplex.set_objective_coeff(var, new_coeff);
        //         }
        //     }
        //     ParametricType::Rhs => {
        //         // Update RHS: b_i(λ) = b_i + λ * e_i
        //         for (&constraint, e_i) in &self.parametric_coeffs {
        //             let base_rhs = self.simplex.get_rhs(constraint);
        //             let new_rhs = base_rhs + lambda * e_i;
        //             self.simplex.set_rhs(constraint, new_rhs);
        //         }
        //     }
        // }
    }

    /// Compute next breakpoint where basis changes.
    fn compute_next_breakpoint(&self, current_lambda: &BigRational) -> BigRational {
        // Simplified: advance by fixed step
        // Full implementation would compute exact breakpoint from dual feasibility
        current_lambda.clone() + BigRational::one()
    }

    /// Get current basis from simplex.
    fn get_current_basis(&self) -> Vec<VarId> {
        // Placeholder: extract basis variables from simplex
        Vec::new()
    }

    /// Create interval with given bounds and basis.
    ///
    /// TODO: Re-enable after SimplexSolver API is available.
    fn create_interval(
        &self,
        lambda_min: BigRational,
        lambda_max: BigRational,
        basis: &[VarId],
    ) -> ParametricInterval {
        // TODO: Compute slope and intercept for z(λ) = a + b*λ
        // let value_at_min = self.simplex.objective_value();

        // Simplified: assume constant objective in interval
        let slope = BigRational::zero();
        let intercept = BigRational::zero(); // was value_at_min

        ParametricInterval {
            lambda_min,
            lambda_max,
            value_slope: slope,
            value_intercept: intercept,
            basis: basis.to_vec(),
        }
    }

    /// Get breakpoints.
    pub fn breakpoints(&self) -> &[Breakpoint] {
        &self.breakpoints
    }

    /// Get intervals.
    pub fn intervals(&self) -> &[ParametricInterval] {
        &self.intervals
    }

    /// Evaluate optimal value at parameter λ.
    pub fn evaluate(&self, lambda: &BigRational) -> Option<BigRational> {
        // Find interval containing λ
        for interval in &self.intervals {
            if lambda >= &interval.lambda_min && lambda < &interval.lambda_max {
                // z(λ) = intercept + slope * λ
                let value = interval.value_intercept.clone() + &interval.value_slope * lambda;
                return Some(value);
            }
        }

        None
    }

    /// Get statistics.
    pub fn stats(&self) -> &ParametricSimplexStats {
        &self.stats
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_parametric_creation() {
        // TODO: Re-enable after SimplexSolver API is available
        // let simplex = SimplexSolver::new();
        let solver = ParametricSimplexSolver::default_config();

        assert_eq!(solver.breakpoints().len(), 0);
    }

    #[test]
    fn test_set_parametric_coeff() {
        // TODO: Re-enable after SimplexSolver API is available
        // let simplex = SimplexSolver::new();
        let mut solver = ParametricSimplexSolver::default_config();

        solver.set_parametric_coeff(0, BigRational::new(2.into(), 1.into()));

        assert!(solver.parametric_coeffs.contains_key(&0));
    }
}