glpk-rust 0.2.1

use libc;

// Include the generated GLPK bindings
#[allow(non_upper_case_globals)]
#[allow(non_camel_case_types)]
#[allow(non_snake_case)]
#[allow(improper_ctypes)]
#[allow(dead_code)]
mod glpk {
    include!(concat!(env!("OUT_DIR"), "/bindings.rs"));
}

use libc::c_int;
use std::collections::HashMap;
use std::fmt;

pub type Bound = (i32, i32);
pub type ID<'a> = &'a str;
pub type Objective<'a> = HashMap<ID<'a>, f64>;

/// Thin wrapper around the numeric GLPK constants that glpk-sys
/// (or your bindgen output) does not re-export.
///
/// Feel free to extend this file if you need more flags later on.
pub mod glp_consts {
    // Row/column bound types
    pub const GLP_FR: i32 = 1;
    pub const GLP_LO: i32 = 2;
    pub const GLP_UP: i32 = 3;
    pub const GLP_DB: i32 = 4;
    pub const GLP_FX: i32 = 5;

    // Variable kinds
    pub const GLP_CV: i32 = 1;
    pub const GLP_IV: i32 = 2;
    pub const GLP_BV: i32 = 3;

    // Optimisation direction
    pub const GLP_MIN: i32 = 1;
    pub const GLP_MAX: i32 = 2;

    // Solution / MIP status
    pub const GLP_UNDEF: i32 = 1;
    pub const GLP_FEAS: i32 = 2;
    pub const GLP_INFEAS: i32 = 3;
    pub const GLP_NOFEAS: i32 = 4;
    pub const GLP_OPT: i32 = 5;
    pub const GLP_UNBND: i32 = 6;
}

#[derive(Debug, PartialEq)]
pub enum Status {
    Undefined = 1,
    Feasible = 2,
    Infeasible = 3,
    NoFeasible = 4,
    Optimal = 5,
    Unbounded = 6,
    SimplexFailed = 7,
    MIPFailed = 8,
    EmptySpace = 9,
}

#[derive(Clone, Hash)]
pub struct Variable<'a> {
    pub id: ID<'a>,
    pub bound: Bound,
}

#[derive(Clone, Hash)]
pub struct IntegerSparseMatrix {
    pub rows: Vec<i32>,
    pub cols: Vec<i32>,
    pub vals: Vec<i32>,
}

#[derive(Hash)]
pub struct SparseLEIntegerPolyhedron<'a> {
    pub a: IntegerSparseMatrix,
    pub b: Vec<Bound>,
    pub variables: Vec<Variable<'a>>,
    pub double_bound: bool,
}

impl<'a> SparseLEIntegerPolyhedron<'a> {
    pub fn get_hash(&self) -> u64 {
        use std::collections::hash_map::DefaultHasher;
        use std::hash::{Hash, Hasher};

        let mut hasher = DefaultHasher::new();
        self.hash(&mut hasher);
        let hash = hasher.finish();
        hash
    }
}

impl fmt::Display for SparseLEIntegerPolyhedron<'_> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let space = 4;
        let nr_rows = self.a.rows.iter().max().unwrap_or(&0) + 1;
        let nr_cols = self.a.cols.iter().max().unwrap_or(&0) + 1;
        // Create 2d array of zeroes
        let mut matrix = vec![vec![0; nr_cols as usize]; nr_rows as usize];
        for ((r, c), v) in self
            .a
            .rows
            .iter()
            .zip(self.a.cols.iter())
            .zip(self.a.vals.iter())
        {
            matrix[*r as usize][*c as usize] = *v;
        }
        for variable in &self.variables {
            write!(
                f,
                "{:>space$}",
                &variable.id.chars().take(3).collect::<String>()
            )?;
        }
        writeln!(f)?;
        for (idx, row) in matrix.iter().enumerate() {
            if self.double_bound {
                write!(f, "{:>space$} <= ", self.b[idx].0)?;
            }
            for val in row.iter() {
                write!(f, "{:>space$}", val)?;
            }
            write!(f, " <= {:>space$}", self.b[idx].1)?;
            writeln!(f)?;
        }
        Ok(())
    }
}

#[derive(Debug)]
pub struct Solution {
    pub status: Status,
    pub objective: f64,
    pub solution: HashMap<String, i32>,
    pub error: Option<String>,
}

pub struct SolutionResponse {
    pub solutions: Vec<Solution>,
}

use std::collections::BTreeMap;

