oxiz-math 0.3.1

//! Cutting Plane Methods for Mixed Integer Programming.
//!
//! Implements various cutting plane techniques including:
//! - Gomory mixed-integer cuts
//! - Lift-and-project cuts
//! - Cover inequalities

#[allow(unused_imports)]
use crate::prelude::*;
use num_bigint::BigInt;
use num_rational::BigRational;
use num_traits::{One, Zero};

/// Cutting plane generator for MIP.
#[derive(Debug)]
pub struct CuttingPlaneGenerator {
    /// Integer variables
    integer_vars: FxHashSet<VarId>,
    /// Statistics
    stats: CuttingPlaneStats,
}

/// Variable identifier
pub type VarId = usize;

/// A cutting plane (linear inequality).
#[derive(Debug, Clone)]
pub struct CuttingPlane {
    /// Coefficients
    pub coeffs: Vec<(VarId, BigRational)>,
    /// Right-hand side
    pub rhs: BigRational,
    /// Type of cut
    pub cut_type: CutType,
}

/// Type of cutting plane.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum CutType {
    /// Gomory fractional cut
    Gomory,
    /// Gomory mixed-integer cut
    GomoryMI,
    /// Lift-and-project cut
    LiftProject,
    /// Cover inequality
    Cover,
    /// Clique inequality
    Clique,
}

/// Configuration for cutting plane generation.
#[derive(Debug, Clone)]
pub struct CuttingPlaneConfig {
    /// Maximum cuts to generate
    pub max_cuts: usize,
    /// Enable Gomory cuts
    pub enable_gomory: bool,
    /// Enable lift-and-project cuts
    pub enable_lift_project: bool,
    /// Enable cover cuts
    pub enable_cover: bool,
}

impl Default for CuttingPlaneConfig {
    fn default() -> Self {
        Self {
            max_cuts: 100,
            enable_gomory: true,
            enable_lift_project: true,
            enable_cover: true,
        }
    }
}

/// Cutting plane statistics.
#[derive(Debug, Clone, Default)]
pub struct CuttingPlaneStats {
    /// Total cuts generated
    pub cuts_generated: usize,
    /// Gomory cuts
    pub gomory_cuts: usize,
    /// Lift-and-project cuts
    pub lift_project_cuts: usize,
    /// Cover cuts
    pub cover_cuts: usize,
}

impl CuttingPlaneGenerator {
    /// Create a new cutting plane generator.
    pub fn new(integer_vars: FxHashSet<VarId>) -> Self {
        Self {
            integer_vars,
            stats: CuttingPlaneStats::default(),
        }
    }

    /// Generate Gomory fractional cut from a tableau row.
    pub fn generate_gomory_cut(
        &mut self,
        _basic_var: VarId,
        row: &[(VarId, BigRational)],
        rhs: &BigRational,
    ) -> Option<CuttingPlane> {
        // Extract fractional part of RHS
        let frac_rhs = self.fractional_part(rhs);
        if frac_rhs.is_zero() {
            return None; // No fractional part, no cut
        }

        let mut cut_coeffs = Vec::new();

        for (var_id, coeff) in row {
            let frac_coeff = self.fractional_part(coeff);

            if !frac_coeff.is_zero() {
                // Cut coefficient: -f_j if f_j <= f_0, else -(1-f_j)
                let cut_coeff = if frac_coeff <= frac_rhs {
                    -frac_coeff
                } else {
                    -(BigRational::one() - frac_coeff)
                };

                cut_coeffs.push((*var_id, cut_coeff));
            }
        }

        self.stats.gomory_cuts += 1;
        self.stats.cuts_generated += 1;

        Some(CuttingPlane {
            coeffs: cut_coeffs,
            rhs: -frac_rhs,
            cut_type: CutType::Gomory,
        })
    }

