lpsolve 1.0.1

//! Wrapper for lpsolve
//!
//! ![Build status](https://gitlab.com/cmr/rust-lpsolve/badges/master/build.svg)
//! [![Crates.io](https://img.shields.io/crates/v/lpsolve.svg)](https://crates.io/crates/lpsolve)
//!
//! lpsolve is a free software (LGPL) solver for mixed integer linear programming problems. The
//! documentation here is nonexistent when it comes to understanding how to model systems or
//! precisely what the consequences of various methods are.  The [upstream
//! documentation](http://lpsolve.sourceforge.net/5.5/) for lpsolve is much more comprehensive.
//!
//! # Quick Start
//!
//! ```ignore
//! use lpsolve::prelude::*;  // Import everything you need
//!
//! // Terse builder API - maximize 30*x1 + 50*x2 subject to constraints
//! let solution = Problem::builder()
//!     .cols(2)
//!     .max(&[30.0, 50.0])           // Shorthand for maximize + objective
//!     .le(&[2.0, 3.0], 100.0)       // 2*x1 + 3*x2 <= 100
//!     .le(&[1.0, 2.0], 60.0)        // x1 + 2*x2 <= 60
//!     .non_negative()               // All variables >= 0
//!     .solve()?;
//!
//! println!("Status: {:?}", solution.status());
//! println!("Objective: {}", solution.objective_value());
//! println!("Variables: {:?}", solution.variables());
//! ```
//!
//! See the [prelude] module for convenient imports and shorthand methods.
//!
//! # Indexing Convention
//!
//! **IMPORTANT**: This library follows lpsolve's 1-based indexing convention:
//!
//! - **Columns (variables)**: Indexed from **1** to `num_cols()` inclusive
//! - **Rows (constraints)**: Indexed from **1** to `num_rows()` inclusive
//!   - Row **0** is reserved for the objective function in some operations
//!
//! When passing coefficient arrays (like in `set_objective_function` or `add_constraint`),
//! the first element (index 0) represents the constant term, and subsequent elements
//! correspond to variables 1, 2, 3, etc.
//!
//! ```ignore
//! // Setting objective: 10*x1 + 20*x2 + 5 (constant)
//! problem.set_objective_function(&[5.0, 10.0, 20.0])?;
//! //                                 ^    ^     ^
//! //                            constant x1    x2
//! ```
//!
//! # Error Handling
//!
//! All getters and setters perform bounds checking and return descriptive errors.
//!
//! The performance of lpsolve is mediocre compared to commercial solvers and some other free
//! software solvers. Performance here is how how long it takes to solve [benchmark
//! models](http://plato.asu.edu/bench.html).
//!
//! If you need help choosing a solver, the following is an excellent report:
//!
//! <http://prod.sandia.gov/techlib/access-control.cgi/2013/138847.pdf>
//!
//! # Stability
//!
//! `lpsolve-sys` is versioned separately from this wrapper. This wrapper is provisionally
//! unstable, but the functions that are currently wrapped are not likely to change.
//!
//! # License
//!
//! This crate and `lpsolve-sys` are licensed under either of
//!
//!  * Apache License, Version 2.0, ([LICENSE-APACHE](LICENSE-APACHE) or http://www.apache.org/licenses/LICENSE-2.0)
//!  * MIT license ([LICENSE-MIT](LICENSE-MIT) or http://opensource.org/licenses/MIT)
//!
//! at your option. However, please note that lpsolve itself is LGPL. The default configuration right
//! now builds a bundled copy of lpsolve and links to it statically.

extern crate libc;
extern crate lpsolve_sys as lp;
#[macro_use]
extern crate bitflags;

pub mod error;
pub mod builder;
pub mod prelude;

pub use error::{LpSolveError, Result};
pub use builder::{ProblemBuilder, Solution, SolverConfig};

use std::ffi::{CStr, CString};
use std::io::Write;
use std::ops::Deref;
use error::{CReturnCode, OpContext, EntityType};

#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Ord, PartialOrd, Hash)]
pub enum Verbosity {
    Critical = 1,
    Severe = 2,
    Important = 3,
    Normal = 4,
    Detailed = 5,
    Full = 6,
}

#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Ord, PartialOrd, Hash)]
pub enum ConstraintType {
    Le = 1,
    Eq = 3,
    Ge = 2,
    Free = 0,
}

#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Ord, PartialOrd, Hash)]
pub enum SOSType {
    Type1 = 1,
    Type2 = 2,
}

#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Ord, PartialOrd, Hash)]
pub enum VarType {
    Binary = 1,
    Float = 0,
}

#[repr(u8)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Ord, PartialOrd, Hash)]
pub enum BoundsMode {
    Restrictive = 1,
    None = 0,
}

#[repr(C)]
#[derive(Debug, PartialEq, Eq, Ord, PartialOrd, Hash, Clone, Copy)]
pub enum SimplexType {
    PrimalPrimal = 5,
    DualPrimal = 6,
    PrimalDual = 9,
    DualDual = 10,
}

#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, Ord, PartialOrd, Hash)]
pub enum SolveStatus {
    OutOfMemory = -2,
    NotRun = -1,
    Optimal = 0,
    Suboptimal = 1,
    Infeasible = 2,
    Unbounded = 3,
    Degenerate = 4,
    NumericalFailure = 5,
    UserAbort = 6,
    Timeout = 7,
    Presolved = 8,
    ProcFail = 9,
    ProcBreak = 11,
    FeasibleFound = 12,
    NoFeasibleFound = 13,
}

bitflags! {
    #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
    pub struct MPSOptions: ::libc::c_int {
        const CRITICAL = 1;
        const SEVERE = 2;
        const IMPORTANT = 3;
        const NORMAL = 4;
        const DETAILED = 5;
        const FULL = 6;
        const FREE = 8;
        const IBM = 16;
        const NEGOBJCONST = 32;
    }
}

bitflags! {
    /// Presolve methods for problem reduction
    #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
    pub struct PresolveMode: ::libc::c_int {
        /// No presolve
        const NONE = 0;
        /// Presolve rows
        const ROWS = 1;
        /// Presolve columns
        const COLS = 2;
        /// Eliminate linearly dependent rows
        const LINDEP = 4;
        /// Aggregate if possible (NOT IMPLEMENTED)
        const AGGREGATE = 8;
        /// Create sparser problem (NOT IMPLEMENTED)
        const SPARSER = 16;
        /// Convert constraints to SOS1 where possible
        const SOS = 32;
        /// Reduce MIP problem
        const REDUCEMIP = 64;
        /// Simplify knapsack-type constraints
        const KNAPSACK = 128;
        /// Eliminate constraints by substitution
        const ELIMEQ2 = 256;
        /// Aggressive bound tightening on implied free variables
        const IMPLIEDFREE = 512;
        /// Reduce GCD
        const REDUCEGCD = 1024;
        /// Fix variables based on probing
        const PROBEFIX = 2048;
        /// Probe to reduce problem
        const PROBEREDUCE = 4096;
        /// Identify and delete redundant rows
        const ROWDOMINATE = 8192;
        /// Identify and delete redundant columns
        const COLDOMINATE = 16384;
        /// Merge rows
        const MERGEROWS = 32768;
        /// Convert implied slack variables to equalities
        const IMPLIEDSLK = 65536;
        /// Fix columns based on dual values
        const COLFIXDUAL = 131072;
        /// Tighten bounds on variables
        const BOUNDS = 262144;
        /// Calculate duals
        const DUALS = 524288;
        /// Calculate sensitivity duals
        const SENSDUALS = 1048576;

        /// Typical presolve (rows + cols)
        const TYPICAL = Self::ROWS.bits() | Self::COLS.bits();
    }
}

