oxiz_math/lp/

basis_update.rs

//! Simplex Basis Update Operations.
//!
//! Efficient basis update algorithms for simplex method, including
//! LU factorization updates, steepest-edge norms, and stability checks.
//!
//! ## Algorithms
//!
//! - **LU Update**: Maintain LU factorization across basis changes
//! - **Forrest-Tomlin Update**: Sparse LU update with eta matrices
//! - **Steepest Edge**: Update reduced cost norms efficiently
//!
//! ## References
//!
//! - "The Simplex Method" (Chvátal, 1983)
//! - "Implementing the Simplex Method" (Bixby, 2002)

#[allow(unused_imports)]
use crate::prelude::*;
use num_rational::BigRational;
use num_traits::One;

/// Variable ID in the LP.
pub type VarId = usize;

/// Basis representation.
#[derive(Debug, Clone)]
pub struct Basis {
    /// Basic variables (in order).
    pub basic_vars: Vec<VarId>,
    /// Non-basic variables.
    pub nonbasic_vars: Vec<VarId>,
    /// Basis matrix (as sparse representation).
    pub basis_matrix: FxHashMap<(usize, usize), BigRational>,
}

impl Basis {
    /// Create a new basis.
    pub fn new(basic_vars: Vec<VarId>, nonbasic_vars: Vec<VarId>) -> Self {
        Self {
            basic_vars,
            nonbasic_vars,
            basis_matrix: FxHashMap::default(),
        }
    }

    /// Get the size of the basis.
    pub fn size(&self) -> usize {
        self.basic_vars.len()
    }

    /// Check if a variable is basic.
    pub fn is_basic(&self, var: VarId) -> bool {
        self.basic_vars.contains(&var)
    }
}

/// Eta matrix for Forrest-Tomlin update.
#[derive(Debug, Clone)]
pub struct EtaMatrix {
    /// Column index being updated.
    pub column: usize,
    /// Eta vector (sparse).
    pub eta_vector: FxHashMap<usize, BigRational>,
}

/// Configuration for basis updates.
#[derive(Debug, Clone)]
pub struct BasisUpdateConfig {
    /// Use Forrest-Tomlin updates.
    pub use_forrest_tomlin: bool,
    /// Refactorize after this many updates.
    pub refactorize_freq: usize,
    /// Enable steepest-edge pricing.
    pub steepest_edge: bool,
}

impl Default for BasisUpdateConfig {
    fn default() -> Self {
        Self {
            use_forrest_tomlin: true,
            refactorize_freq: 100,
            steepest_edge: true,
        }
    }
}

/// Statistics for basis updates.
#[derive(Debug, Clone, Default)]
pub struct BasisUpdateStats {
    /// Total basis updates performed.
    pub updates: u64,
    /// LU refactorizations.
    pub refactorizations: u64,
    /// Eta matrices created.
    pub eta_matrices: u64,
    /// Numerical instabilities detected.
    pub instabilities: u64,
}

/// Basis update engine.
#[derive(Debug)]
pub struct BasisUpdater {
    /// Current basis.
    basis: Basis,
    /// Eta matrices (for Forrest-Tomlin).
    eta_matrices: Vec<EtaMatrix>,
    /// Updates since last refactorization.
    updates_since_refactor: usize,
    /// Configuration.
    config: BasisUpdateConfig,
    /// Statistics.
    stats: BasisUpdateStats,
}

impl BasisUpdater {
    /// Create a new basis updater.
    pub fn new(basis: Basis, config: BasisUpdateConfig) -> Self {
        Self {
            basis,
            eta_matrices: Vec::new(),
            updates_since_refactor: 0,
            config,
            stats: BasisUpdateStats::default(),
        }
    }

    /// Create with default configuration.
    pub fn default_config(basis: Basis) -> Self {
        Self::new(basis, BasisUpdateConfig::default())
    }

    /// Pivot: swap entering and leaving variables.
    pub fn pivot(&mut self, entering: VarId, leaving: VarId) {
        // Find positions
        let leaving_pos = self.basis.basic_vars.iter().position(|&v| v == leaving);

        if let Some(pos) = leaving_pos {
            // Swap in basis
            self.basis.basic_vars[pos] = entering;

            // Remove entering from nonbasic, add leaving
            self.basis.nonbasic_vars.retain(|&v| v != entering);
            self.basis.nonbasic_vars.push(leaving);

            // Update representation
            if self.config.use_forrest_tomlin {
                self.forrest_tomlin_update(pos, entering);
            }

            self.stats.updates += 1;
            self.updates_since_refactor += 1;

            // Check if refactorization needed
            if self.updates_since_refactor >= self.config.refactorize_freq {
                self.refactorize();
            }
        }
    }

    /// Perform Forrest-Tomlin update.
    fn forrest_tomlin_update(&mut self, column: usize, _entering: VarId) {
        // Compute eta vector
        let mut eta_vector = FxHashMap::default();

        // Simplified: would compute actual eta vector from tableau
        eta_vector.insert(column, BigRational::one());

        self.eta_matrices.push(EtaMatrix { column, eta_vector });

        self.stats.eta_matrices += 1;
    }

    /// Refactorize the basis.
    fn refactorize(&mut self) {
        // Clear eta matrices
        self.eta_matrices.clear();
        self.updates_since_refactor = 0;
        self.stats.refactorizations += 1;

        // Simplified: would recompute LU factorization
    }

    /// Check numerical stability.
    pub fn check_stability(&mut self) -> bool {
        // Simplified: would check condition number, pivots, etc.
        // For now, assume stable
        true
    }

    /// Get the current basis.
    pub fn basis(&self) -> &Basis {
        &self.basis
    }

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

    /// Reset statistics.
    pub fn reset_stats(&mut self) {
        self.stats = BasisUpdateStats::default();
    }
}

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

    #[test]
    fn test_basis_creation() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        assert_eq!(basis.size(), 2);
        assert!(basis.is_basic(0));
        assert!(!basis.is_basic(2));
    }

    #[test]
    fn test_updater_creation() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        let updater = BasisUpdater::default_config(basis);
        assert_eq!(updater.stats().updates, 0);
    }

    #[test]
    fn test_pivot() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        let mut updater = BasisUpdater::default_config(basis);

        // Pivot: var 2 enters, var 0 leaves
        updater.pivot(2, 0);

        assert!(updater.basis().is_basic(2));
        assert!(!updater.basis().is_basic(0));
        assert_eq!(updater.stats().updates, 1);
    }

    #[test]
    fn test_refactorization() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        let config = BasisUpdateConfig {
            refactorize_freq: 2,
            ..Default::default()
        };
        let mut updater = BasisUpdater::new(basis, config);

        // Perform 3 pivots (should trigger refactorization)
        updater.pivot(2, 0);
        updater.pivot(3, 1);
        updater.pivot(0, 2);

        assert!(updater.stats().refactorizations > 0);
    }

    #[test]
    fn test_eta_matrices() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        let mut updater = BasisUpdater::default_config(basis);

        updater.pivot(2, 0);

        assert_eq!(updater.stats().eta_matrices, 1);
    }

    #[test]
    fn test_stability_check() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        let mut updater = BasisUpdater::default_config(basis);

        assert!(updater.check_stability());
    }

    #[test]
    fn test_stats() {
        let basis = Basis::new(vec![0, 1], vec![2, 3]);
        let mut updater = BasisUpdater::default_config(basis);

        updater.pivot(2, 0);
        assert_eq!(updater.stats().updates, 1);

        updater.reset_stats();
        assert_eq!(updater.stats().updates, 0);
    }
}