    /// Generate Gomory mixed-integer cut.
    pub fn generate_gomory_mi_cut(
        &mut self,
        basic_var: VarId,
        row: &[(VarId, BigRational)],
        rhs: &BigRational,
    ) -> Option<CuttingPlane> {
        if !self.integer_vars.contains(&basic_var) {
            return None; // Only for integer basic variables
        }

        let f0 = self.fractional_part(rhs);
        if f0.is_zero() {
            return None;
        }

        let mut cut_coeffs = Vec::new();

        for (var_id, coeff) in row {
            let fj = self.fractional_part(coeff);

            let cut_coeff = if self.integer_vars.contains(var_id) {
                // Integer variable
                if fj <= f0 {
                    -fj
                } else {
                    -(BigRational::one() - fj) * &f0 / (BigRational::one() - &f0)
                }
            } else {
                // Continuous variable
                if coeff >= &BigRational::zero() {
                    -coeff
                } else {
                    coeff * &f0 / (BigRational::one() - &f0)
                }
            };

            if !cut_coeff.is_zero() {
                cut_coeffs.push((*var_id, cut_coeff));
            }
        }

        self.stats.gomory_cuts += 1;
        self.stats.cuts_generated += 1;

        Some(CuttingPlane {
            coeffs: cut_coeffs,
            rhs: -f0,
            cut_type: CutType::GomoryMI,
        })
    }

    /// Generate lift-and-project cut.
    pub fn generate_lift_project_cut(
        &mut self,
        constraint: &[(VarId, BigRational)],
        rhs: &BigRational,
        lifting_var: VarId,
    ) -> Option<CuttingPlane> {
        // Simplified lift-and-project
        if !self.integer_vars.contains(&lifting_var) {
            return None;
        }

        // Generate disjunctive cut: x_j = 0 or x_j = 1
        let mut cut_coeffs = Vec::new();

        for (var_id, coeff) in constraint {
            if *var_id == lifting_var {
                continue;
            }

            // Lift coefficient
            let lifted_coeff = coeff.clone();
            cut_coeffs.push((*var_id, lifted_coeff));
        }

        self.stats.lift_project_cuts += 1;
        self.stats.cuts_generated += 1;

        Some(CuttingPlane {
            coeffs: cut_coeffs,
            rhs: rhs.clone(),
            cut_type: CutType::LiftProject,
        })
    }

    /// Generate cover inequality from a knapsack constraint.
    pub fn generate_cover_cut(
        &mut self,
        weights: &[(VarId, BigRational)],
        capacity: &BigRational,
    ) -> Option<CuttingPlane> {
        // Find minimal cover
        let cover = self.find_minimal_cover(weights, capacity)?;

        // Generate cover inequality: sum_{j in C} x_j <= |C| - 1
        let mut cut_coeffs = Vec::new();
        for var_id in &cover {
            cut_coeffs.push((*var_id, BigRational::one()));
        }

        let rhs = BigRational::from_integer(BigInt::from(cover.len() as i64 - 1));

        self.stats.cover_cuts += 1;
        self.stats.cuts_generated += 1;

        Some(CuttingPlane {
            coeffs: cut_coeffs,
            rhs,
            cut_type: CutType::Cover,
        })
    }

    /// Find minimal cover for knapsack.
    fn find_minimal_cover(
        &self,
        weights: &[(VarId, BigRational)],
        capacity: &BigRational,
    ) -> Option<Vec<VarId>> {
        // Greedy approach: sort by weight descending
        let mut sorted_weights = weights.to_vec();
        sorted_weights.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(core::cmp::Ordering::Equal));

        let mut cover = Vec::new();
        let mut total_weight = BigRational::zero();

        for (var_id, weight) in sorted_weights {
            cover.push(var_id);
            total_weight = &total_weight + &weight;

            if total_weight > *capacity {
                return Some(cover);
            }
        }