bitflags! {
    /// Scaling modes for numerical stability
    /// Note: The base scaling method (0-7) and modifiers (8+) are combined
    #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
    pub struct ScalingMode: ::libc::c_int {
        /// No scaling
        const NONE = 0;
        /// Scale to convergence using largest absolute value
        const EXTREME = 1;
        /// Scale based on range of values
        const RANGE = 2;
        /// Scale by mean value
        const MEAN = 3;
        /// Scale by geometric mean
        const GEOMETRIC = 4;
        /// Future use 1
        const FUTURE1 = 5;
        /// Future use 2
        const FUTURE2 = 6;
        /// Curtis-Reid scaling
        const CURTISREID = 7;

        // Scaling modifiers (can be OR'd with base methods)
        /// Use quadratic scaling
        const QUADRATIC = 8;
        /// Use logarithmic scaling
        const LOGARITHMIC = 16;
        /// Scale to power of 2
        const POWER2 = 32;
        /// Ensure no scaled number exceeds 1
        const EQUILIBRATE = 64;
        /// Apply to integer columns
        const INTEGERS = 128;
        /// Update incrementally every solve
        const DYNUPDATE = 256;
        /// Only scale rows
        const ROWSONLY = 512;
        /// Only scale columns
        const COLSONLY = 1024;

        // Useful combinations
        /// Typical scaling (usually GEOMETRIC + EQUILIBRATE + INTEGERS)
        const TYPICAL = Self::GEOMETRIC.bits() | Self::EQUILIBRATE.bits() | Self::INTEGERS.bits();
        /// All methods
        const ALL = Self::EXTREME.bits() | Self::RANGE.bits() | Self::MEAN.bits() | Self::GEOMETRIC.bits();
    }
}

/// Pivoting rules for simplex algorithm
///
/// The pivot rule determines which variable enters the basis at each iteration.
/// Different rules offer tradeoffs between iteration count and computation per iteration.
#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum PivotRule {
    /// Select the first column that improves the objective (Dantzig's rule)
    ///
    /// **Fastest per iteration** but may require many more iterations.
    /// Use for small problems or when you need quick approximate solutions.
    FirstIndex = 0,

    /// Steepest edge pricing (default)
    ///
    /// **Best overall performance** for most problems. Balances iteration count
    /// and per-iteration cost. Uses geometric information to select pivots that
    /// give the steepest descent/ascent in objective space.
    SteepestEdge = 1,

    /// Devex pricing (Devex algorithm)
    ///
    /// **Good middle ground** between FirstIndex and SteepestEdge. Cheaper per
    /// iteration than steepest edge but typically requires fewer iterations than
    /// FirstIndex. Good for medium-sized problems.
    Devex = 2,

    /// Random pivot selection
    ///
    /// **Experimental use only**. Rarely useful except for testing or specific
    /// degenerate problems where it can help break cycling.
    Random = 3,
}

/// Branch and bound rules for MIP (Mixed Integer Programming)
///
/// When branching on a fractional variable, determines whether to try
/// rounding up or down first. Can significantly impact solve time.
#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum BranchRule {
    /// Branch up first (ceiling) - default
    ///
    /// For a fractional value like 2.7, creates child nodes with x ≥ 3 first.
    /// Often better when the optimal solution tends toward upper bounds.
    Ceiling = 0,

    /// Branch down first (floor)
    ///
    /// For a fractional value like 2.7, creates child nodes with x ≤ 2 first.
    /// Can be better when the optimal solution tends toward lower bounds.
    Floor = 1,

    /// Automatic selection based on problem heuristics
    ///
    /// **Recommended for most problems**. Uses pseudocost information and other
    /// heuristics to decide branching direction dynamically. Typically outperforms
    /// fixed strategies on non-trivial MIP problems.
    Automatic = 2,
}

bitflags! {
    /// Anti-degeneracy strategies
    #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
    pub struct AntiDegen: ::libc::c_int {
        /// No anti-degeneracy
        const NONE = 0;
        /// Fixed variable perturbation
        const FIXEDVARS = 1;
        /// Column check
        const COLUMNCHECK = 2;
        /// Stalling detection
        const STALLING = 4;
        /// Numerical failure detection
        const NUMFAILURE = 8;
        /// Lost feasibility detection
        const LOSTFEAS = 16;
        /// Infeasibility detection
        const INFEASIBLE = 32;
        /// Dynamic strategy
        const DYNAMIC = 64;
        /// Apply during branch-and-bound
        const DURINGBB = 128;
        /// RHS perturbation
        const RHSPERTURB = 256;
        /// Bound flipping
        const BOUNDFLIP = 512;

        /// Default anti-degeneracy (FIXEDVARS | STALLING)
        const DEFAULT = Self::FIXEDVARS.bits() | Self::STALLING.bits();
    }
}

bitflags! {
    /// Improvement strategies
    #[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
    pub struct ImproveMode: ::libc::c_int {
        /// No iterative improvement
        const NONE = 0;
        /// Improve solution
        const SOLUTION = 1;
        /// Improve dual feasibility
        const DUALFEAS = 2;
        /// Improve theta gap
        const THETAGAP = 4;
        /// Use branch-and-bound simplex
        const BBSIMPLEX = 8;

        /// Default improvement (DUALFEAS | THETAGAP)
        const DEFAULT = Self::DUALFEAS.bits() | Self::THETAGAP.bits();
        /// Inverse improvement (SOLUTION | THETAGAP)
        const INVERSE = Self::SOLUTION.bits() | Self::THETAGAP.bits();
    }
}

/// Basis crash modes for initial basis construction
///
/// Determines how to construct an initial basis before starting the simplex algorithm.
/// A good initial basis can significantly reduce the number of iterations needed.
#[repr(C)]
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum BasisCrash {
    /// No crash - use all-slack basis
    ///
    /// Starts with all slack variables in the basis (identity matrix).
    /// **Safest** but may require more iterations. Use when crash heuristics fail
    /// or for very small problems where crash overhead isn't worth it.
    None = 0,

    /// Most feasible basis
    ///
    /// Attempts to find a basis that is most likely to be feasible.
    /// **Good general purpose choice**. Reduces Phase 1 iterations at the cost
    /// of some setup time. Recommended for most problems.
    MostFeasible = 1,

    /// Least degenerate basis
    ///
    /// Constructs a basis to minimize degeneracy (multiple optimal solutions at a vertex).
    /// Can help when the problem is known to have significant degeneracy, but adds
    /// overhead. Typically only useful for specialized problems.
    LeastDegen = 2,
}

/// A linear programming problem.
pub struct Problem {
    lprec: *mut lp::lprec,
}

macro_rules! cptr {
    ($e:expr) => {
        if $e.is_null() {
            None
        } else {
            Some(Problem { lprec: $e })
        }
    };
}

unsafe extern "C" fn write_modeldata(
    val: *mut libc::c_void,
    buf: *mut libc::c_char,
) -> libc::c_int {
    use std::panic::AssertUnwindSafe;
    use std::panic::catch_unwind;

    // Catch panics to prevent unwinding through C code
    let result = catch_unwind(AssertUnwindSafe(|| unsafe {
        // SAFETY: val is a pointer to a pointer to a trait object (&mut dyn Write)
        // We dereference once to get the fat pointer, then dereference again to get the trait object
        let writer_ptr = val as *mut *mut dyn Write;
        let writer = &mut **writer_ptr;
        let buf = CStr::from_ptr(buf);
        writer.write_all(buf.to_bytes())
    }));

    match result {
        Ok(Ok(())) => 0, // No panic, write succeeded
        Ok(Err(_)) => 1, // No panic, write failed
        Err(_) => 1,     // Panic occurred, treat as write failure
    }
}

impl Problem {
    /// Create a new problem builder
    pub fn builder() -> ProblemBuilder {
        ProblemBuilder::new()
    }

    /// Helper to check buffer size and bounds in one go
    fn check_buffer_and_bounds(
        &self,
        buffer_len: usize,
        required: usize,
        index: libc::c_int,
        min: libc::c_int,
        max: libc::c_int,
        context: OpContext,
    ) -> Result<()> {
        if buffer_len < required {
            return Err(LpSolveError::BufferTooSmall {
                required,
                provided: buffer_len,
                context,
            });
        }
        if index < min || index > max {
            return Err(LpSolveError::IndexOutOfBounds {
                index,
                min,
                max,
                context,
            });
        }
        Ok(())
    }

    /// Initialize an empty problem with space for `rows` and `cols`.
    ///
    /// # Note
    /// When using `add_constraint` to build the problem dynamically, some versions
    /// of the underlying lpsolve library may have issues when `rows` is 0. Consider
    /// using [`Problem::builder()`] instead, which handles initialization correctly.
    pub fn new(rows: libc::c_int, cols: libc::c_int) -> Option<Problem> {
        let ptr = unsafe { lp::make_lp(rows, cols) };
        cptr!(ptr)
    }