#[derive(Debug)]
pub enum MatrixValidationError {
    LengthMismatch {
        rows: usize,
        cols: usize,
        vals: usize,
    },
    EmptyMatrix,
    NegativeRowIndex {
        index: usize,
        row: i32,
    },
    NegativeColIndex {
        index: usize,
        col: i32,
    },
    RowOutOfBounds {
        index: usize,
        row: i32,
        n_rows: usize,
    },
    ColOutOfBounds {
        index: usize,
        col: i32,
        n_cols: usize,
    },
}

impl fmt::Display for MatrixValidationError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            MatrixValidationError::LengthMismatch { rows, cols, vals } => {
                write!(f, "Matrix array length mismatch: rows={}, cols={}, vals={}", rows, cols, vals)
            }
            MatrixValidationError::EmptyMatrix => {
                write!(f, "Matrix is empty or all coefficients cancelled out")
            }
            MatrixValidationError::NegativeRowIndex { index, row } => {
                write!(f, "Negative row index at position {}: {}", index, row)
            }
            MatrixValidationError::NegativeColIndex { index, col } => {
                write!(f, "Negative column index at position {}: {}", index, col)
            }
            MatrixValidationError::RowOutOfBounds { index, row, n_rows } => {
                write!(f, "Row index {} at position {} exceeds bounds (n_rows={})", row, index, n_rows)
            }
            MatrixValidationError::ColOutOfBounds { index, col, n_cols } => {
                write!(f, "Column index {} at position {} exceeds bounds (n_cols={})", col, index, n_cols)
            }
        }
    }
}

impl std::error::Error for MatrixValidationError {}

#[derive(Debug)]
pub enum SolverError {
    /// Matrix validation errors (dimension mismatches, out of bounds, etc.)
    MatrixValidation(MatrixValidationError),
    /// The constraint matrix is empty
    EmptyConstraintMatrix,
    /// Number of variables doesn't match maximum column index
    VariableColumnMismatch {
        n_variables: usize,
        max_col_index: usize,
    },
    /// Number of bounds doesn't match maximum row index
    RowBoundMismatch {
        n_bounds: usize,
        max_row_index: usize,
    },
    /// Unknown solver status code
    UnknownStatus(i32),
}

impl fmt::Display for SolverError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            SolverError::MatrixValidation(e) => write!(f, "Matrix validation error: {}", e),
            SolverError::EmptyConstraintMatrix => {
                write!(f, "The constraint matrix A cannot be empty")
            }
            SolverError::VariableColumnMismatch { n_variables, max_col_index } => {
                write!(f, "The number of variables ({}) must be at least as large as the maximum column index ({})", n_variables, max_col_index)
            }
            SolverError::RowBoundMismatch { n_bounds, max_row_index } => {
                write!(f, "The number of bounds ({}) must be at least as large as the maximum row index ({})", n_bounds, max_row_index)
            }
            SolverError::UnknownStatus(code) => {
                write!(f, "Unknown solver status code: {}", code)
            }
        }
    }
}

impl std::error::Error for SolverError {}

impl From<MatrixValidationError> for SolverError {
    fn from(err: MatrixValidationError) -> Self {
        SolverError::MatrixValidation(err)
    }
}

pub fn validate_and_deduplicate_matrix(
    rows: &[i32],
    cols: &[i32],
    vals: &[i32],
    n_rows: usize,
    n_cols: usize,
) -> Result<(Vec<i32>, Vec<i32>, Vec<i32>), MatrixValidationError> {
    if rows.len() != cols.len() || rows.len() != vals.len() {
        return Err(MatrixValidationError::LengthMismatch {
            rows: rows.len(),
            cols: cols.len(),
            vals: vals.len(),
        });
    }

    if rows.is_empty() {
        return Err(MatrixValidationError::EmptyMatrix);
    }

    let mut deduped: BTreeMap<(i32, i32), i32> = BTreeMap::new();

    for (i, ((&row, &col), &val)) in rows.iter().zip(cols.iter()).zip(vals.iter()).enumerate() {
        if row < 0 {
            return Err(MatrixValidationError::NegativeRowIndex { index: i, row });
        }

        if col < 0 {
            return Err(MatrixValidationError::NegativeColIndex { index: i, col });
        }

        if row as usize >= n_rows {
            return Err(MatrixValidationError::RowOutOfBounds {
                index: i,
                row,
                n_rows,
            });
        }

        if col as usize >= n_cols {
            return Err(MatrixValidationError::ColOutOfBounds {
                index: i,
                col,
                n_cols,
            });
        }

        *deduped.entry((row, col)).or_insert(0) += val;
    }

    let mut clean_rows = Vec::with_capacity(deduped.len());
    let mut clean_cols = Vec::with_capacity(deduped.len());
    let mut clean_vals = Vec::with_capacity(deduped.len());

    for ((row, col), val) in deduped {
        // Remove coefficients that cancelled out, e.g. +3 and -3.
        if val != 0 {
            clean_rows.push(row);
            clean_cols.push(col);
            clean_vals.push(val);
        }
    }

    if clean_rows.is_empty() {
        return Err(MatrixValidationError::EmptyMatrix);
    }

    Ok((clean_rows, clean_cols, clean_vals))
}

