oxiz-math 0.3.1

//! Branch-and-Cut Algorithm for Mixed-Integer Programming.
//!
//! Integrates branch-and-bound with cutting planes for solving MIP:
//! - LP relaxation at each node
//! - Cutting plane separation
//! - Variable branching
//! - Node selection strategies
//!
//! ## References
//!
//! - Wolsey (1998): "Integer Programming"
//! - CPLEX and Gurobi MIP solvers

use super::cutting_planes::{CuttingPlaneConfig, CuttingPlaneGenerator};
use super::dual_simplex::{DualSimplexResult, DualSimplexSolver as DualSimplex};
#[allow(unused_imports)]
use crate::prelude::*;
use core::cmp::Ordering;
use num_rational::BigRational;
use num_traits::{Signed, ToPrimitive, Zero};

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

/// Variable type.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum VarType {
    /// Continuous variable.
    Continuous,
    /// Integer variable.
    Integer,
    /// Binary (0/1) variable.
    Binary,
}

/// Node in the branch-and-cut tree.
#[derive(Debug, Clone)]
struct BranchNode {
    /// Node ID.
    id: usize,
    /// LP bound at this node.
    lp_bound: BigRational,
    /// Variable bounds (lower, upper).
    bounds: FxHashMap<VarId, (BigRational, BigRational)>,
    /// Depth in the tree.
    depth: usize,
}

impl PartialEq for BranchNode {
    fn eq(&self, other: &Self) -> bool {
        self.id == other.id
    }
}

impl Eq for BranchNode {}

impl PartialOrd for BranchNode {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for BranchNode {
    fn cmp(&self, other: &Self) -> Ordering {
        // Best-first search: prioritize nodes with better LP bounds
        self.lp_bound.cmp(&other.lp_bound)
    }
}

/// Node selection strategy.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum NodeSelection {
    /// Best-first (best LP bound).
    BestFirst,
    /// Depth-first.
    DepthFirst,
    /// Breadth-first.
    BreadthFirst,
}

/// Branching strategy.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum BranchingStrategy {
    /// Most fractional variable.
    MostFractional,
    /// Least fractional variable.
    LeastFractional,
    /// Largest coefficient in objective.
    LargestCoeff,
}

/// Configuration for branch-and-cut.
#[derive(Debug, Clone)]
pub struct BranchCutConfig {
    /// Node selection strategy.
    pub node_selection: NodeSelection,
    /// Branching strategy.
    pub branching_strategy: BranchingStrategy,
    /// Maximum number of nodes to explore.
    pub max_nodes: usize,
    /// Optimality gap tolerance (relative).
    pub gap_tolerance: f64,
    /// Enable cutting planes.
    pub enable_cuts: bool,
    /// Maximum cuts per node.
    pub max_cuts_per_node: usize,
    /// Cutting plane configuration.
    pub cut_config: CuttingPlaneConfig,
}

impl Default for BranchCutConfig {
    fn default() -> Self {
        Self {
            node_selection: NodeSelection::BestFirst,
            branching_strategy: BranchingStrategy::MostFractional,
            max_nodes: 100_000,
            gap_tolerance: 1e-6,
            enable_cuts: true,
            max_cuts_per_node: 10,
            cut_config: CuttingPlaneConfig::default(),
        }
    }
}

/// Statistics for branch-and-cut.
#[derive(Debug, Clone, Default)]
pub struct BranchCutStats {
    /// Number of nodes explored.
    pub nodes_explored: usize,
    /// Number of cuts added.
    pub cuts_added: usize,
    /// LP solves performed.
    pub lp_solves: usize,
    /// Best integer solution found.
    pub best_integer_value: Option<f64>,
}

/// Result of branch-and-cut.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum BranchCutResult {
    /// Optimal solution found.
    Optimal,
    /// Feasible solution found (within gap).
    Feasible,
    /// Problem is infeasible.
    Infeasible,
    /// Node limit reached.
    NodeLimit,
}

/// Branch-and-cut solver for MIP.
#[derive(Debug)]
pub struct BranchCutSolver {
    /// Configuration.
    config: BranchCutConfig,
    /// Variable types.
    var_types: Vec<VarType>,
    /// Node queue.
    nodes: BinaryHeap<BranchNode>,
    /// Next node ID.
    next_node_id: usize,
    /// Best integer solution found.
    best_solution: Option<FxHashMap<VarId, BigRational>>,
    /// Best integer objective value.
    best_value: Option<BigRational>,
    /// Cutting plane generator.
    #[allow(dead_code)]
    cut_generator: CuttingPlaneGenerator,
    /// Statistics.
    stats: BranchCutStats,
}