    /// Reads an lp-format model from `path`.
    pub fn read_lp<P: Deref<Target = CStr>, C: Deref<Target = CStr>>(
        path: &P,
        verbosity: Verbosity,
        initial_name: &C,
    ) -> Option<Problem> {
        let ptr = unsafe {
            lp::read_LP(
                path.as_ptr() as *mut _,
                verbosity as libc::c_int,
                initial_name.as_ptr() as *mut _,
            )
        };
        cptr!(ptr)
    }

    /// Read an mps-format model from `path` using the "free" formatting.
    pub fn read_freemps<P: Deref<Target = CStr>>(path: &P, options: MPSOptions) -> Option<Problem> {
        let ptr = unsafe { lp::read_freeMPS(path.as_ptr() as *mut _, options.bits()) };
        cptr!(ptr)
    }

    /// Read an mps-format model from `path` using the fixed formatting.
    pub fn read_fixedmps<P: Deref<Target = CStr>>(
        path: &P,
        options: MPSOptions,
    ) -> Option<Problem> {
        let ptr = unsafe { lp::read_MPS(path.as_ptr() as *mut _, options.bits()) };
        cptr!(ptr)
    }

    /// Write an lp-format model into `out`.
    pub fn write_lp(&self, out: &mut dyn Write) -> Result<()> {
        // Create a raw pointer to the trait object
        let out_ptr = out as *mut dyn Write;
        let out_ptr_ptr = &out_ptr as *const *mut dyn Write;
        unsafe {
            lp::write_lpex(
                self.lprec,
                out_ptr_ptr as *mut libc::c_void,
                write_modeldata,
            )
        }.check_with(LpSolveError::WriteError)
    }

    /// Write an mps-format model into `out` using the fixed formatting.
    pub fn write_fixedmps(&self, out: &mut dyn Write) -> Result<()> {
        self.write_mps(out, 1)
    }

    /// Write an mps-format model into `out` using the "free" formatting.
    pub fn write_freemps(&self, out: &mut dyn Write) -> Result<()> {
        self.write_mps(out, 2)
    }

    /// Write an mps-format model into `out` using `formatting`.
    ///
    /// `formatting` must be 1 for fixed or 2 for free.
    pub fn write_mps(&self, out: &mut dyn Write, formatting: libc::c_int) -> Result<()> {
        debug_assert!(formatting == 1 || formatting == 2);
        // Create a raw pointer to the trait object
        let out_ptr = out as *mut dyn Write;
        let out_ptr_ptr = &out_ptr as *const *mut dyn Write;
        unsafe {
            lp::MPS_writefileex(
                self.lprec,
                formatting,
                out_ptr_ptr as *mut libc::c_void,
                write_modeldata,
            )
        }.check_with(LpSolveError::WriteError)
    }

