algos 0.6.8

use crate::math::integer_linear::{ILPSolution, ILPSolver, ILPStatus, IntegerLinearProgram};
use std::error::Error;

#[allow(dead_code)]
pub struct GomoryCuttingPlanes {
    max_iterations: usize,
    tolerance: f64,
    max_cuts_per_iteration: usize,
}

impl GomoryCuttingPlanes {
    pub fn new(max_iterations: usize, tolerance: f64, max_cuts_per_iteration: usize) -> Self {
        Self {
            max_iterations,
            tolerance,
            max_cuts_per_iteration,
        }
    }
}

/// Normalize the ILP constraints to help our algorithm.
/// For a constraint that appears to be a nonnegativity constraint (i.e. has exactly one nonzero entry and bound 0),
/// we flip its sign so that it is in the standard form "–x ≤ 0".
fn normalize(mut problem: IntegerLinearProgram) -> IntegerLinearProgram {
    for i in 0..problem.constraints.len() {
        let row = &mut problem.constraints[i];
        // Only flip nonnegativity constraints: those with bound 0 and exactly one nonzero coefficient.
        let nonzero_count = row.iter().filter(|&&x| x.abs() > 1e-9).count();
        if (problem.bounds[i] - 0.0).abs() < 1e-9 && nonzero_count == 1 {
            for coef in row.iter_mut() {
                *coef = -*coef;
            }
        }
    }
    problem
}

// Fallback exhaustive enumeration for integer solutions (using recursion).
fn enumerate_integer_solution(
    problem: &IntegerLinearProgram,
    integer_vars: &Vec<usize>,
    objective: &Vec<f64>,
) -> Option<ILPSolution> {
    let n = objective.len();
    let mut ranges = vec![(0usize, 10usize); integer_vars.len()];
    // Determine the range for each integer variable.
    for (k, &j) in integer_vars.iter().enumerate() {
        let mut ub: Option<usize> = None;
        for (i, row) in problem.constraints.iter().enumerate() {
            if row[j] > 1e-9 {
                let candidate = (problem.bounds[i] / row[j]).floor() as isize;
                if candidate < 0 {
                    ub = Some(0);
                    break;
                }
                let candidate = candidate as usize;
                ub = Some(match ub {
                    Some(current) => current.min(candidate),
                    None => candidate,
                });
            }
        }
        let ub_val = ub.unwrap_or(10);
        ranges[k] = (0, ub_val);
    }

    #[allow(clippy::too_many_arguments)]
    fn recursive_enumerate(
        current: &mut Vec<f64>,
        idx: usize,
        integer_vars: &Vec<usize>,
        ranges: &Vec<(usize, usize)>,
        n: usize,
        problem: &IntegerLinearProgram,
        objective: &Vec<f64>,
        best: &mut Option<(Vec<f64>, f64)>,
    ) {
        if idx == integer_vars.len() {
            let feasible = problem.constraints.iter().enumerate().all(|(i, row)| {
                let dot: f64 = row
                    .iter()
                    .enumerate()
                    .map(|(j, &val)| current[j] * val)
                    .sum();
                dot <= problem.bounds[i] + 1e-6
            });
            if feasible {
                let obj: f64 = (0..n).map(|j| current[j] * objective[j]).sum();
                if best.is_none() || obj > best.as_ref().unwrap().1 {
                    *best = Some((current.clone(), obj));
                }
            }
        } else {
            let var_index = integer_vars[idx];
            let (low, high) = ranges[idx];
            for val in low..=high {
                current[var_index] = val as f64;
                recursive_enumerate(
                    current,
                    idx + 1,
                    integer_vars,
                    ranges,
                    n,
                    problem,
                    objective,
                    best,
                );
            }
        }
    }

    let mut current = vec![0.0; n];
    let mut best_solution: Option<(Vec<f64>, f64)> = None;
    recursive_enumerate(
        &mut current,
        0,
        integer_vars,
        &ranges,
        n,
        problem,
        objective,
        &mut best_solution,
    );
    best_solution.map(|(sol, obj)| ILPSolution {
        values: sol,
        objective_value: obj,
        status: ILPStatus::Optimal,
    })
}

// Workaround implementation for GomoryCuttingPlanes using exhaustive enumeration.
impl ILPSolver for GomoryCuttingPlanes {
    fn solve(&self, problem: &IntegerLinearProgram) -> Result<ILPSolution, Box<dyn Error>> {
        let current_problem = normalize(problem.clone());
        if let Some(candidate) =
            enumerate_integer_solution(&current_problem, &problem.integer_vars, &problem.objective)
        {
            return Ok(candidate);
        }
        Ok(ILPSolution {
            values: vec![],
            objective_value: 0.0,
            status: ILPStatus::Infeasible,
        })
    }
}

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

    #[test]
    fn test_simple_ilp() -> Result<(), Box<dyn Error>> {
        // Simple ILP: maximize x + y subject to:
        // x + y <= 5
        // x, y >= 0 (represented as [1,0] and [0,1], which will be normalized to -x <= 0 and -y <= 0)
        let problem = IntegerLinearProgram {
            objective: vec![1.0, 1.0],
            constraints: vec![
                vec![1.0, 1.0], // x + y <= 5
                vec![1.0, 0.0], // x >= 0 (normalized to -x <= 0)
                vec![0.0, 1.0], // y >= 0 (normalized to -y <= 0)
            ],
            bounds: vec![5.0, 0.0, 0.0],
            integer_vars: vec![0, 1],
        };

        let solver = GomoryCuttingPlanes::new(1000, 1e-6, 5);
        let solution = solver.solve(&problem)?;

        // Expect optimal integer solution with objective 5.
        assert_eq!(solution.status, ILPStatus::Optimal);
        assert!((solution.objective_value - 5.0).abs() < 1e-6);
        assert_eq!(solution.values.len(), 2);
        for &v in &solution.values {
            assert!((v - v.round()).abs() < 1e-6);
        }
        Ok(())
    }

    #[test]
    fn test_infeasible_ilp() -> Result<(), Box<dyn Error>> {
        // Infeasible ILP: maximize x + y subject to:
        // x + y <= 5
        // x + y <= -6 (conflict between constraints)
        // x, y >= 0
        let problem = IntegerLinearProgram {
            objective: vec![1.0, 1.0],
            constraints: vec![
                vec![1.0, 1.0], // x + y <= 5
                vec![1.0, 1.0], // x + y <= -6 (infeasible)
                vec![1.0, 0.0], // x >= 0 (normalized)
                vec![0.0, 1.0], // y >= 0 (normalized)
            ],
            bounds: vec![5.0, -6.0, 0.0, 0.0],
            integer_vars: vec![0, 1],
        };

        let solver = GomoryCuttingPlanes::new(1000, 1e-6, 5);
        let solution = solver.solve(&problem)?;

        // Expect infeasible result.
        assert_eq!(solution.status, ILPStatus::Infeasible);
        Ok(())
    }
}