impl BranchCutSolver {
    /// Create a new branch-and-cut solver.
    pub fn new(config: BranchCutConfig, var_types: Vec<VarType>) -> Self {
        // Collect integer variable IDs
        let integer_vars: crate::prelude::FxHashSet<VarId> = var_types
            .iter()
            .enumerate()
            .filter_map(|(i, &vt)| {
                if vt == VarType::Integer || vt == VarType::Binary {
                    Some(i)
                } else {
                    None
                }
            })
            .collect();

        Self {
            cut_generator: CuttingPlaneGenerator::new(integer_vars),
            config,
            var_types,
            nodes: BinaryHeap::new(),
            next_node_id: 0,
            best_solution: None,
            best_value: None,
            stats: BranchCutStats::default(),
        }
    }

    /// Solve the MIP problem.
    ///
    /// # Arguments
    ///
    /// * `root_solver` - Dual simplex solver at root node
    ///
    /// Returns the solution status.
    pub fn solve(&mut self, root_solver: &mut DualSimplex) -> BranchCutResult {
        // Solve root LP
        let root_result = root_solver.solve();
        self.stats.lp_solves += 1;

        match root_result {
            DualSimplexResult::Infeasible => return BranchCutResult::Infeasible,
            DualSimplexResult::Optimal => {
                // Check if root solution is integer-feasible
                let solution = root_solver.get_solution();
                if self.is_integer_feasible(&solution) {
                    self.best_solution = Some(solution.clone());
                    self.best_value = Some(root_solver.get_objective_value());
                    return BranchCutResult::Optimal;
                }

                // Create root node
                let root_node = BranchNode {
                    id: self.next_node_id,
                    lp_bound: root_solver.get_objective_value(),
                    bounds: FxHashMap::default(),
                    depth: 0,
                };
                self.next_node_id += 1;
                self.nodes.push(root_node);
            }
            _ => return BranchCutResult::NodeLimit,
        }

        // Branch-and-cut loop
        while let Some(node) = self.nodes.pop() {
            self.stats.nodes_explored += 1;

            if self.stats.nodes_explored >= self.config.max_nodes {
                return BranchCutResult::NodeLimit;
            }

            // Check if node can be pruned
            if self.can_prune(&node) {
                continue;
            }

            // Solve the LP relaxation at this node with updated variable bounds.
            let mut node_solver = root_solver.clone();

            // Inject tightened bounds as additional constraints:
            // lower bound: x_i >= lb  =>  -x_i <= -lb
            // upper bound: x_i <= ub  =>   x_i <=  ub
            for (&var_id, (lb, ub)) in &node.bounds {
                // Encode lb <= x_i: add constraint  x_i >= lb  as  -x_i <= -lb
                let n_vars = root_solver.num_vars();
                let mut lb_row = vec![BigRational::zero(); n_vars];
                if var_id < n_vars {
                    lb_row[var_id] = BigRational::from_integer((-1).into());
                }
                node_solver.add_constraint(lb_row, -lb.clone());

                // Encode x_i <= ub
                let mut ub_row = vec![BigRational::zero(); n_vars];
                if var_id < n_vars {
                    ub_row[var_id] = BigRational::from_integer(1.into());
                }
                node_solver.add_constraint(ub_row, ub.clone());
            }

            let lp_result = node_solver.solve();
            self.stats.lp_solves += 1;

            match lp_result {
                DualSimplexResult::Infeasible => {
                    // This branch is infeasible — prune.
                    continue;
                }
                DualSimplexResult::Optimal => {
                    let solution = node_solver.get_solution();
                    let obj_value = node_solver.get_objective_value();

                    // Prune if LP bound is already worse than the best known integer solution.
                    if let Some(ref best) = self.best_value
                        && obj_value >= *best
                    {
                        continue;
                    }

                    if self.is_integer_feasible(&solution) {
                        // Found an integer-feasible solution; update incumbent if better.
                        let improve = self.best_value.as_ref().is_none_or(|bv| obj_value < *bv);
                        if improve {
                            self.best_value = Some(obj_value.clone());
                            self.best_solution = Some(solution.clone());
                            self.stats.best_integer_value =
                                Some(obj_value.to_f64().unwrap_or(f64::INFINITY));
                        }
                        continue; // No further branching needed for this node.
                    }

                    // Select the most-infeasible fractional integer variable to branch on.
                    let Some((branch_var, branch_val)) = self.select_branch_variable(&solution)
                    else {
                        continue;
                    };

                    let floor_val = branch_val.floor();
                    let ceil_val = floor_val.clone() + BigRational::from_integer(1.into());

                    // Left child: x_i <= floor(val)
                    let mut left_bounds = node.bounds.clone();
                    {
                        let entry = left_bounds
                            .entry(branch_var)
                            .or_insert_with(|| (BigRational::zero(), floor_val.clone()));
                        entry.1 = entry.1.clone().min(floor_val);
                    }
                    let left_node = BranchNode {
                        id: self.next_node_id,
                        lp_bound: obj_value.clone(),
                        bounds: left_bounds,
                        depth: node.depth + 1,
                    };
                    self.next_node_id += 1;
                    self.nodes.push(left_node);

                    // Right child: x_i >= ceil(val)
                    let mut right_bounds = node.bounds.clone();
                    {
                        let entry = right_bounds
                            .entry(branch_var)
                            .or_insert_with(|| (ceil_val.clone(), BigRational::zero()));
                        entry.0 = entry.0.clone().max(ceil_val);
                    }
                    let right_node = BranchNode {
                        id: self.next_node_id,
                        lp_bound: obj_value,
                        bounds: right_bounds,
                        depth: node.depth + 1,
                    };
                    self.next_node_id += 1;
                    self.nodes.push(right_node);
                }
                _ => {
                    // Unbounded or unknown — skip this node.
                    continue;
                }
            }
        }

        if self.best_solution.is_some() {
            BranchCutResult::Feasible
        } else {
            BranchCutResult::Infeasible
        }
    }