    /// Reserve enough memory for `rows` and `cols`.
    ///
    /// If `rows` or `cols` are less than the current number of rows or columns, the additional
    /// rows and columns will be deleted.
    pub fn resize(&mut self, rows: libc::c_int, cols: libc::c_int) -> Result<()> {
        unsafe { lp::resize_lp(self.lprec, rows, cols) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: None,
                context: OpContext::Resize,
            })
    }

    /// Add a column to the model.
    ///
    /// If `values` is empty, an all-zero column will be added.
    ///
    /// The slice must have exactly `num_rows() + 1` elements (including objective coefficient).
    pub fn add_column(&mut self, values: &mut [f64]) -> Result<()> {
        let expected = (self.num_rows() as usize).saturating_add(1);
        if values.len() != expected {
            return Err(LpSolveError::DimensionMismatch {
                expected,
                actual: values.len(),
                context: OpContext::AddColumn,
            });
        }
        unsafe { lp::add_column(self.lprec, values.as_mut_ptr()) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(self.num_cols() + 1),
                context: OpContext::AddColumn,
            })
    }

    /// Add a column to the model, scattering `values` by `indices`.
    ///
    /// The values for the column are taken from `values`. The value from `values[i]` will be
    /// placed into row `indices[i]`. Both slices must have the same length.
    pub fn add_column_scatter(&mut self, values: &mut [f64], indices: &mut [libc::c_int]) -> Result<()> {
        if values.len() != indices.len() {
            return Err(LpSolveError::DimensionMismatch {
                expected: values.len(),
                actual: indices.len(),
                context: OpContext::AddColumn,
            });
        }
        unsafe {
            lp::add_columnex(
                self.lprec,
                values.len() as libc::c_int,
                values.as_mut_ptr(),
                indices.as_mut_ptr(),
            )
        }.check_with(LpSolveError::ColumnOperationFailed {
            column: Some(self.num_cols() + 1),
            context: OpContext::AddColumn,
        })
    }

    /// Read a column from the model.
    ///
    /// The slice must have space for at least `num_rows() + 1` elements.
    pub fn get_column(&self, values: &mut [f64], column: libc::c_int) -> Result<()> {
        self.check_buffer_and_bounds(
            values.len(),
            (self.num_rows() as usize).saturating_add(1),
            column,
            1,
            self.num_cols(),
            OpContext::GetColumn,
        )?;
        unsafe { lp::get_column(self.lprec, column, values.as_mut_ptr()) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row: 0,
                col: column,
                op: error::MatrixOp::GetColumn,
            })
    }

    /// Read a row from the model.
    ///
    /// The slice must have space for at least `num_cols() + 1` elements.
    pub fn get_row(&self, values: &mut [f64], row: libc::c_int) -> Result<()> {
        self.check_buffer_and_bounds(
            values.len(),
            (self.num_cols() as usize).saturating_add(1),
            row,
            1,
            self.num_rows(),
            OpContext::GetRow,
        )?;
        unsafe { lp::get_row(self.lprec, row, values.as_mut_ptr()) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row,
                col: 0,
                op: error::MatrixOp::GetRow,
            })
    }

    /// Set (replace) an entire row in the model.
    ///
    /// The slice must have at least `num_cols() + 1` elements.
    pub fn set_row(&mut self, row: libc::c_int, values: &mut [f64]) -> Result<()> {
        self.check_buffer_and_bounds(
            values.len(),
            (self.num_cols() as usize).saturating_add(1),
            row,
            1,
            self.num_rows(),
            OpContext::SetRow,
        )?;
        unsafe { lp::set_row(self.lprec, row, values.as_mut_ptr()) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row,
                col: 0,
                op: error::MatrixOp::SetRow,
            })
    }

    /// Set (replace) an entire row using sparse format.
    ///
    /// Both slices must have the same length.
    pub fn set_rowex(&mut self, row: libc::c_int, values: &mut [f64], indices: &mut [libc::c_int]) -> Result<()> {
        if row < 1 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 1,
                max: self.num_rows(),
                context: OpContext::SetRow,
            });
        }
        if values.len() != indices.len() {
            return Err(LpSolveError::DimensionMismatch {
                expected: values.len(),
                actual: indices.len(),
                context: OpContext::SetRow,
            });
        }
        unsafe {
            lp::set_rowex(
                self.lprec,
                row,
                values.len() as libc::c_int,
                values.as_mut_ptr(),
                indices.as_mut_ptr(),
            )
        }.check_with(LpSolveError::MatrixOperationFailed {
            row,
            col: 0,
            op: error::MatrixOp::SetRow,
        })
    }

    /// Set (replace) an entire column in the model.
    ///
    /// The slice must have at least `num_rows() + 1` elements.
    pub fn set_column(&mut self, col: libc::c_int, values: &mut [f64]) -> Result<()> {
        self.check_buffer_and_bounds(
            values.len(),
            (self.num_rows() as usize).saturating_add(1),
            col,
            1,
            self.num_cols(),
            OpContext::SetColumn,
        )?;
        unsafe { lp::set_column(self.lprec, col, values.as_mut_ptr()) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row: 0,
                col,
                op: error::MatrixOp::SetColumn,
            })
    }

    /// Set (replace) an entire column using sparse format.
    ///
    /// Both slices must have the same length.
    pub fn set_columnex(
        &mut self,
        col: libc::c_int,
        values: &mut [f64],
        indices: &mut [libc::c_int],
    ) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetColumn,
            });
        }
        if values.len() != indices.len() {
            return Err(LpSolveError::DimensionMismatch {
                expected: values.len(),
                actual: indices.len(),
                context: OpContext::SetColumn,
            });
        }
        unsafe {
            lp::set_columnex(
                self.lprec,
                col,
                values.len() as libc::c_int,
                values.as_mut_ptr(),
                indices.as_mut_ptr(),
            )
        }.check_with(LpSolveError::MatrixOperationFailed {
            row: 0,
            col,
            op: error::MatrixOp::SetColumn,
        })
    }

    /// Add a constraint to the model.
    ///
    /// The constraint is that `coeffs * vars OP target`, where `OP` is specified by `kind`.
    ///
    /// For optimal performance, use the `matrix_builder` method and add the objective function
    /// first. This method is otherwise very slow for large models.
    ///
    /// The slice must have exactly `num_cols() + 1` elements (including constant term).
    pub fn add_constraint(&mut self, coeffs: &mut [f64], target: f64, kind: ConstraintType) -> Result<()> {
        let expected = (self.num_cols() as usize).saturating_add(1);
        if coeffs.len() != expected {
            return Err(LpSolveError::DimensionMismatch {
                expected,
                actual: coeffs.len(),
                context: OpContext::AddConstraint,
            });
        }
        unsafe {
            lp::add_constraint(
                self.lprec,
                coeffs.as_mut_ptr(),
                kind as libc::c_int,
                target,
            )
        }.check_with(LpSolveError::ConstraintAdditionFailed { row: None })
    }

    /// Add a [Special Ordered Set](http://lpsolve.sourceforge.net/5.5/SOS.htm) constraint.
    ///
    /// The `weights` are scattered by `variables`, that is, `weights[i]` will be specified for
    /// column `variables[i]`. Both slices must have the same length.
    pub fn add_sos_constraint(
        &mut self,
        name: &CStr,
        sostype: SOSType,
        priority: libc::c_int,
        weights: &mut [f64],
        variables: &mut [libc::c_int],
    ) -> Result<()> {
        if weights.len() != variables.len() {
            return Err(LpSolveError::DimensionMismatch {
                expected: weights.len(),
                actual: variables.len(),
                context: OpContext::AddConstraint,
            });
        }
        let result = unsafe {
            lp::add_SOS(
                self.lprec,
                name.as_ptr() as *mut _,
                sostype as libc::c_int,
                priority,
                weights.len() as libc::c_int,
                variables.as_mut_ptr(),
                weights.as_mut_ptr(),
            )
        };
        if result == 0 {
            Err(LpSolveError::SOSConstraintError)
        } else {
            Ok(())
        }
    }

    /// Delete a column from the model.
    ///
    /// The other columns are shifted leftward. `col` cannot be 0, as that column represents the
    /// RHS, which must always be present.
    pub fn del_column(&mut self, col: libc::c_int) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::DeleteColumn,
            });
        }
        unsafe { lp::del_column(self.lprec, col) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::DeleteColumn,
            })
    }

    /// Delete a constraint from the model.
    ///
    /// The other constraints are shifted upwards. `row` cannot be 0, as that row represents the
    /// objective function, which must always be present.
    pub fn del_constraint(&mut self, row: libc::c_int) -> Result<()> {
        if row < 1 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 1,
                max: self.num_rows(),
                context: OpContext::DeleteRow,
            });
        }
        unsafe { lp::del_constraint(self.lprec, row) }
            .check_with(LpSolveError::ConstraintAdditionFailed {
                row: Some(row),
            })
    }

    /// Returns `true` if the specified column can be negative, `false` otherwise.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [0, num_cols()].
    pub fn is_negative(&self, col: libc::c_int) -> Result<bool> {
        if col < 0 || col > self.num_cols() {
            Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 0,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            })
        } else {
            Ok(unsafe { lp::is_negative(self.lprec, col) } == 1)
        }
    }

    /// Set a variable to be either binary or floating point.
    pub fn set_variable_type(&mut self, col: libc::c_int, vartype: VarType) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetVariableType,
            });
        }
        unsafe { lp::set_binary(self.lprec, col, vartype as libc::c_uchar) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetVariableType,
            })
    }

    /// Set a variable to be binary (convenience method)
    pub fn set_binary(&mut self, col: libc::c_int, is_binary: bool) -> Result<()> {
        self.set_variable_type(col, if is_binary { VarType::Binary } else { VarType::Float })
    }

    /// Get the type of a variable.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [0, num_cols()].
    pub fn get_variable_type(&self, col: libc::c_int) -> Result<VarType> {
        if col < 0 || col > self.num_cols() {
            Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 0,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            })
        } else {
            let res = if unsafe { lp::is_binary(self.lprec, col) } == 1 {
                VarType::Binary
            } else {
                VarType::Float
            };
            Ok(res)
        }
    }

    /// Set the upper and lower bounds of a variable.
    pub fn set_bounds(&mut self, col: libc::c_int, lower: f64, upper: f64) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetBounds,
            });
        }
        unsafe { lp::set_bounds(self.lprec, col, lower, upper) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetBounds,
            })
    }

    /// Set the bounds mode to 'tighten'.
    ///
    /// If the bounds mode is `true`, then when `set_bounds`, `set_lower_bound`, or
    /// `set_upper_bound` is used and the provided bound is less restrictive than the current bound
    /// (ie, its absolute value is larger), then the request will be ignored. However, if the new
    /// bound is more restrictive (ie, its absolute value is smaller) the request will go through.
    ///
    /// If the bounds mode is `false`, the bounds will always be set as specified.
    pub fn set_bounds_mode(&mut self, tighten: BoundsMode) {
        unsafe { lp::set_bounds_tighter(self.lprec, tighten as libc::c_uchar) };
    }

    /// Get the bounds mode.
    ///
    /// See `set_bounds_mode` for what this value means.
    pub fn get_bounds_mode(&self) -> BoundsMode {
        if unsafe { lp::get_bounds_tighter(self.lprec) } == 1 {
            BoundsMode::Restrictive
        } else {
            BoundsMode::None
        }
    }

    /// Set the type of a constraint.
    pub fn set_constraint_type(&mut self, row: libc::c_int, contype: ConstraintType) -> Result<()> {
        if row < 1 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 1,
                max: self.num_rows(),
                context: OpContext::SetRow,
            });
        }
        unsafe { lp::set_constr_type(self.lprec, row, contype as libc::c_int) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row,
                col: 0,
                op: error::MatrixOp::SetRow,
            })
    }

    /// Get the type of a constraint.
    ///
    /// # Errors
    /// Returns an error if `row` is out of bounds [0, num_rows()] or if the constraint type is invalid.
    pub fn get_constraint_type(&self, row: libc::c_int) -> Result<ConstraintType> {
        if row < 0 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 0,
                max: self.num_rows(),
                context: OpContext::GetRow,
            });
        }
        let res = unsafe { lp::get_constr_type(self.lprec, row) };
        match res {
            1 => Ok(ConstraintType::Le),
            2 => Ok(ConstraintType::Ge),
            3 => Ok(ConstraintType::Eq),
            _ => Err(LpSolveError::MatrixOperationFailed {
                row,
                col: 0,
                op: error::MatrixOp::GetRow,
            }),
        }
    }

    /// Set a variable to be unbounded.
    pub fn set_unbounded(&mut self, col: libc::c_int) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetUnbounded,
            });
        }
        unsafe { lp::set_unbounded(self.lprec, col) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetUnbounded,
            })
    }

    /// Check if a variable is unbounded.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [0, num_cols()].
    pub fn is_unbounded(&self, col: libc::c_int) -> Result<bool> {
        if col < 0 || col > self.num_cols() {
            Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 0,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            })
        } else {
            Ok(unsafe { lp::is_unbounded(self.lprec, col) } == 1)
        }
    }

    /// Set the practical value for "infinite"
    ///
    /// This is the bound of the absolute value of all variables. If the absolute value of a
    /// variable is larger than this, it is considered to have diverged.
    pub fn set_infinite(&mut self, infinity: f64) {
        unsafe { lp::set_infinite(self.lprec, infinity) };
    }

    /// Get the value of "infinite"
    ///
    /// See set_infinite for more details.
    pub fn get_infinite(&self) -> f64 {
        unsafe { lp::get_infinite(self.lprec) }
    }

    /// Set a variable's integer type.
    pub fn set_integer(&mut self, col: libc::c_int, must_be_integer: bool) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetInteger,
            });
        }
        unsafe { lp::set_int(self.lprec, col, if must_be_integer { 1 } else { 0 }) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetInteger,
            })
    }

    /// Check if a variable is an integer.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [0, num_cols()].
    pub fn is_integer(&self, col: libc::c_int) -> Result<bool> {
        if col < 0 || col > self.num_cols() {
            Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 0,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            })
        } else {
            Ok(unsafe { lp::is_int(self.lprec, col) } == 1)
        }
    }

    /// Set a variable to be semi-continuous.
    ///
    /// A semi-continuous variable can be either 0 or within its bounds, but nothing in between 0 and the lower bound.
    pub fn set_semicont(&mut self, col: libc::c_int, must_be_semicont: bool) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetSemicont,
            });
        }
        unsafe { lp::set_semicont(self.lprec, col, if must_be_semicont { 1 } else { 0 }) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetSemicont,
            })
    }

    /// Check if a variable is semi-continuous.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [0, num_cols()].
    pub fn is_semicont(&self, col: libc::c_int) -> Result<bool> {
        if col < 0 || col > self.num_cols() {
            Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 0,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            })
        } else {
            Ok(unsafe { lp::is_semicont(self.lprec, col) } == 1)
        }
    }

    /// Set variable weights for branching decisions in MIP problems.
    ///
    /// The weights slice must have at least `num_cols() + 1` elements (index 0 is ignored).
    pub fn set_var_weights(&mut self, weights: &mut [f64]) -> Result<()> {
        let required = (self.num_cols() as usize).saturating_add(1);
        if weights.len() < required {
            return Err(LpSolveError::BufferTooSmall {
                required,
                provided: weights.len(),
                context: OpContext::SetColumn,
            });
        }
        unsafe { lp::set_var_weights(self.lprec, weights.as_mut_ptr()) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: None,
                context: OpContext::SetColumn,
            })
    }

    /// Get the branching priority for a variable.
    pub fn get_var_priority(&self, col: libc::c_int) -> libc::c_int {
        unsafe { lp::get_var_priority(self.lprec, col) }
    }

    /// Check if a variable is a member of a Special Ordered Set.
    pub fn is_sos_var(&self, col: libc::c_int) -> bool {
        unsafe { lp::is_SOS_var(self.lprec, col) == 1 }
    }

    /// Sets the objective function.
    ///
    /// The slice must have exactly `num_cols() + 1` elements (including constant term at index 0).
    pub fn set_objective_function(&mut self, coeffs: &mut [f64]) -> Result<()> {
        let expected = (self.num_cols() as usize).saturating_add(1);
        if coeffs.len() != expected {
            return Err(LpSolveError::DimensionMismatch {
                expected,
                actual: coeffs.len(),
                context: OpContext::SetObjective,
            });
        }
        unsafe { lp::set_obj_fn(self.lprec, coeffs.as_mut_ptr()) }
            .check_with(LpSolveError::ObjectiveFunctionError)
    }

    /// Scatters `coeffs` into the objective function coefficients with `indices`.
    ///
    /// Both slices must have the same length.
    pub fn scatter_objective_function(&mut self, coeffs: &mut [f64], indices: &mut [libc::c_int]) -> Result<()> {
        if coeffs.len() != indices.len() {
            return Err(LpSolveError::DimensionMismatch {
                expected: coeffs.len(),
                actual: indices.len(),
                context: OpContext::SetObjective,
            });
        }
        unsafe {
            lp::set_obj_fnex(
                self.lprec,
                coeffs.len() as libc::c_int,
                coeffs.as_mut_ptr(),
                indices.as_mut_ptr(),
            )
        }.check_with(LpSolveError::ObjectiveFunctionError)
    }

    /// Sets the range of a constraint.
    ///
    /// If the constraint type is `<=`, then the "actual" constraint will be `RHS - range <= v <= RHS`.
    ///
    /// If the constraint type is `>=`, then the "actual" constraint will be `RHS <= v <= RHS + range`.
    ///
    /// This puts a bound on the constraint and can give the solver more freedom, and is more
    /// efficient than adding an extra constraint.
    pub fn set_constraint_range(&mut self, row: libc::c_int, range: f64) -> Result<()> {
        if row < 1 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 1,
                max: self.num_rows(),
                context: OpContext::SetRow,
            });
        }
        unsafe { lp::set_rh_range(self.lprec, row, range) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row,
                col: 0,
                op: error::MatrixOp::SetRow,
            })
    }

    /// Get the range on a constraint if one is set.
    ///
    /// Returns `Ok(None)` if no range is set (default infinite range).
    ///
    /// # Errors
    /// Returns an error if `row` is out of bounds [0, num_rows()].
    pub fn get_constraint_range(&self, row: libc::c_int) -> Result<Option<f64>> {
        if row < 0 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 0,
                max: self.num_rows(),
                context: OpContext::GetRow,
            });
        }
        let delta = unsafe { lp::get_rh_range(self.lprec, row) };
        if delta == self.get_infinite() {
            Ok(None)
        } else {
            Ok(Some(delta))
        }
    }

    /// Set the right-hand side value of a constraint.
    pub fn set_rh(&mut self, row: libc::c_int, value: f64) -> Result<()> {
        if row < 1 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 1,
                max: self.num_rows(),
                context: OpContext::SetRow,
            });
        }
        unsafe { lp::set_rh(self.lprec, row, value) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row,
                col: 0,
                op: error::MatrixOp::SetRow,
            })
    }

    /// Get the right-hand side value of a constraint.
    ///
    /// # Errors
    /// Returns an error if `row` is out of bounds [1, num_rows()].
    pub fn get_rh(&self, row: libc::c_int) -> Result<f64> {
        if row < 1 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 1,
                max: self.num_rows(),
                context: OpContext::GetRow,
            });
        }
        Ok(unsafe { lp::get_rh(self.lprec, row) })
    }

    /// Set the timeout for the solver in seconds.
    ///
    /// A timeout of 0 means no timeout (unlimited time).
    /// When the timeout is reached, solve() will return with SolveStatus::Timeout.
    pub fn set_timeout(&mut self, seconds: libc::c_long) {
        unsafe { lp::set_timeout(self.lprec, seconds) };
    }

    /// Get the current timeout setting in seconds.
    ///
    /// Returns 0 if no timeout is set.
    pub fn get_timeout(&self) -> libc::c_long {
        unsafe { lp::get_timeout(self.lprec) }
    }

    /// Set the optimization sense to maximize the objective function.
    pub fn set_maximize(&mut self) {
        unsafe { lp::set_maxim(self.lprec) };
    }

    /// Set the optimization sense to minimize the objective function.
    pub fn set_minimize(&mut self) {
        unsafe { lp::set_minim(self.lprec) };
    }

    /// Check if the problem is set to maximize (true) or minimize (false).
    pub fn is_maximize(&self) -> bool {
        unsafe { lp::is_maxim(self.lprec) == 1 }
    }

    /// Check if the problem is set to minimize (true) or maximize (false).
    pub fn is_minimize(&self) -> bool {
        !self.is_maximize()
    }

    /// Solve the model.
    pub fn solve(&mut self) -> SolveStatus {
        use SolveStatus::*;
        match unsafe { lp::solve(self.lprec) } {
            -2 => OutOfMemory,
            -1 => NotRun,
            0 => Optimal,
            1 => Suboptimal,
            2 => Infeasible,
            3 => Unbounded,
            4 => Degenerate,
            5 => NumericalFailure,
            6 => UserAbort,
            7 => Timeout,
            9 => Presolved,
            10 => ProcFail,
            11 => ProcBreak,
            12 => FeasibleFound,
            13 => NoFeasibleFound,
            status => panic!("unknown solve status {}", status),
        }
    }

    /// Read out the values assigned to variables from the most recent `solve`.
    ///
    /// Returns `None` if `vars` does not have at least as many elements as the underlying model
    /// has columns. Otherwise, returns `Some` with the slice truncated to the number of columns.
    pub fn get_solution_variables<'a>(&self, vars: &'a mut [f64]) -> Option<&'a mut [f64]> {
        let cols = self.num_cols();
        if vars.len() < cols as usize {
            None
        } else {
            unsafe { lp::get_variables(self.lprec, vars.as_mut_ptr()) };
            Some(&mut vars[..cols as usize])
        }
    }

    /// Get the objective function value from the most recent `solve`.
    pub fn get_objective(&self) -> f64 {
        unsafe { lp::get_objective(self.lprec) }
    }

    /// Get the total number of iterations performed during the last solve.
    pub fn get_total_iter(&self) -> libc::c_long {
        unsafe { lp::get_total_iter(self.lprec) }
    }

    /// Get the total number of branch-and-bound nodes evaluated during the last solve.
    pub fn get_total_nodes(&self) -> libc::c_long {
        unsafe { lp::get_total_nodes(self.lprec) }
    }

    /// Get the time elapsed during the last solve, in seconds.
    pub fn time_elapsed(&self) -> f64 {
        unsafe { lp::time_elapsed(self.lprec) }
    }

    /// Get the dual solution (shadow prices) from the most recent `solve`.
    ///
    /// The slice must have space for at least `num_rows() + num_cols()` elements.
    /// Returns `None` if the slice is too small, otherwise returns the dual values.
    /// Elements [0..num_rows()] correspond to constraint shadow prices.
    /// Elements [num_rows()..num_rows()+num_cols()] correspond to reduced costs.
    pub fn get_dual_solution<'a>(&self, duals: &'a mut [f64]) -> Option<&'a mut [f64]> {
        // Avoid integer overflow when adding dimensions
        let rows = self.num_rows() as usize;
        let cols = self.num_cols() as usize;
        let total = rows.saturating_add(cols);
        if duals.len() < total {
            None
        } else {
            unsafe { lp::get_dual_solution(self.lprec, duals.as_mut_ptr()) };
            Some(&mut duals[..total])
        }
    }

    /// Get lambda multipliers (Lagrangian multipliers) from the most recent `solve`.
    ///
    /// The slice must have space for at least `num_rows() + 1` elements.
    pub fn get_lambda<'a>(&self, lambda: &'a mut [f64]) -> Option<&'a mut [f64]> {
        let rows = (self.num_rows() as usize).saturating_add(1);
        if lambda.len() < rows {
            None
        } else {
            unsafe { lp::get_lambda(self.lprec, lambda.as_mut_ptr()) };
            Some(&mut lambda[..rows])
        }
    }

    /// Get the constraint values from the most recent `solve`.
    ///
    /// The slice must have space for at least `num_rows() + 1` elements.
    pub fn get_constraints<'a>(&self, constraints: &'a mut [f64]) -> Option<&'a mut [f64]> {
        let rows = (self.num_rows() as usize).saturating_add(1);
        if constraints.len() < rows {
            None
        } else {
            unsafe { lp::get_constraints(self.lprec, constraints.as_mut_ptr()) };
            Some(&mut constraints[..rows])
        }
    }

    /// Get sensitivity analysis for right-hand side values.
    ///
    /// Returns sensitivity ranges for constraint RHS values.
    /// Each of `duals`, `dualsfrom`, and `dualstill` must have space for at least `num_rows() + 1` elements.
    pub fn get_sensitivity_rhs(
        &self,
        duals: &mut [f64],
        dualsfrom: &mut [f64],
        dualstill: &mut [f64],
    ) -> Result<()> {
        let required = (self.num_rows() as usize).saturating_add(1);
        let context = OpContext::GetSensitivityRhs;

        for (buffer_len, _name) in [(duals.len(), "duals"), (dualsfrom.len(), "dualsfrom"), (dualstill.len(), "dualstill")] {
            if buffer_len < required {
                return Err(LpSolveError::BufferTooSmall {
                    required,
                    provided: buffer_len,
                    context,
                });
            }
        }

        unsafe {
            lp::get_sensitivity_rhs(
                self.lprec,
                duals.as_mut_ptr(),
                dualsfrom.as_mut_ptr(),
                dualstill.as_mut_ptr(),
            )
        }.check_with(LpSolveError::MatrixOperationFailed {
            row: 0,
            col: 0,
            op: error::MatrixOp::Get,
        })
    }

    /// Get sensitivity analysis for objective function coefficients.
    ///
    /// Returns sensitivity ranges for objective coefficients.
    /// Each of `objfrom` and `objtill` must have space for at least `num_cols() + 1` elements.
    pub fn get_sensitivity_obj(&self, objfrom: &mut [f64], objtill: &mut [f64]) -> Result<()> {
        let required = (self.num_cols() as usize).saturating_add(1);
        let context = OpContext::GetSensitivityObj;

        for (buffer_len, _name) in [(objfrom.len(), "objfrom"), (objtill.len(), "objtill")] {
            if buffer_len < required {
                return Err(LpSolveError::BufferTooSmall {
                    required,
                    provided: buffer_len,
                    context,
                });
            }
        }

        unsafe {
            lp::get_sensitivity_obj(self.lprec, objfrom.as_mut_ptr(), objtill.as_mut_ptr())
        }.check_with(LpSolveError::MatrixOperationFailed {
            row: 0,
            col: 0,
            op: error::MatrixOp::Get,
        })
    }

    /// Set the presolve mode using type-safe flags.
    ///
    /// Presolve reduces the problem size before solving.
    pub fn set_presolve(&mut self, mode: PresolveMode) {
        unsafe { lp::set_presolve(self.lprec, mode.bits(), -1) };
    }

    /// Get the current presolve settings.
    pub fn get_presolve(&self) -> PresolveMode {
        let bits = unsafe { lp::get_presolve(self.lprec) };
        PresolveMode::from_bits_truncate(bits)
    }

    /// Set the scaling mode.
    ///
    /// Scaling improves numerical stability.
    pub fn set_scaling(&mut self, mode: ScalingMode) {
        unsafe { lp::set_scaling(self.lprec, mode.bits()) };
    }

    /// Get the current scaling mode.
    pub fn get_scaling(&self) -> ScalingMode {
        let bits = unsafe { lp::get_scaling(self.lprec) };
        ScalingMode::from_bits_truncate(bits)
    }

    /// Set the simplex algorithm type.
    pub fn set_simplex_type(&mut self, simplex_type: SimplexType) {
        unsafe { lp::set_simplextype(self.lprec, simplex_type as libc::c_int) };
    }

    /// Get the current simplex algorithm type.
    pub fn get_simplex_type(&self) -> libc::c_int {
        unsafe { lp::get_simplextype(self.lprec) }
    }

    /// Set whether to prefer the dual simplex method.
    pub fn set_prefer_dual(&mut self, prefer_dual: bool) {
        unsafe { lp::set_preferdual(self.lprec, if prefer_dual { 1 } else { 0 }) };
    }

    /// Set the pivoting strategy.
    ///
    /// Affects how the simplex algorithm chooses pivot elements.
    pub fn set_pivoting(&mut self, rule: PivotRule) {
        unsafe { lp::set_pivoting(self.lprec, rule as libc::c_int) };
    }

    /// Get the current pivoting strategy.
    pub fn get_pivoting(&self) -> PivotRule {
        let value = unsafe { lp::get_pivoting(self.lprec) };
        match value {
            0 => PivotRule::FirstIndex,
            1 => PivotRule::SteepestEdge,
            2 => PivotRule::Devex,
            3 => PivotRule::Random,
            _ => PivotRule::FirstIndex, // Default fallback
        }
    }

    /// Set the integer tolerance (epsilon for integer).
    ///
    /// Variables within this tolerance of an integer are considered integer.
    pub fn set_epsint(&mut self, epsint: f64) {
        unsafe { lp::set_epsint(self.lprec, epsint) };
    }

    /// Get the integer tolerance.
    pub fn get_epsint(&self) -> f64 {
        unsafe { lp::get_epsint(self.lprec) }
    }

    /// Set the feasibility tolerance (epsilon for bounds).
    pub fn set_epsb(&mut self, epsb: f64) {
        unsafe { lp::set_epsb(self.lprec, epsb) };
    }

    /// Get the feasibility tolerance.
    pub fn get_epsb(&self) -> f64 {
        unsafe { lp::get_epsb(self.lprec) }
    }

    /// Set the optimality tolerance (epsilon for dual).
    pub fn set_epsd(&mut self, epsd: f64) {
        unsafe { lp::set_epsd(self.lprec, epsd) };
    }

    /// Get the optimality tolerance.
    pub fn get_epsd(&self) -> f64 {
        unsafe { lp::get_epsd(self.lprec) }
    }

    /// Set the pivot zero tolerance (epsilon for pivot).
    pub fn set_epspivot(&mut self, epspivot: f64) {
        unsafe { lp::set_epspivot(self.lprec, epspivot) };
    }

    /// Get the pivot zero tolerance.
    pub fn get_epspivot(&self) -> f64 {
        unsafe { lp::get_epspivot(self.lprec) }
    }

    /// Set the MIP gap tolerance.
    ///
    /// The solver stops when the gap between the best integer solution and the
    /// best possible solution is less than this value (absolute or relative).
    pub fn set_mip_gap(&mut self, absolute: bool, gap: f64) {
        unsafe { lp::set_mip_gap(self.lprec, if absolute { 1 } else { 0 }, gap) };
    }

    /// Get the MIP gap for the given mode.
    ///
    /// If `absolute` is true, returns the absolute MIP gap, otherwise returns the relative gap.
    pub fn get_mip_gap(&self, absolute: bool) -> f64 {
        unsafe { lp::get_mip_gap(self.lprec, if absolute { 1 } else { 0 }) }
    }

    /// Set the branch-and-bound rule for MIP problems.
    pub fn set_bb_rule(&mut self, rule: BranchRule) {
        unsafe { lp::set_bb_rule(self.lprec, rule as libc::c_int) };
    }

    /// Get the branch-and-bound rule.
    pub fn get_bb_rule(&self) -> BranchRule {
        let value = unsafe { lp::get_bb_rule(self.lprec) };
        match value {
            0 => BranchRule::Ceiling,
            1 => BranchRule::Floor,
            2 => BranchRule::Automatic,
            _ => BranchRule::Ceiling, // Default fallback
        }
    }

    /// Set the maximum depth for branch-and-bound.
    ///
    /// 0 means no limit.
    pub fn set_bb_depth_limit(&mut self, depth: libc::c_int) {
        unsafe { lp::set_bb_depthlimit(self.lprec, depth) };
    }

    /// Get the branch-and-bound depth limit.
    pub fn get_bb_depth_limit(&self) -> libc::c_int {
        unsafe { lp::get_bb_depthlimit(self.lprec) }
    }

    /// Set verbosity level for solver output.
    pub fn set_verbose(&mut self, level: Verbosity) {
        unsafe { lp::set_verbose(self.lprec, level as libc::c_int) };
    }

    /// Get the current verbosity level.
    pub fn get_verbose(&self) -> libc::c_int {
        unsafe { lp::get_verbose(self.lprec) }
    }

    /// Reset the basis to the default (all slack basis).
    pub fn default_basis(&mut self) {
        unsafe { lp::default_basis(self.lprec) };
    }

    /// Reset the basis and reinitialize.
    pub fn reset_basis(&mut self) {
        unsafe { lp::reset_basis(self.lprec) };
    }

    /// Read a basis from a file.
    ///
    /// The basis file is typically created by `write_basis`.
    pub fn read_basis(&mut self, filename: &CStr) -> Result<()> {
        unsafe {
            lp::read_basis(
                self.lprec,
                filename.as_ptr() as *mut _,
                std::ptr::null_mut(),
            )
        }.check_with(LpSolveError::BasisOperationFailed)
    }

    /// Write the current basis to a file.
    ///
    /// This allows warm-starting future solves with the same or similar problems.
    pub fn write_basis(&self, filename: &CStr) -> Result<()> {
        unsafe { lp::write_basis(self.lprec, filename.as_ptr() as *mut _) }
            .check_with(LpSolveError::BasisOperationFailed)
    }

    /// Set the basis crash mode.
    ///
    /// The crash mode determines how the initial basis is constructed.
    pub fn set_basiscrash(&mut self, mode: BasisCrash) {
        unsafe { lp::set_basiscrash(self.lprec, mode as libc::c_int) };
    }

    /// Get the current basis crash mode.
    pub fn get_basiscrash(&self) -> BasisCrash {
        let value = unsafe { lp::get_basiscrash(self.lprec) };
        match value {
            0 => BasisCrash::None,
            1 => BasisCrash::MostFeasible,
            2 => BasisCrash::LeastDegen,
            _ => BasisCrash::None, // Default fallback
        }
    }

    /// Enable efficient row addition mode.
    ///
    /// When true, constraints can be added more efficiently by deferring matrix reorganization.
    /// Call with false when done adding constraints.
    pub fn set_add_rowmode(&mut self, enabled: bool) -> Result<()> {
        unsafe { lp::set_add_rowmode(self.lprec, if enabled { 1 } else { 0 }) }
            .check_with(LpSolveError::InAddRowMode)
    }

    /// Check if add row mode is currently enabled.
    pub fn is_add_rowmode(&self) -> bool {
        unsafe { lp::is_add_rowmode(self.lprec) == 1 }
    }

    /// Set individual upper and lower bounds separately.
    ///
    /// This is an alternative to `set_bounds` when you only want to change one bound.
    pub fn set_upbo(&mut self, col: libc::c_int, value: f64) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetUpbo,
            });
        }
        unsafe { lp::set_upbo(self.lprec, col, value) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetUpbo,
            })
    }

    /// Get the upper bound of a variable.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [1, num_cols()].
    pub fn get_upbo(&self, col: libc::c_int) -> Result<f64> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            });
        }
        Ok(unsafe { lp::get_upbo(self.lprec, col) })
    }

    /// Set the lower bound of a variable.
    pub fn set_lowbo(&mut self, col: libc::c_int, value: f64) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetLowbo,
            });
        }
        unsafe { lp::set_lowbo(self.lprec, col, value) }
            .check_with(LpSolveError::ColumnOperationFailed {
                column: Some(col),
                context: OpContext::SetLowbo,
            })
    }

    /// Get the lower bound of a variable.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [1, num_cols()].
    pub fn get_lowbo(&self, col: libc::c_int) -> Result<f64> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            });
        }
        Ok(unsafe { lp::get_lowbo(self.lprec, col) })
    }

    /// Get a single matrix element.
    ///
    /// Row 0 is the objective function. Rows 1..num_rows() are constraints.
    /// Columns are 1-indexed from 1 to num_cols().
    ///
    /// # Errors
    /// Returns an error if indices are out of bounds.
    pub fn get_mat(&self, row: libc::c_int, col: libc::c_int) -> Result<f64> {
        if row < 0 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 0,
                max: self.num_rows(),
                context: OpContext::MatrixGet,
            });
        }
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::MatrixGet,
            });
        }
        Ok(unsafe { lp::get_mat(self.lprec, row, col) })
    }

    /// Set a single matrix element.
    ///
    /// Rows are 1-indexed from 1 to num_rows(). Columns are 1-indexed from 1 to num_cols().
    pub fn set_mat(&mut self, row: libc::c_int, col: libc::c_int, value: f64) -> Result<()> {
        if row < 0 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 0,
                max: self.num_rows(),
                context: OpContext::MatrixSet,
            });
        }
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::MatrixSet,
            });
        }
        unsafe { lp::set_mat(self.lprec, row, col, value) }
            .check_with(LpSolveError::MatrixOperationFailed {
                row,
                col,
                op: error::MatrixOp::Set,
            })
    }

    /// Get the number of non-zero elements in the matrix.
    pub fn get_nonzeros(&self) -> libc::c_int {
        unsafe { lp::get_nonzeros(self.lprec) }
    }

    /// Check if a given solution is feasible.
    ///
    /// The values slice must follow lpsolve's 1-indexed convention:
    /// - Must have exactly `num_cols() + 1` elements
    /// - Element at index 0 is unused (typically 0.0)
    /// - Elements at indices 1..=num_cols() contain the variable values
    ///
    /// # Errors
    /// Returns an error if `values.len()` is not exactly `num_cols() + 1`.
    pub fn is_feasible(&self, values: &mut [f64], threshold: f64) -> Result<bool> {
        let required = (self.num_cols() as usize).saturating_add(1);
        if values.len() != required {
            return Err(LpSolveError::DimensionMismatch {
                expected: required,
                actual: values.len(),
                context: OpContext::Validate,
            });
        }
        Ok(unsafe { lp::is_feasible(self.lprec, values.as_mut_ptr(), threshold) == 1 })
    }

    /// Get the number of binary variables in the model.
    ///
    /// If `working` is true, counts variables in the working basis, otherwise counts all binary variables.
    pub fn bin_count(&self, working: bool) -> libc::c_int {
        unsafe { lp::bin_count(self.lprec, if working { 1 } else { 0 }) }
    }

    /// Get the number of integer variables (including binary) in the model.
    pub fn mip_count(&self) -> libc::c_int {
        unsafe { lp::MIP_count(self.lprec) }
    }

    /// Get the number of Special Ordered Sets in the model.
    pub fn sos_count(&self) -> libc::c_int {
        unsafe { lp::SOS_count(self.lprec) }
    }

    /// Set anti-degeneracy strategy.
    ///
    /// Anti-degeneracy helps avoid cycling in the simplex algorithm.
    pub fn set_anti_degen(&mut self, strategy: AntiDegen) {
        unsafe { lp::set_anti_degen(self.lprec, strategy.bits()) };
    }

    /// Get the current anti-degeneracy strategy.
    pub fn get_anti_degen(&self) -> AntiDegen {
        let bits = unsafe { lp::get_anti_degen(self.lprec) };
        AntiDegen::from_bits_truncate(bits)
    }

    /// Set improvement strategy.
    ///
    /// Controls iterative improvement of the solution.
    pub fn set_improve(&mut self, mode: ImproveMode) {
        unsafe { lp::set_improve(self.lprec, mode.bits()) };
    }

    /// Get the current improvement strategy.
    pub fn get_improve(&self) -> ImproveMode {
        let bits = unsafe { lp::get_improve(self.lprec) };
        ImproveMode::from_bits_truncate(bits)
    }

    /// Set an objective value bound.
    ///
    /// The solver stops if it finds a solution with objective value better than this bound.
    pub fn set_obj_bound(&mut self, bound: f64) {
        unsafe { lp::set_obj_bound(self.lprec, bound) };
    }

    /// Get the objective bound.
    pub fn get_obj_bound(&self) -> f64 {
        unsafe { lp::get_obj_bound(self.lprec) }
    }

    /// Set a break-at-value for the objective function.
    ///
    /// The solver stops when it reaches this objective value.
    pub fn set_break_at_value(&mut self, value: f64) {
        unsafe { lp::set_break_at_value(self.lprec, value) };
    }

    /// Get the break-at-value.
    pub fn get_break_at_value(&self) -> f64 {
        unsafe { lp::get_break_at_value(self.lprec) }
    }

    /// Set whether to break at first feasible solution in MIP.
    pub fn set_break_at_first(&mut self, break_at_first: bool) {
        unsafe { lp::set_break_at_first(self.lprec, if break_at_first { 1 } else { 0 }) };
    }

    /// Check if break-at-first is enabled.
    pub fn is_break_at_first(&self) -> bool {
        unsafe { lp::is_break_at_first(self.lprec) == 1 }
    }

    /// Remove scaling from the model.
    ///
    /// Call this after solve() if you need to access the unscaled model.
    pub fn unscale(&mut self) {
        unsafe { lp::unscale(self.lprec) };
    }

    /// Construct a wrapper for a pre-existing `lprec`.
    ///
    /// # Safety
    ///
    /// The caller must ensure that:
    /// - `lprec` is a valid pointer to an initialized `lprec` structure
    /// - `lprec` has not been freed or deleted
    /// - The returned `Problem` becomes the sole owner of the `lprec` and will free it on drop
    /// - No other code will call `delete_lp` on this pointer
    pub unsafe fn from_lprec(lprec: *mut lp::lprec) -> Problem {
        debug_assert!(!lprec.is_null(), "from_lprec called with null pointer");
        Problem { lprec }
    }

    /// Get the raw `lprec` pointer.
    ///
    /// # Safety Note
    ///
    /// This is the internal pointer managed by this `Problem` instance.
    /// **DO NOT** call `delete_lp` on it - the `Problem` will free it on drop.
    /// Only use this for FFI interop or inspection.
    ///
    /// If you need to transfer ownership, use `into_lprec()` instead.
    pub fn to_lprec(&self) -> *mut lp::lprec {
        self.lprec
    }

    /// Consume the Problem and return the raw lprec pointer without freeing it.
    ///
    /// The caller becomes responsible for calling `delete_lp` on the returned pointer
    /// when done with it.
    ///
    /// # Example
    /// ```ignore
    /// let problem = Problem::new(2, 3).unwrap();
    /// let ptr = problem.into_lprec();
    /// // ... use ptr ...
    /// unsafe { lp::delete_lp(ptr); }  // Must manually free
    /// ```
    pub fn into_lprec(self) -> *mut lp::lprec {
        let ptr = self.lprec;
        std::mem::forget(self);  // Prevent drop
        ptr
    }

    pub fn num_cols(&self) -> libc::c_int {
        unsafe { lp::get_Ncolumns(self.lprec) }
    }

    pub fn num_rows(&self) -> libc::c_int {
        unsafe { lp::get_Nrows(self.lprec) }
    }

    /// Set the name of a row (constraint).
    ///
    /// Row 0 is the objective function. Rows 1..num_rows() are constraints.
    pub fn set_row_name(&mut self, row: libc::c_int, name: &CStr) -> Result<()> {
        if row < 0 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 0,
                max: self.num_rows(),
                context: OpContext::SetRowName,
            });
        }
        unsafe { lp::set_row_name(self.lprec, row, name.as_ptr() as *mut _) }
            .check_with(LpSolveError::NamingError {
                entity: EntityType::Row,
                index: row,
            })
    }

    /// Get the name of a row (constraint).
    ///
    /// Returns `Ok(None)` if the row has no name set.
    ///
    /// # Errors
    /// Returns an error if `row` is out of bounds [0, num_rows()].
    pub fn get_row_name(&self, row: libc::c_int) -> Result<Option<CString>> {
        if row < 0 || row > self.num_rows() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: row,
                min: 0,
                max: self.num_rows(),
                context: OpContext::GetRow,
            });
        }
        let ptr = unsafe { lp::get_row_name(self.lprec, row) };
        if ptr.is_null() {
            Ok(None)
        } else {
            Ok(unsafe { Some(CStr::from_ptr(ptr).to_owned()) })
        }
    }

    /// Set the name of a column (variable).
    ///
    /// Columns are numbered 1..num_cols().
    pub fn set_col_name(&mut self, col: libc::c_int, name: &CStr) -> Result<()> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::SetColName,
            });
        }
        unsafe { lp::set_col_name(self.lprec, col, name.as_ptr() as *mut _) }
            .check_with(LpSolveError::NamingError {
                entity: EntityType::Column,
                index: col,
            })
    }

    /// Get the name of a column (variable).
    ///
    /// Returns `Ok(None)` if the column has no name set.
    ///
    /// # Errors
    /// Returns an error if `col` is out of bounds [1, num_cols()].
    pub fn get_col_name(&self, col: libc::c_int) -> Result<Option<CString>> {
        if col < 1 || col > self.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.num_cols(),
                context: OpContext::GetColumn,
            });
        }
        let ptr = unsafe { lp::get_col_name(self.lprec, col) };
        if ptr.is_null() {
            Ok(None)
        } else {
            Ok(unsafe { Some(CStr::from_ptr(ptr).to_owned()) })
        }
    }

    /// Set the problem name.
    pub fn set_problem_name(&mut self, name: &CStr) -> Result<()> {
        unsafe { lp::set_lp_name(self.lprec, name.as_ptr() as *mut _) }
            .check_with(LpSolveError::NamingError {
                entity: EntityType::Problem,
                index: 0,
            })
    }

    /// Get the problem name.
    ///
    /// Returns `Ok(None)` if no name is set.
    pub fn get_problem_name(&self) -> Result<Option<CString>> {
        let ptr = unsafe { lp::get_lp_name(self.lprec) };
        if ptr.is_null() {
            Ok(None)
        } else {
            Ok(unsafe { Some(CStr::from_ptr(ptr).to_owned()) })
        }
    }
}

impl Drop for Problem {
    fn drop(&mut self) {
        unsafe { lp::delete_lp(self.lprec) }
    }
}

impl Clone for Problem {
    /// Clone the problem.
    ///
    /// # Panics
    ///
    /// Panics if memory allocation fails (out of memory).
    /// Use `try_clone()` for a non-panicking alternative.
    fn clone(&self) -> Problem {
        self.try_clone().expect("OOM when trying to copy_lp")
    }
}

impl Problem {
    /// Try to clone the problem, returning an error instead of panicking on OOM.
    pub fn try_clone(&self) -> Result<Problem> {
        let ptr = unsafe { lp::copy_lp(self.lprec) };
        if ptr.is_null() {
            Err(LpSolveError::OutOfMemory)
        } else {
            Ok(Problem { lprec: ptr })
        }
    }
}

unsafe impl Send for Problem {}

#[cfg(test)]
mod tests {
    use crate::{Problem, SolveStatus};
    #[test]
    fn smoke() {
        let mut lp = Problem::new(0, 0).unwrap();
        assert_eq!(lp.solve(), SolveStatus::NotRun);
    }
}