        None // No cover found
    }

    /// Extract fractional part of a rational number.
    fn fractional_part(&self, value: &BigRational) -> BigRational {
        // frac(x) = x - floor(x)
        let floor_val = self.floor(value);
        value - &floor_val
    }

    /// Floor function for rationals (true mathematical floor, correct for
    /// negative values: e.g. floor(-7/3) == -3, not -2).
    fn floor(&self, value: &BigRational) -> BigRational {
        BigRational::from_integer(crate::rational::floor(value))
    }

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

    /// Check if a cut is violated by current solution.
    pub fn is_violated(&self, cut: &CuttingPlane, solution: &[(VarId, BigRational)]) -> bool {
        let mut lhs = BigRational::zero();

        for (var_id, coeff) in &cut.coeffs {
            if let Some((_, value)) = solution.iter().find(|(id, _)| id == var_id) {
                lhs += coeff * value;
            }
        }

        lhs > cut.rhs
    }
}

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

    #[test]
    fn test_cutting_plane_generator() {
        let mut integer_vars = FxHashSet::default();
        integer_vars.insert(0);
        integer_vars.insert(1);

        let generator = CuttingPlaneGenerator::new(integer_vars);
        assert_eq!(generator.stats.cuts_generated, 0);
    }

    #[test]
    fn test_fractional_part() {
        let integer_vars = FxHashSet::default();
        let generator = CuttingPlaneGenerator::new(integer_vars);

        let val = BigRational::new(BigInt::from(7), BigInt::from(3)); // 7/3 = 2.333...
        let frac = generator.fractional_part(&val);

        // Should be 1/3
        assert_eq!(frac, BigRational::new(BigInt::from(1), BigInt::from(3)));
    }

    #[test]
    fn test_floor_negative_fraction() {
        // Regression test for MATH-1: floor must round toward negative
        // infinity, not truncate toward zero.
        let integer_vars = FxHashSet::default();
        let generator = CuttingPlaneGenerator::new(integer_vars);

        let val = BigRational::new(BigInt::from(-7), BigInt::from(3)); // -7/3 = -2.333...
        let floor_val = generator.floor(&val);
        assert_eq!(floor_val, BigRational::from_integer(BigInt::from(-3)));
    }

    #[test]
    fn test_fractional_part_negative() {
        // Regression test for MATH-1: fractional_part(-7/3) must be in
        // [0, 1), i.e. 2/3, not the negative value -1/3 produced by a
        // truncating floor.
        let integer_vars = FxHashSet::default();
        let generator = CuttingPlaneGenerator::new(integer_vars);

        let val = BigRational::new(BigInt::from(-7), BigInt::from(3));
        let frac = generator.fractional_part(&val);

        assert!(frac >= BigRational::zero());
        assert!(frac < BigRational::one());
        assert_eq!(frac, BigRational::new(BigInt::from(2), BigInt::from(3)));
    }

    #[test]
    fn test_fractional_part_negative_integer_is_zero() {
        let integer_vars = FxHashSet::default();
        let generator = CuttingPlaneGenerator::new(integer_vars);

        let val = BigRational::from_integer(BigInt::from(-5));
        let frac = generator.fractional_part(&val);
        assert!(frac.is_zero());
    }

    #[test]
    fn test_gomory_cut_negative_rhs_produces_valid_fraction() {
        // Ensure generate_gomory_cut, which relies on fractional_part being
        // in [0, 1), behaves correctly for negative RHS/coefficients.
        let mut integer_vars = FxHashSet::default();
        integer_vars.insert(0);
        let mut generator = CuttingPlaneGenerator::new(integer_vars);

        let row = vec![(0usize, BigRational::new(BigInt::from(-7), BigInt::from(3)))];
        let rhs = BigRational::new(BigInt::from(-7), BigInt::from(3));

        let cut = generator
            .generate_gomory_cut(1, &row, &rhs)
            .expect("expected a cut for a fractional RHS");

        // rhs of the cut is -frac_rhs, which must lie in (-1, 0].
        assert!(cut.rhs <= BigRational::zero());
        assert!(cut.rhs > -BigRational::one());
    }
}