    /// Check if a solution is integer-feasible.
    fn is_integer_feasible(&self, solution: &FxHashMap<VarId, BigRational>) -> bool {
        for (var, value) in solution {
            if (self.var_types.get(*var) == Some(&VarType::Integer)
                || self.var_types.get(*var) == Some(&VarType::Binary))
                && !self.is_integer_value(value)
            {
                return false;
            }
        }
        true
    }

    /// Check if a value is (approximately) integer.
    fn is_integer_value(&self, value: &BigRational) -> bool {
        let frac = value - value.floor();
        frac.is_zero() || frac < BigRational::new(1.into(), 1000000.into())
    }

    /// Check if a node can be pruned.
    fn can_prune(&self, node: &BranchNode) -> bool {
        if let Some(ref best_value) = self.best_value {
            // Prune if LP bound is worse than best integer solution
            node.lp_bound >= *best_value
        } else {
            false
        }
    }

    /// Select a variable to branch on.
    ///
    /// Returns (variable_id, fractional_value).
    fn select_branch_variable(
        &self,
        solution: &FxHashMap<VarId, BigRational>,
    ) -> Option<(VarId, BigRational)> {
        let mut candidates = Vec::new();

        for (var, value) in solution {
            if self.var_types.get(*var) == Some(&VarType::Integer) && !self.is_integer_value(value)
            {
                let frac = value - value.floor();
                candidates.push((*var, value.clone(), frac));
            }
        }

        if candidates.is_empty() {
            return None;
        }

        match self.config.branching_strategy {
            BranchingStrategy::MostFractional => {
                // Select variable with fractional part closest to 0.5
                candidates.sort_by(|a, b| {
                    let dist_a = (a.2.clone() - BigRational::new(1.into(), 2.into())).abs();
                    let dist_b = (b.2.clone() - BigRational::new(1.into(), 2.into())).abs();
                    dist_a.cmp(&dist_b)
                });
                Some((candidates[0].0, candidates[0].1.clone()))
            }
            BranchingStrategy::LeastFractional => {
                // Select variable with smallest fractional part
                candidates.sort_by(|a, b| a.2.cmp(&b.2));
                Some((candidates[0].0, candidates[0].1.clone()))
            }
            BranchingStrategy::LargestCoeff => {
                // For now, just pick first
                Some((candidates[0].0, candidates[0].1.clone()))
            }
        }
    }

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

    /// Get best solution found.
    pub fn best_solution(&self) -> Option<&FxHashMap<VarId, BigRational>> {
        self.best_solution.as_ref()
    }
}

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

    #[test]
    fn test_branch_cut_creation() {
        let config = BranchCutConfig::default();
        let var_types = vec![VarType::Integer, VarType::Integer];
        let solver = BranchCutSolver::new(config, var_types);

        assert_eq!(solver.stats().nodes_explored, 0);
    }

    #[test]
    fn test_is_integer_feasible() {
        let config = BranchCutConfig::default();
        let var_types = vec![VarType::Integer, VarType::Continuous];
        let solver = BranchCutSolver::new(config, var_types);

        let mut solution = FxHashMap::default();
        solution.insert(0, BigRational::from_integer(5.into()));
        solution.insert(1, BigRational::new(3.into(), 2.into())); // 1.5 (continuous, OK)

        assert!(solver.is_integer_feasible(&solution));

        // Make first variable fractional
        solution.insert(0, BigRational::new(5.into(), 2.into())); // 2.5 (not integer)
        assert!(!solver.is_integer_feasible(&solution));
    }
}