// 1    GLP_FR      Free variable (−∞ < x < ∞)
// 2    GLP_LO      Variable with lower bound (x ≥ l)
// 3    GLP_UP      Variable with upper bound (x ≤ u)
// 4    GLP_DB      Double-bounded variable (l ≤ x ≤ u)
// 5    GLP_FX      Fixed variable (x = l = u)

// Value	Symbol	Meaning	Example
// 1	GLP_FR	Free (no bounds)	−∞<x<∞
// 2	GLP_LO	Lower bound only	x≥l
// 3	GLP_UP	Upper bound only	𝑥≤𝑢
// 4	GLP_DB	Double bounded (both bounds)	l≤x≤u
// 5	GLP_FX	Fixed (both bounds are equal)	x=l=u

// Variable kind
// 1 - continuous, 2 - integer, 3 - binary

// Solution status
// GLP_UNDEF   1  /* Solution is undefined */
// GLP_FEAS    2  /* Solution is feasible */
// GLP_INFEAS  3  /* No feasible solution exists */
// GLP_NOFEAS  4  /* No feasible solution exists (dual) */
// GLP_OPT     5  /* Optimal solution found */
// GLP_UNBND   6  /* Problem is unbounded */
/// Solves a set of integer linear programming (ILP) problems defined by the given polytope and objectives.
/// The function ignores objective variables not present in the polytope. All other variables not provided in objective but present in polytope will be set to 0.
///
/// # Arguments
///
/// * `polytope` - A reference to a `Polytope` struct that defines the constraints and variables of the ILP.
/// * `objectives` - A vector of `Objective` structs that define the objective functions to be optimized.
/// * `term_out` - A boolean flag indicating whether to enable terminal output for the GLPK solver.
///
/// # Returns
///
/// Returns `Ok(Vec<Solution>)` with a vector of `Solution` structs, each representing the result
/// of solving the ILP for one of the objectives.
///
/// Returns `Err(SolverError)` if there are validation errors with the input data.
///
/// # Errors
///
/// This function will return an error if:
/// * The constraint matrix is empty (`SolverError::EmptyConstraintMatrix`)
/// * The lengths of rows, columns, and values arrays don't match (`SolverError::MatrixValidation`)
/// * The number of variables doesn't match the maximum column index (`SolverError::VariableColumnMismatch`)
/// * The number of bounds doesn't match the maximum row index (`SolverError::RowBoundMismatch`)
/// * An unknown solver status is encountered (`SolverError::UnknownStatus`)
/// * Matrix validation fails due to invalid indices (`SolverError::MatrixValidation`)
///
/// # Examples
///
/// ```
/// use glpk_rust::*;
/// use std::collections::HashMap;
///
/// let variables = vec![
///     Variable { id: "x1", bound: (0, 5) },
///     Variable { id: "x2", bound: (0, 3) },
/// ];
///
/// let mut polytope = SparseLEIntegerPolyhedron {
///     a: IntegerSparseMatrix {
///         rows: vec![0, 1],
///         cols: vec![0, 1],
///         vals: vec![1, 1],
///     },
///     b: vec![(0, 10), (0, 8)],
///     variables,
///     double_bound: true,
/// };
///
/// let mut objective = HashMap::new();
/// objective.insert("x1", 1.0);
/// objective.insert("x2", 1.0);
/// let objectives = vec![objective];
///
/// let solutions = solve_ilps(&mut polytope, objectives, false, false, false)?;
/// for solution in solutions {
///     println!("{:?}", solution);
/// }
/// # Ok::<(), SolverError>(())
/// ```
pub fn solve_ilps<'a>(
    polytope: &mut SparseLEIntegerPolyhedron<'a>,
    objectives: Vec<Objective<'a>>,
    maximize: bool,
    presolve: bool,
    term_out: bool,
) -> Result<Vec<Solution>, SolverError> {
    // Initialize an empty vector to store solutions
    let mut solutions: Vec<Solution> = Vec::new();

    // Check that n_rows in shape is at least as large as the max row index in A
    let n_cols = polytope.variables.len();

    // If the polytope is empty, return an error
    if polytope.a.rows.is_empty() || polytope.a.cols.is_empty() {
        return Err(SolverError::EmptyConstraintMatrix);
    }

    // Check if rows, columns and values are the same length
    if (polytope.a.rows.len() != polytope.a.cols.len())
        || (polytope.a.rows.len() != polytope.a.vals.len())
        || (polytope.a.cols.len() != polytope.a.vals.len())
    {
        return Err(SolverError::MatrixValidation(
            MatrixValidationError::LengthMismatch {
                rows: polytope.a.rows.len(),
                cols: polytope.a.cols.len(),
                vals: polytope.a.vals.len(),
            },
        ));
    }

    // Get the maximum row and column indices from the constraint matrix A
    let poly_n_cols = (*polytope.a.cols.iter().max().unwrap() + 1) as usize;
    if n_cols < poly_n_cols {
        return Err(SolverError::VariableColumnMismatch {
            n_variables: n_cols,
            max_col_index: poly_n_cols,
        });
    }

    // Check that the number of rows is equal to the length of b
    let max_row_idx = (*polytope.a.rows.iter().max().unwrap() + 1) as usize;
    if max_row_idx > polytope.b.len() {
        return Err(SolverError::RowBoundMismatch {
            n_bounds: polytope.b.len(),
            max_row_index: max_row_idx,
        });
    }

    unsafe {
        // Enable or disable terminal output
        glpk::glp_term_out(term_out as c_int);

        // Create the problem
        let direction = if maximize { 2 } else { 1 };
        let lp = glpk::glp_create_prob();
        glpk::glp_set_obj_dir(lp, direction);

        // Add constraints (rows)
        glpk::glp_add_rows(lp, polytope.b.len() as i32);
        for (i, &b) in polytope.b.iter().enumerate() {
            let double_bound = if polytope.double_bound {
                glp_consts::GLP_DB
            } else {
                glp_consts::GLP_UP
            };
            glpk::glp_set_row_bnds(lp, (i + 1) as i32, double_bound, b.0 as f64, b.1 as f64);
        }

        // Add variables
        glpk::glp_add_cols(lp, polytope.variables.len() as i32);
        for (i, var) in polytope.variables.iter().enumerate() {
            // Set col bounds
            if var.bound.0 == var.bound.1 {
                glpk::glp_set_col_bnds(
                    lp,
                    (i + 1) as i32,
                    glp_consts::GLP_FX,
                    var.bound.0 as f64,
                    var.bound.1 as f64,
                );
            } else {
                glpk::glp_set_col_bnds(
                    lp,
                    (i + 1) as i32,
                    glp_consts::GLP_DB,
                    var.bound.0 as f64,
                    var.bound.1 as f64,
                );
            }

            // Set col kind - always integer variables
            // From GLPK docs: "Setting a column to GLP_BV has the same effect as if it were set to GLP_IV, its lower bound were
            // set 0, and its upper bound were set to 1."
            glpk::glp_set_col_kind(lp, (i + 1) as i32, glp_consts::GLP_IV);
        }

        let (clean_rows, clean_cols, clean_vals) = validate_and_deduplicate_matrix(
            &polytope.a.rows,
            &polytope.a.cols,
            &polytope.a.vals,
            polytope.b.len(),
            polytope.variables.len(),
        )?;
        
        // Set the constraint matrix
        let rows: Vec<i32> = std::iter::once(0)
            .chain(clean_rows.iter().map(|x| *x + 1))
            .collect();
        let cols: Vec<i32> = std::iter::once(0)
            .chain(clean_cols.iter().map(|x| *x + 1))
            .collect();
        let vals_f64: Vec<f64> = std::iter::once(0.0)
            .chain(clean_vals.iter().map(|x| *x as f64))
            .collect();

        // ne: The number of non-zero elements in the constraint matrix (not the number of rows or columns).
        // ia: An array of row indices (1-based) for each non-zero element.
        // ja: An array of column indices (1-based) for each non-zero element.
        // ar: An array of values corresponding to each non-zero element.
        glpk::glp_load_matrix(
            lp,
            (vals_f64.len() - 1) as i32,
            rows.as_ptr(),
            cols.as_ptr(),
            vals_f64.as_ptr(),
        );

        // Solve for multiple objectives
        for obj in objectives.iter() {
            // Setup empty solution
            let mut solution = Solution {
                status: Status::Undefined,
                objective: 0.0,
                solution: HashMap::new(),
                error: None,
            };

            // Update the objective function
            for (j, v) in polytope.variables.iter().enumerate() {
                let coef = obj.get(&v.id).unwrap_or(&0.0);
                glpk::glp_set_obj_coef(lp, (j + 1) as i32, *coef as f64);
            }

            // First, solve the LP relaxation to provide an initial basis for the MIP solver
            let mut simplex_params = glpk::glp_smcp::default();
            glpk::glp_init_smcp(&mut simplex_params);
            simplex_params.msg_lev = if term_out { 3 } else { 0 };

            let simplex_ret = glpk::glp_simplex(lp, &simplex_params);
            if simplex_ret != 0 {
                solution.status = Status::SimplexFailed;
                solution.error = Some(format!(
                    "GLPK simplex solver failed with code: {}",
                    simplex_ret
                ));
                solutions.push(solution);
                continue;
            }

            // Now solve the integer problem using the LP relaxation basis
            let mut mip_params = glpk::glp_iocp::default();
            glpk::glp_init_iocp(&mut mip_params);
            mip_params.presolve = if presolve { 1 } else { 0 };
            mip_params.msg_lev = if term_out { 3 } else { 0 };
            let mip_ret = glpk::glp_intopt(lp, &mip_params);

            if mip_ret != 0 {
                solution.status = Status::MIPFailed;
                solution.error = Some(format!("GLPK MIP solver failed with code: {}", mip_ret));
                solutions.push(solution);
                continue;
            }

            let status = glpk::glp_mip_status(lp);
            match status {
                1 => {
                    solution.status = Status::Undefined;
                    solution.error = Some("Solution is undefined".to_string());
                }
                2 => {
                    solution.status = Status::Feasible;
                    solution.objective = glpk::glp_mip_obj_val(lp);
                    for (j, var) in polytope.variables.iter().enumerate() {
                        let x = glpk::glp_mip_col_val(lp, (j + 1) as i32);
                        solution.solution.insert(var.id.to_string(), x as i32);
                    }
                }
                3 => {
                    solution.status = Status::Infeasible;
                    solution.error = Some("Infeasible solution exists".to_string());
                }
                4 => {
                    solution.status = Status::NoFeasible;
                    solution.error = Some("No feasible solution exists".to_string());
                }
                5 => {
                    solution.status = Status::Optimal;
                    solution.objective = glpk::glp_mip_obj_val(lp);
                    for (j, var) in polytope.variables.iter().enumerate() {
                        let x = glpk::glp_mip_col_val(lp, (j + 1) as i32);
                        solution.solution.insert(var.id.to_string(), x as i32);
                    }
                }
                6 => {
                    solution.status = Status::Unbounded;
                    solution.error = Some("Problem is unbounded".to_string());
                }
                x => {
                    // Clean up before returning error
                    glpk::glp_delete_prob(lp);
                    return Err(SolverError::UnknownStatus(x));
                }
            }
            solutions.push(solution);
        }

        // Clean up
        glpk::glp_delete_prob(lp);
        glpk::glp_free_env();

        Ok(solutions)
    }
}

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

    #[test]
    fn test_status_enum() {
        assert_eq!(Status::Optimal as i32, 5);
        assert_eq!(Status::Infeasible as i32, 3);
    }

    #[test]
    fn test_variable_creation() {
        let var = Variable {
            id: "x1",
            bound: (0, 10),
        };
        assert_eq!(var.id, "x1");
        assert_eq!(var.bound, (0, 10));
    }

    #[test]
    fn test_sparse_le_polyhedron_creation() {
        let variables = vec![
            Variable {
                id: "x1",
                bound: (0, 5),
            },
            Variable {
                id: "x2",
                bound: (0, 3),
            },
        ];

        let polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 1],
                cols: vec![0, 1],
                vals: vec![2, 3],
            },
            b: vec![(0, 10), (0, 15)],
            variables,
            double_bound: true,
        };

        assert_eq!(polytope.a.vals, vec![2, 3]);
        assert_eq!(polytope.b, vec![(0, 10), (0, 15)]);
        assert!(polytope.double_bound);
        assert_eq!(polytope.variables.len(), 2);
    }

    #[test]
    fn test_simple_conjunction() {
        let variables = vec![
            Variable {
                id: "x1",
                bound: (0, 1),
            },
            Variable {
                id: "x2",
                bound: (0, 1),
            },
            Variable {
                id: "x3",
                bound: (0, 1),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0, 0],
                cols: vec![0, 1, 2],
                vals: vec![1, 1, 1],
            },
            b: vec![(0, 3)],
            variables,
            double_bound: false,
        };
        let objectives = vec![HashMap::from_iter(vec![
            ("x1", 1.0),
            ("x2", 1.0),
            ("x3", 1.0),
        ])];
        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        let solution = &solutions[0];
        assert_eq!(solution.status, Status::Optimal);
        assert_eq!(solution.solution.get("x1"), Some(&1));
        assert_eq!(solution.solution.get("x2"), Some(&1));
        assert_eq!(solution.solution.get("x3"), Some(&1));
    }

    #[test]
    fn test_empty_matrix_handling() {
        let variables = vec![];
        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![],
                cols: vec![],
                vals: vec![],
            },
            b: vec![],
            variables,
            double_bound: false,
        };

        let objectives = vec![HashMap::new()];
        let result = solve_ilps(&mut polytope, objectives, false, false, false);

        assert!(result.is_err());
        assert!(matches!(result.unwrap_err(), SolverError::EmptyConstraintMatrix));
    }

    #[test]
    fn test_single_variable_optimization() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, false, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
        assert_eq!(solutions[0].solution.get("x"), Some(&0));
    }

    #[test]
    fn test_binary_variables() {
        let variables = vec![
            Variable {
                id: "b1",
                bound: (0, 1),
            },
            Variable {
                id: "b2",
                bound: (0, 1),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0],
                cols: vec![0, 1],
                vals: vec![1, 1],
            },
            b: vec![(1, 1)],
            variables,
            double_bound: false,
        };

        let mut objective = HashMap::new();
        objective.insert("b1", 1.0);
        objective.insert("b2", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
        let b1_val = solutions[0].solution.get("b1").unwrap();
        let b2_val = solutions[0].solution.get("b2").unwrap();
        assert_eq!(b1_val + b2_val, 1);
    }

    #[test]
    fn test_fixed_variables() {
        let variables = vec![
            Variable {
                id: "fixed",
                bound: (5, 5),
            },
            Variable {
                id: "free",
                bound: (0, 10),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0],
                cols: vec![0, 1],
                vals: vec![1, 1],
            },
            b: vec![(0, 15)],
            variables,
            double_bound: false,
        };

        let mut objective = HashMap::new();
        objective.insert("free", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
        assert_eq!(solutions[0].solution.get("fixed"), Some(&5));
        assert_eq!(solutions[0].solution.get("free"), Some(&10));
    }

    #[test]
    fn test_infeasible_problem() {
        let variables = vec![Variable {
            id: "x",
            bound: (5, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 3)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, false, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert!(matches!(solutions[0].status, Status::MIPFailed));
    }

    #[test]
    fn test_maximize_vs_minimize() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        let objectives = vec![objective.clone()];

        // Test minimization
        let min_solutions = solve_ilps(&mut polytope, objectives.clone(), false, false, false).unwrap();
        assert_eq!(min_solutions[0].solution.get("x"), Some(&0));

        // Test maximization
        let max_solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();
        assert_eq!(max_solutions[0].solution.get("x"), Some(&5));
    }

    #[test]
    fn test_multiple_objectives() {
        let variables = vec![
            Variable {
                id: "x",
                bound: (0, 10),
            },
            Variable {
                id: "y",
                bound: (0, 10),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0],
                cols: vec![0, 1],
                vals: vec![1, 1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let mut obj1 = HashMap::new();
        obj1.insert("x", 1.0);

        let mut obj2 = HashMap::new();
        obj2.insert("y", 1.0);

        let objectives = vec![obj1, obj2];

        let solutions = solve_ilps(&mut polytope, objectives, false, false, false).unwrap();

        assert_eq!(solutions.len(), 2);
        assert_eq!(solutions[0].status, Status::Optimal);
        assert_eq!(solutions[1].status, Status::Optimal);
    }

    #[test]
    fn test_objective_with_missing_variables() {
        let variables = vec![
            Variable {
                id: "x",
                bound: (0, 10),
            },
            Variable {
                id: "y",
                bound: (0, 10),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0],
                cols: vec![0, 1],
                vals: vec![1, 1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        objective.insert("z", 5.0); // Variable not in polytope
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, false, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
        assert_eq!(solutions[0].solution.get("y"), Some(&0));
    }

    #[test]
    fn test_negative_coefficients() {
        let variables = vec![
            Variable {
                id: "x",
                bound: (0, 10),
            },
            Variable {
                id: "y",
                bound: (0, 10),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0],
                cols: vec![0, 1],
                vals: vec![-1, 1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        objective.insert("y", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, false, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
    }

    #[test]
    fn test_zero_coefficient_matrix() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![0],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        let objectives = vec![objective];

        // Zero coefficients are removed during deduplication, resulting in empty matrix
        let result = solve_ilps(&mut polytope, objectives, false, false, false);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            SolverError::MatrixValidation(MatrixValidationError::EmptyMatrix)
        ));
    }

    #[test]
    fn test_large_bounds() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 1000000),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 500000)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
        assert_eq!(solutions[0].solution.get("x"), Some(&500000));
    }

    #[test]
    fn test_json_data_polyhedron() {
        let variables = vec![
            Variable {
                id: "a",
                bound: (0, 1),
            },
            Variable {
                id: "b",
                bound: (0, 1),
            },
            Variable {
                id: "18a7bec7bbb9fe127d6107f77af0f11b24a6a846",
                bound: (0, 1),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0, 0, 1, 1, 1, 2],
                cols: vec![0, 1, 2, 0, 1, 2, 2],
                vals: vec![-1, -1, 2, 1, 1, -2, -1],
            },
            b: vec![(0, 1), (0, 0), (0, -1)],
            variables,
            double_bound: false,
        };

        let mut objective = HashMap::new();
        objective.insert("a", 1.0);
        objective.insert("b", -1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        let solution = &solutions[0];
        assert_eq!(solution.status, Status::Optimal);
        assert_eq!(solution.objective, 1.0);
        assert_eq!(solution.solution.get("a"), Some(&1));
        assert_eq!(solution.solution.get("b"), Some(&0));
    }

    #[test]
    fn test_complex_constraint_system() {
        let variables = vec![
            Variable {
                id: "x1",
                bound: (0, 10),
            },
            Variable {
                id: "x2",
                bound: (0, 10),
            },
            Variable {
                id: "x3",
                bound: (0, 10),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0, 0, 1, 1, 1],
                cols: vec![0, 1, 2, 0, 1, 2],
                vals: vec![1, 2, 3, 2, 1, 1],
            },
            b: vec![(0, 15), (0, 10)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x1", 3.0);
        objective.insert("x2", 2.0);
        objective.insert("x3", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
    }

    #[test]
    fn test_double_bound_vs_single_bound() {
        let variables = vec![Variable {
            id: "x",
            bound: (2, 8),
        }];

        let mut polytope_double = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 10)],
            variables: variables.clone(),
            double_bound: true,
        };

        let mut polytope_single = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 10)],
            variables,
            double_bound: false,
        };

        let mut objective = HashMap::new();
        objective.insert("x", 1.0);
        let objectives = vec![objective];

        let solutions_double = solve_ilps(
            &mut polytope_double,
            objectives.clone(),
            false,
            false,
            false,
        ).unwrap();
        let solutions_single = solve_ilps(&mut polytope_single, objectives, false, false, false).unwrap();

        assert_eq!(solutions_double.len(), 1);
        assert_eq!(solutions_single.len(), 1);
        assert_eq!(solutions_double[0].status, Status::Optimal);
        assert_eq!(solutions_single[0].status, Status::Optimal);
    }

    #[test]
    fn test_sparse_matrix_with_gaps() {
        let variables = vec![
            Variable {
                id: "x1",
                bound: (0, 10),
            },
            Variable {
                id: "x2",
                bound: (0, 10),
            },
            Variable {
                id: "x3",
                bound: (0, 10),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 2],
                cols: vec![0, 2],
                vals: vec![1, 1],
            },
            b: vec![(0, 5), (0, 5), (0, 5)],
            variables,
            double_bound: true,
        };

        let mut objective = HashMap::new();
        objective.insert("x1", 1.0);
        objective.insert("x2", 1.0);
        objective.insert("x3", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, false, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
    }

    #[test]
    fn test_mismatched_matrix_dimensions() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 1],
                cols: vec![0],
                vals: vec![1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let objectives = vec![HashMap::new()];
        let result = solve_ilps(&mut polytope, objectives, false, false, false);

        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            SolverError::MatrixValidation(MatrixValidationError::LengthMismatch { .. })
        ));
    }

    #[test]
    fn test_variable_column_mismatch() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0],
                cols: vec![1],
                vals: vec![1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let objectives = vec![HashMap::new()];
        let result = solve_ilps(&mut polytope, objectives, false, false, false);

        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            SolverError::VariableColumnMismatch { .. }
        ));
    }

    #[test]
    fn test_row_bound_mismatch() {
        let variables = vec![Variable {
            id: "x",
            bound: (0, 10),
        }];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 1],
                cols: vec![0, 0],
                vals: vec![1, 1],
            },
            b: vec![(0, 5)],
            variables,
            double_bound: true,
        };

        let objectives = vec![HashMap::new()];
        let result = solve_ilps(&mut polytope, objectives, false, false, false);

        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            SolverError::RowBoundMismatch { .. }
        ));
    }

    #[test]
    fn test_fixed_boolean_variables() {
        // Test that boolean variables can be fixed to specific values
        // Using bounds that won't be classified as binary to avoid GLPK binary variable behavior
        let variables = vec![
            Variable {
                id: "fixed_zero",
                bound: (0, 0),
            },
            Variable {
                id: "fixed_one",
                bound: (1, 1),
            },
            Variable {
                id: "free_binary",
                bound: (0, 1),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0, 0],
                cols: vec![0, 1, 2],
                vals: vec![1, 1, 1],
            },
            b: vec![(0, 10)],
            variables,
            double_bound: false,
        };

        let mut objective = HashMap::new();
        objective.insert("fixed_zero", 1.0);
        objective.insert("fixed_one", 1.0);
        objective.insert("free_binary", 1.0);
        let objectives = vec![objective];

        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);

        // Fixed to 0 should remain 0 even when classified as binary
        assert_eq!(solutions[0].solution.get("fixed_zero"), Some(&0));

        // Fixed to 1 should remain 1 even when classified as binary
        assert_eq!(solutions[0].solution.get("fixed_one"), Some(&1));

        // Free binary variable should be able to take value 0 or 1
        // With maximization, it should be 1
        assert_eq!(solutions[0].solution.get("free_binary"), Some(&1));
    }

    #[test]
    fn test_atmost_one_of_constraints_with_bounds_fixed_in_variables() {
        // 0 <= [  1,   1,   1,   3] <=   4
        // 0 <= [ -1,  -1,  -1,  -5] <=  -2
        let variables = vec![
            Variable {
                id: "a",
                bound: (0, 1),
            },
            Variable {
                id: "b",
                bound: (0, 1),
            },
            Variable {
                id: "c",
                bound: (0, 1),
            },
            Variable {
                id: "A",
                bound: (1, 1),
            }, // Fixed to 1 means at most one of a,b,c can be 1
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0, 0, 0, 1, 1, 1, 1],
                cols: vec![3, 0, 1, 2, 3, 0, 1, 2],
                vals: vec![3, 1, 1, 1, -5, -1, -1, -1],
            },
            b: vec![(0, 4), (0, -2)],
            variables,
            double_bound: false,
        };

        // Try to take a b and c
        let objective = HashMap::from([("b", 1.0), ("a", 1.0), ("c", 1.0)]);
        let objectives = vec![objective];
        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();
        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);

        // Check that at most one of a,b,c is 1
        let sum_abc = solutions[0].solution.get("a").unwrap()
            + solutions[0].solution.get("b").unwrap()
            + solutions[0].solution.get("c").unwrap();
        assert!(sum_abc <= 1);
    }

    #[test]
    fn test_that_shape_too_small() {
        let variables = vec![
            Variable {
                id: "a",
                bound: (0, 1),
            },
            Variable {
                id: "b",
                bound: (0, 1),
            },
            Variable {
                id: "c",
                bound: (0, 1),
            },
            Variable {
                id: "R",
                bound: (0, 1),
            },
        ];

        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0],
                cols: vec![0, 1],
                vals: vec![0, 0],
            },
            b: vec![(0, 1)],
            variables,
            double_bound: false,
        };
        let objective = HashMap::from([("a", 1.0), ("b", 1.0), ("c", 1.0), ("R", 1.0)]);
        let objectives = vec![objective];

        // Zero coefficients are removed during deduplication, resulting in empty matrix
        let result = solve_ilps(&mut polytope, objectives, true, false, false);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            SolverError::MatrixValidation(MatrixValidationError::EmptyMatrix)
        ));
    }

    #[test]
    fn test_duplicate_matrix_indices() {
        let variables = vec![
            Variable {
                id: "x1",
                bound: (0, 10),
            },
            Variable {
                id: "x2",
                bound: (0, 10),
            },
        ];

        // Create a matrix with duplicate (row, col) indices
        // These should be automatically merged: 1 + 2 = 3 at position (0, 0)
        let mut polytope = SparseLEIntegerPolyhedron {
            a: IntegerSparseMatrix {
                rows: vec![0, 0], // Same row
                cols: vec![0, 0], // Same column - will be deduplicated!
                vals: vec![1, 2], // Will be summed to 3
            },
            b: vec![(0, 10)],
            variables,
            double_bound: false,
        };

        let objective = HashMap::from([("x1", 1.0)]);
        let objectives = vec![objective];
        let solutions = solve_ilps(&mut polytope, objectives, true, false, false).unwrap();

        // Should succeed with deduplication
        assert_eq!(solutions.len(), 1);
        assert_eq!(solutions[0].status, Status::Optimal);
    }
}