lpsolve 1.0.1

use crate::{
    Problem, ConstraintType, SolveStatus, SimplexType,
    PresolveMode, ScalingMode, PivotRule, BranchRule, AntiDegen,
    ImproveMode, BasisCrash, Verbosity, MPSOptions,
    error::{LpSolveError, Result, FileFormat, OpContext}
};
use std::ffi::CString;
use std::cell::OnceCell;
use std::path::Path;

/// Builder for creating and configuring LP problems with a fluent API
///
/// # Example
/// ```ignore
/// let solution = Problem::builder()
///     .rows(2)
///     .cols(3)
///     .maximize()
///     .objective(&[10.0, 20.0, 30.0])
///     .constraint(&[1.0, 2.0, 3.0], ConstraintType::Le, 100.0)?
///     .constraint(&[4.0, 5.0, 6.0], ConstraintType::Ge, 50.0)?
///     .integer_variable(1)
///     .bounds(2, 0.0, 50.0)?
///     .solve()?;
///
/// println!("Optimal value: {}", solution.objective_value());
/// ```
pub struct ProblemBuilder {
    rows: Option<i32>,
    cols: Option<i32>,
    constraints: Vec<PendingConstraint>,
    objective: Option<Vec<f64>>,
    maximize: bool,
    variable_configs: Vec<VariableConfig>,
    solver_config: SolverConfig,
    name: Option<CString>,
}

struct PendingConstraint {
    coefficients: Vec<f64>,
    rhs: f64,
    constraint_type: ConstraintType,
    name: Option<CString>,
}

struct VariableConfig {
    column: i32,
    is_integer: bool,
    is_binary: bool,
    is_unbounded: bool,
    is_semicont: bool,
    lower_bound: Option<f64>,
    upper_bound: Option<f64>,
    name: Option<CString>,
}

/// Configuration for the solver
#[derive(Default)]
pub struct SolverConfig {
    pub presolve: Option<PresolveMode>,
    pub scaling: Option<ScalingMode>,
    pub simplex_type: Option<SimplexType>,
    pub pivot_rule: Option<PivotRule>,
    pub branch_rule: Option<BranchRule>,
    pub anti_degen: Option<AntiDegen>,
    pub improve: Option<ImproveMode>,
    pub basis_crash: Option<BasisCrash>,
    pub timeout: Option<i64>,
    pub verbosity: Option<Verbosity>,
    pub mip_gap_abs: Option<f64>,
    pub mip_gap_rel: Option<f64>,
    pub epsint: Option<f64>,
    pub epsb: Option<f64>,
    pub epsd: Option<f64>,
    pub epspivot: Option<f64>,
    pub bb_depth_limit: Option<i32>,
    pub obj_bound: Option<f64>,
    pub break_at_value: Option<f64>,
    pub break_at_first: Option<bool>,
    pub prefer_dual: Option<bool>,
}

impl ProblemBuilder {
    /// Create a new problem builder
    pub fn new() -> Self {
        Self {
            rows: None,
            cols: None,
            constraints: Vec::new(),
            objective: None,
            maximize: false,
            variable_configs: Vec::new(),
            solver_config: SolverConfig::default(),
            name: None,
        }
    }

    /// Load a problem from an LP format file
    ///
    /// Returns the loaded Problem directly, ready to be configured further or solved.
    pub fn from_lp_file<P: AsRef<Path>>(path: P, verbosity: Verbosity) -> Result<Problem> {
        let path_cstr = CString::new(path.as_ref().to_string_lossy().as_bytes())
            .map_err(|_| LpSolveError::FileReadError { format: FileFormat::LP })?;
        let initial_name = CString::new("").unwrap();

        Problem::read_lp(&path_cstr, verbosity, &initial_name)
            .ok_or(LpSolveError::FileReadError { format: FileFormat::LP })
    }

    /// Load a problem from a free MPS format file
    ///
    /// Returns the loaded Problem directly, ready to be configured further or solved.
    pub fn from_freemps_file<P: AsRef<Path>>(path: P, options: MPSOptions) -> Result<Problem> {
        let path_cstr = CString::new(path.as_ref().to_string_lossy().as_bytes())
            .map_err(|_| LpSolveError::FileReadError { format: FileFormat::FreeMPS })?;

        Problem::read_freemps(&path_cstr, options)
            .ok_or(LpSolveError::FileReadError { format: FileFormat::FreeMPS })
    }

    /// Load a problem from a fixed MPS format file
    ///
    /// Returns the loaded Problem directly, ready to be configured further or solved.
    pub fn from_fixedmps_file<P: AsRef<Path>>(path: P, options: MPSOptions) -> Result<Problem> {
        let path_cstr = CString::new(path.as_ref().to_string_lossy().as_bytes())
            .map_err(|_| LpSolveError::FileReadError { format: FileFormat::FixedMPS })?;

        Problem::read_fixedmps(&path_cstr, options)
            .ok_or(LpSolveError::FileReadError { format: FileFormat::FixedMPS })
    }

    /// Set the initial number of rows (constraints)
    pub fn rows(mut self, rows: i32) -> Self {
        self.rows = Some(rows);
        self
    }

    /// Set the initial number of columns (variables)
    pub fn cols(mut self, cols: i32) -> Self {
        self.cols = Some(cols);
        self
    }

    /// Set the problem name
    ///
    /// Returns an error if the name contains null bytes.
    pub fn name<S: AsRef<str>>(mut self, name: S) -> Result<Self> {
        self.name = Some(CString::new(name.as_ref())
            .map_err(|_| LpSolveError::InvalidName { context: "problem name" })?);
        Ok(self)
    }

    /// Set the objective function coefficients
    ///
    /// The slice should NOT include the constant term (index 0).
    /// For a problem with n variables, pass n coefficients.
    ///
    /// # Panics
    /// Panics if `cols()` has been set and the coefficient count doesn't match.
    /// Use `try_objective()` for a non-panicking alternative.
    pub fn objective(self, coeffs: &[f64]) -> Self {
        self.try_objective(coeffs).expect("objective coefficient count mismatch")
    }

    /// Set the objective function coefficients with validation
    ///
    /// Returns an error immediately if `cols()` has been set and dimensions don't match.
    pub fn try_objective(mut self, coeffs: &[f64]) -> Result<Self> {
        // Early validation if cols is known
        if let Some(cols) = self.cols {
            if coeffs.len() != cols as usize {
                return Err(LpSolveError::DimensionMismatch {
                    expected: cols as usize,
                    actual: coeffs.len(),
                    context: OpContext::SetObjective,
                });
            }
        }
        self.objective = Some(coeffs.to_vec());
        Ok(self)
    }

    /// Set optimization to maximize (default is minimize)
    pub fn maximize(mut self) -> Self {
        self.maximize = true;
        self
    }

    /// Set optimization to minimize (this is the default)
    pub fn minimize(mut self) -> Self {
        self.maximize = false;
        self
    }

    /// Add a constraint to the problem
    ///
    /// The coefficients should NOT include the RHS constant (index 0).
    ///
    /// # Panics
    /// Panics if `cols()` has been set and the coefficient count doesn't match.
    /// Use `try_constraint()` for a non-panicking alternative.
    pub fn constraint(
        self,
        coeffs: &[f64],
        constraint_type: ConstraintType,
        rhs: f64,
    ) -> Self {
        self.try_constraint(coeffs, constraint_type, rhs)
            .expect("constraint coefficient count mismatch")
    }

    /// Add a constraint to the problem with validation
    ///
    /// Returns an error immediately if `cols()` has been set and dimensions don't match.
    pub fn try_constraint(
        mut self,
        coeffs: &[f64],
        constraint_type: ConstraintType,
        rhs: f64,
    ) -> Result<Self> {
        // Early validation if cols is known
        if let Some(cols) = self.cols {
            if coeffs.len() != cols as usize {
                return Err(LpSolveError::DimensionMismatch {
                    expected: cols as usize,
                    actual: coeffs.len(),
                    context: OpContext::AddConstraint,
                });
            }
        }
        self.constraints.push(PendingConstraint {
            coefficients: coeffs.to_vec(),
            rhs,
            constraint_type,
            name: None,
        });
        Ok(self)
    }

    /// Add a named constraint
    ///
    /// Returns an error if the name contains null bytes or if dimensions don't match.
    pub fn named_constraint<S: AsRef<str>>(
        mut self,
        name: S,
        coeffs: &[f64],
        constraint_type: ConstraintType,
        rhs: f64,
    ) -> Result<Self> {
        // Early validation if cols is known
        if let Some(cols) = self.cols {
            if coeffs.len() != cols as usize {
                return Err(LpSolveError::DimensionMismatch {
                    expected: cols as usize,
                    actual: coeffs.len(),
                    context: OpContext::AddConstraint,
                });
            }
        }
        self.constraints.push(PendingConstraint {
            coefficients: coeffs.to_vec(),
            rhs,
            constraint_type,
            name: Some(CString::new(name.as_ref())
                .map_err(|_| LpSolveError::InvalidName { context: "constraint name" })?),
        });
        Ok(self)
    }

    /// Configure a variable as binary
    pub fn binary_variable(mut self, col: i32) -> Self {
        self.ensure_var_config(col).is_binary = true;
        self
    }

    /// Configure a variable as integer
    pub fn integer_variable(mut self, col: i32) -> Self {
        self.ensure_var_config(col).is_integer = true;
        self
    }

    /// Configure a variable as unbounded
    pub fn unbounded_variable(mut self, col: i32) -> Self {
        self.ensure_var_config(col).is_unbounded = true;
        self
    }

    /// Configure a variable as semi-continuous
    pub fn semicontinuous_variable(mut self, col: i32) -> Self {
        self.ensure_var_config(col).is_semicont = true;
        self
    }

    /// Set bounds for a variable
    pub fn bounds(mut self, col: i32, lower: f64, upper: f64) -> Self {
        let config = self.ensure_var_config(col);
        config.lower_bound = Some(lower);
        config.upper_bound = Some(upper);
        self
    }

    /// Set lower bound for a variable
    pub fn lower_bound(mut self, col: i32, lower: f64) -> Self {
        self.ensure_var_config(col).lower_bound = Some(lower);
        self
    }

    /// Set upper bound for a variable
    pub fn upper_bound(mut self, col: i32, upper: f64) -> Self {
        self.ensure_var_config(col).upper_bound = Some(upper);
        self
    }

    /// Name a variable
    ///
    /// Returns an error if the name contains null bytes.
    pub fn variable_name<S: AsRef<str>>(mut self, col: i32, name: S) -> Result<Self> {
        self.ensure_var_config(col).name = Some(CString::new(name.as_ref())
            .map_err(|_| LpSolveError::InvalidName { context: "variable name" })?);
        Ok(self)
    }

    /// Name all variables at once
    ///
    /// The slice should have one name per column. Empty strings will be skipped.
    /// Returns an error if any name contains null bytes.
    pub fn variable_names<S: AsRef<str>>(mut self, names: &[S]) -> Result<Self> {
        for (i, name) in names.iter().enumerate() {
            let name_str = name.as_ref();
            if !name_str.is_empty() {
                let col = (i + 1) as i32; // Variables are 1-indexed
                self.ensure_var_config(col).name = Some(CString::new(name_str)
                    .map_err(|_| LpSolveError::InvalidName { context: "variable name" })?);
            }
        }
        Ok(self)
    }

    /// Name a constraint (must be called after the constraint is added)
    ///
    /// Returns an error if the name contains null bytes or if the row is out of bounds.
    pub fn constraint_name<S: AsRef<str>>(mut self, row: i32, name: S) -> Result<Self> {
        if let Some(constraint) = self.constraints.get_mut((row - 1) as usize) {
            constraint.name = Some(CString::new(name.as_ref())
                .map_err(|_| LpSolveError::InvalidName { context: "constraint name" })?);
        }
        Ok(self)
    }

    /// Name all constraints at once
    ///
    /// The slice should have one name per constraint. Empty strings will be skipped.
    /// This should be called after all constraints are added.
    /// Returns an error if any name contains null bytes.
    pub fn constraint_names<S: AsRef<str>>(mut self, names: &[S]) -> Result<Self> {
        for (i, name) in names.iter().enumerate() {
            let name_str = name.as_ref();
            if !name_str.is_empty() && i < self.constraints.len() {
                self.constraints[i].name = Some(CString::new(name_str)
                    .map_err(|_| LpSolveError::InvalidName { context: "constraint name" })?);
            }
        }
        Ok(self)
    }

    /// Get or create a variable configuration, maintaining sorted order for O(log n) lookups.
    ///
    /// The variable_configs Vec is kept sorted by column number to enable binary search.
    fn ensure_var_config(&mut self, col: i32) -> &mut VariableConfig {
        match self.variable_configs.binary_search_by_key(&col, |v| v.column) {
            Ok(idx) => &mut self.variable_configs[idx],
            Err(insert_pos) => {
                // Insert new config at the correct position to maintain sorted order
                self.variable_configs.insert(insert_pos, VariableConfig {
                    column: col,
                    is_integer: false,
                    is_binary: false,
                    is_unbounded: false,
                    is_semicont: false,
                    lower_bound: None,
                    upper_bound: None,
                    name: None,
                });
                &mut self.variable_configs[insert_pos]
            }
        }
    }

    // Solver configuration methods

    /// Set presolve mode
    pub fn presolve(mut self, mode: PresolveMode) -> Self {
        self.solver_config.presolve = Some(mode);
        self
    }

    /// Set scaling mode
    pub fn scaling(mut self, mode: ScalingMode) -> Self {
        self.solver_config.scaling = Some(mode);
        self
    }

    /// Set timeout in seconds
    pub fn timeout(mut self, seconds: i64) -> Self {
        self.solver_config.timeout = Some(seconds);
        self
    }

    /// Set verbosity level
    pub fn verbosity(mut self, level: Verbosity) -> Self {
        self.solver_config.verbosity = Some(level);
        self
    }

    /// Set branch-and-bound rule for MIP problems
    pub fn branch_rule(mut self, rule: BranchRule) -> Self {
        self.solver_config.branch_rule = Some(rule);
        self
    }

    /// Configure for typical LP problem
    pub fn typical_lp(self) -> Self {
        self.presolve(PresolveMode::TYPICAL)
            .scaling(ScalingMode::TYPICAL)
    }

    /// Configure for typical MIP problem
    pub fn typical_mip(self) -> Self {
        self.presolve(PresolveMode::TYPICAL | PresolveMode::REDUCEMIP)
            .scaling(ScalingMode::TYPICAL)
            .branch_rule(BranchRule::Automatic)
    }

    // ===== Terse/Shorthand API Methods =====

    /// Shorthand for `.maximize().objective(coeffs)`
    ///
    /// # Example
    /// ```ignore
    /// Problem::builder().cols(2).max(&[30.0, 50.0])
    /// ```
    pub fn max(self, coeffs: &[f64]) -> Self {
        self.maximize().objective(coeffs)
    }

    /// Shorthand for `.minimize().objective(coeffs)`
    ///
    /// # Example
    /// ```ignore
    /// Problem::builder().cols(2).min(&[1.0, 2.0])
    /// ```
    pub fn min(self, coeffs: &[f64]) -> Self {
        self.minimize().objective(coeffs)
    }

    /// Shorthand for adding a less-than-or-equal constraint
    ///
    /// # Example
    /// ```ignore
    /// builder.le(&[2.0, 3.0], 100.0)  // 2x1 + 3x2 <= 100
    /// ```
    pub fn le(self, coeffs: &[f64], rhs: f64) -> Self {
        self.constraint(coeffs, ConstraintType::Le, rhs)
    }

    /// Shorthand for adding a greater-than-or-equal constraint
    ///
    /// # Example
    /// ```ignore
    /// builder.ge(&[1.0, 1.0], 5.0)  // x1 + x2 >= 5
    /// ```
    pub fn ge(self, coeffs: &[f64], rhs: f64) -> Self {
        self.constraint(coeffs, ConstraintType::Ge, rhs)
    }

    /// Shorthand for adding an equality constraint
    ///
    /// # Example
    /// ```ignore
    /// builder.eq(&[1.0, -1.0], 0.0)  // x1 - x2 = 0
    /// ```
    pub fn eq(self, coeffs: &[f64], rhs: f64) -> Self {
        self.constraint(coeffs, ConstraintType::Eq, rhs)
    }

    /// Set all variables to have a lower bound of 0 (non-negative variables)
    ///
    /// This is a common requirement in LP/MIP problems.
    ///
    /// # Example
    /// ```ignore
    /// builder.cols(5).non_negative()  // All 5 variables >= 0
    /// ```
    pub fn non_negative(mut self) -> Self {
        let cols = self.cols.unwrap_or(0);
        for col in 1..=cols {
            self = self.lower_bound(col, 0.0);
        }
        self
    }

    /// Set multiple variables as binary (0 or 1)
    ///
    /// # Example
    /// ```ignore
    /// builder.binary_vars(&[1, 2, 3])  // First 3 vars are binary
    /// ```
    pub fn binary_vars(mut self, cols: &[i32]) -> Self {
        for &col in cols {
            self = self.binary_variable(col);
        }
        self
    }

    /// Set multiple variables as integer
    ///
    /// # Example
    /// ```ignore
    /// builder.integer_vars(&[1, 2])  // x1 and x2 must be integers
    /// ```
    pub fn integer_vars(mut self, cols: &[i32]) -> Self {
        for &col in cols {
            self = self.integer_variable(col);
        }
        self
    }

    /// Set the same bounds for all variables
    ///
    /// # Example
    /// ```ignore
    /// builder.cols(3).bounds_all(0.0, 100.0)  // All vars in [0, 100]
    /// ```
    pub fn bounds_all(mut self, lower: f64, upper: f64) -> Self {
        let cols = self.cols.unwrap_or(0);
        for col in 1..=cols {
            self = self.bounds(col, lower, upper);
        }
        self
    }

    /// Set the same lower bound for all variables
    ///
    /// # Example
    /// ```ignore
    /// builder.cols(3).lower_bound_all(0.0)  // Equivalent to non_negative()
    /// ```
    pub fn lower_bound_all(mut self, lower: f64) -> Self {
        let cols = self.cols.unwrap_or(0);
        for col in 1..=cols {
            self = self.lower_bound(col, lower);
        }
        self
    }

    /// Set the same upper bound for all variables
    ///
    /// # Example
    /// ```ignore
    /// builder.cols(3).upper_bound_all(100.0)  // All vars <= 100
    /// ```
    pub fn upper_bound_all(mut self, upper: f64) -> Self {
        let cols = self.cols.unwrap_or(0);
        for col in 1..=cols {
            self = self.upper_bound(col, upper);
        }
        self
    }

    /// Set all variables as non-negative integers
    ///
    /// This is a common requirement in many optimization problems.
    ///
    /// # Example
    /// ```ignore
    /// builder.cols(5).non_negative_integers()  // All 5 variables are integers >= 0
    /// ```
    pub fn non_negative_integers(mut self) -> Self {
        let cols = self.cols.unwrap_or(0);
        for col in 1..=cols {
            self = self.integer_variable(col).lower_bound(col, 0.0);
        }
        self
    }

    /// Clone solver configuration from another problem
    ///
    /// This is useful when you want to solve a similar problem with the same settings.
    ///
    /// # Example
    /// ```ignore
    /// let new_problem = Problem::builder()
    ///     .cols(3)
    ///     .clone_config_from(&existing_problem)
    ///     .solve()?;
    /// ```
    pub fn clone_config_from(mut self, other: &Problem) -> Self {
        self.solver_config.presolve = Some(other.get_presolve());
        self.solver_config.scaling = Some(other.get_scaling());
        self.solver_config.pivot_rule = Some(other.get_pivoting());
        self.solver_config.branch_rule = Some(other.get_bb_rule());
        self.solver_config.anti_degen = Some(other.get_anti_degen());
        self.solver_config.improve = Some(other.get_improve());
        self.solver_config.basis_crash = Some(other.get_basiscrash());
        self.solver_config.timeout = Some(other.get_timeout());
        // Note: get_verbose() returns i32, not Verbosity, so we skip it
        self.solver_config.epsint = Some(other.get_epsint());
        self.solver_config.epsb = Some(other.get_epsb());
        self.solver_config.epsd = Some(other.get_epsd());
        self.solver_config.epspivot = Some(other.get_epspivot());
        self.solver_config.bb_depth_limit = Some(other.get_bb_depth_limit());
        self.solver_config.obj_bound = Some(other.get_obj_bound());
        self.solver_config.break_at_value = Some(other.get_break_at_value());
        self.solver_config.break_at_first = Some(other.is_break_at_first());
        self.solver_config.mip_gap_abs = Some(other.get_mip_gap(true));
        self.solver_config.mip_gap_rel = Some(other.get_mip_gap(false));
        self
    }

    /// Validate builder consistency
    fn validate(&self) -> Result<()> {
        let cols = self.cols.unwrap_or(0);

        // Validate objective function dimension
        if let Some(ref obj) = self.objective {
            if obj.len() != cols as usize {
                return Err(LpSolveError::DimensionMismatch {
                    expected: cols as usize,
                    actual: obj.len(),
                    context: OpContext::SetObjective,
                });
            }
        }

        // Validate constraint dimensions
        for constraint in &self.constraints {
            if constraint.coefficients.len() != cols as usize {
                return Err(LpSolveError::DimensionMismatch {
                    expected: cols as usize,
                    actual: constraint.coefficients.len(),
                    context: OpContext::AddConstraint,
                });
            }
        }

        // Validate variable configurations are within bounds
        for config in &self.variable_configs {
            if config.column < 1 || config.column > cols {
                return Err(LpSolveError::IndexOutOfBounds {
                    index: config.column,
                    min: 1,
                    max: cols,
                    context: OpContext::SetColumn,
                });
            }
        }

        Ok(())
    }

    /// Build the problem but don't solve it yet
    pub fn build(self) -> Result<Problem> {
        // Validate before building
        self.validate()?;

        let rows = self.rows.unwrap_or(0);
        let cols = self.cols.unwrap_or(0);

        let mut problem = Problem::new(rows, cols)
            .ok_or(LpSolveError::CreationFailed)?;

        // Set problem name if provided
        if let Some(ref name) = self.name {
            problem.set_problem_name(name)?;
        }

        // Set optimization direction
        if self.maximize {
            problem.set_maximize();
        } else {
            problem.set_minimize();
        }

        // Set objective function if provided
        if let Some(ref obj) = self.objective {
            // Need to prepend 0.0 for the constant term
            let mut full_obj = vec![0.0];
            full_obj.extend_from_slice(obj);
            problem.set_objective_function(&mut full_obj)?;
        }

        // Apply solver configuration first (before consuming self)
        self.apply_solver_config(&mut problem)?;

        // Add constraints
        for (i, constraint) in self.constraints.into_iter().enumerate() {
            // Need to prepend 0.0 for the constant term
            let mut full_coeffs = vec![0.0];
            full_coeffs.extend_from_slice(&constraint.coefficients);

            problem.add_constraint(&mut full_coeffs, constraint.rhs, constraint.constraint_type)?;

            if let Some(name) = constraint.name {
                problem.set_row_name((i + 1) as i32, &name)?;
            }
        }

        // Configure variables
        for config in self.variable_configs.into_iter() {
            if config.is_binary {
                problem.set_binary(config.column, true)?;
            } else if config.is_integer {
                problem.set_integer(config.column, true)?;
            }

            if config.is_unbounded {
                problem.set_unbounded(config.column)?;
            }

            if config.is_semicont {
                problem.set_semicont(config.column, true)?;
            }

            if let (Some(lower), Some(upper)) = (config.lower_bound, config.upper_bound) {
                problem.set_bounds(config.column, lower, upper)?;
            } else if let Some(lower) = config.lower_bound {
                problem.set_lowbo(config.column, lower)?;
            } else if let Some(upper) = config.upper_bound {
                problem.set_upbo(config.column, upper)?;
            }

            if let Some(name) = config.name {
                problem.set_col_name(config.column, &name)?;
            }
        }

        Ok(problem)
    }

    /// Build and solve the problem, returning the solution
    pub fn solve(self) -> Result<Solution> {
        let mut problem = self.build()?;
        let status = problem.solve();

        Ok(Solution::new(problem, status))
    }

    fn apply_solver_config(&self, problem: &mut Problem) -> Result<()> {
        let config = &self.solver_config;

        if let Some(presolve) = config.presolve {
            problem.set_presolve(presolve);
        }

        if let Some(scaling) = config.scaling {
            problem.set_scaling(scaling);
        }

        if let Some(simplex_type) = config.simplex_type {
            problem.set_simplex_type(simplex_type);
        }

        if let Some(pivot_rule) = config.pivot_rule {
            problem.set_pivoting(pivot_rule);
        }

        if let Some(branch_rule) = config.branch_rule {
            problem.set_bb_rule(branch_rule);
        }

        if let Some(anti_degen) = config.anti_degen {
            problem.set_anti_degen(anti_degen);
        }

        if let Some(improve) = config.improve {
            problem.set_improve(improve);
        }

        if let Some(basis_crash) = config.basis_crash {
            problem.set_basiscrash(basis_crash);
        }

        if let Some(timeout) = config.timeout {
            problem.set_timeout(timeout);
        }

        if let Some(verbosity) = config.verbosity {
            problem.set_verbose(verbosity);
        }

        if let Some(gap) = config.mip_gap_abs {
            problem.set_mip_gap(true, gap);
        }

        if let Some(gap) = config.mip_gap_rel {
            problem.set_mip_gap(false, gap);
        }

        if let Some(epsint) = config.epsint {
            problem.set_epsint(epsint);
        }

        if let Some(epsb) = config.epsb {
            problem.set_epsb(epsb);
        }

        if let Some(epsd) = config.epsd {
            problem.set_epsd(epsd);
        }

        if let Some(epspivot) = config.epspivot {
            problem.set_epspivot(epspivot);
        }

        if let Some(limit) = config.bb_depth_limit {
            problem.set_bb_depth_limit(limit);
        }

        if let Some(bound) = config.obj_bound {
            problem.set_obj_bound(bound);
        }

        if let Some(value) = config.break_at_value {
            problem.set_break_at_value(value);
        }

        if let Some(break_first) = config.break_at_first {
            problem.set_break_at_first(break_first);
        }

        if let Some(prefer_dual) = config.prefer_dual {
            problem.set_prefer_dual(prefer_dual);
        }

        Ok(())
    }
}

impl Default for ProblemBuilder {
    fn default() -> Self {
        Self::new()
    }
}

/// Represents a solved problem with convenient access to solution data
pub struct Solution {
    problem: Problem,
    status: SolveStatus,
    /// Cached variable values (computed once, reused for all variable() calls)
    variables: Option<Vec<f64>>,
    /// Cached dual values (computed on first access via interior mutability)
    duals: OnceCell<Option<Vec<f64>>>,
    /// Cached constraint values (computed on first access via interior mutability)
    constraints: OnceCell<Option<Vec<f64>>>,
}

impl Solution {
    pub(crate) fn new(problem: Problem, status: SolveStatus) -> Self {
        // Pre-compute variables if feasible to avoid repeated allocations
        let variables = if matches!(
            status,
            SolveStatus::Optimal | SolveStatus::Suboptimal | SolveStatus::FeasibleFound
        ) {
            let n = problem.num_cols() as usize;
            let mut vars = vec![0.0; n];
            if problem.get_solution_variables(&mut vars).is_some() {
                Some(vars)
            } else {
                None
            }
        } else {
            None
        };

        Self {
            problem,
            status,
            variables,
            duals: OnceCell::new(),
            constraints: OnceCell::new(),
        }
    }

    /// Get the solve status
    pub fn status(&self) -> SolveStatus {
        self.status
    }

    /// Check if the solution is optimal
    pub fn is_optimal(&self) -> bool {
        self.status == SolveStatus::Optimal
    }

    /// Check if a feasible solution was found
    pub fn is_feasible(&self) -> bool {
        matches!(
            self.status,
            SolveStatus::Optimal | SolveStatus::Suboptimal | SolveStatus::FeasibleFound
        )
    }

    /// Get the objective function value
    pub fn objective_value(&self) -> f64 {
        self.problem.get_objective()
    }

    /// Get the values of all variables
    ///
    /// Returns None if no feasible solution exists
    pub fn variables(&self) -> Option<&[f64]> {
        self.variables.as_deref()
    }

    /// Get the value of a specific variable (1-indexed)
    ///
    /// This is now very efficient as variables are cached when Solution is created.
    pub fn variable(&self, col: i32) -> Option<f64> {
        let vars = self.variables.as_ref()?;
        if col < 1 || col > self.problem.num_cols() {
            return None;
        }
        Some(vars[(col - 1) as usize])
    }

    /// Get the value of a specific variable (1-indexed) with proper error handling
    ///
    /// This is now very efficient as variables are cached when Solution is created.
    ///
    /// # Errors
    /// Returns an error if:
    /// - No feasible solution exists
    /// - The column index is out of bounds
    pub fn get_variable(&self, col: i32) -> Result<f64> {
        if !self.is_feasible() {
            return Err(LpSolveError::NoSolutionAvailable);
        }
        let vars = self.variables.as_ref()
            .ok_or(LpSolveError::NoSolutionAvailable)?;
        if col < 1 || col > self.problem.num_cols() {
            return Err(LpSolveError::IndexOutOfBounds {
                index: col,
                min: 1,
                max: self.problem.num_cols(),
                context: OpContext::GetColumn,
            });
        }
        Ok(vars[(col - 1) as usize])
    }

    /// Access variables through a visitor pattern
    ///
    /// Since variables are now cached, this is just as efficient as direct access
    pub fn visit_variables<F>(&self, mut visitor: F) -> Option<()>
    where
        F: FnMut(usize, f64),
    {
        let vars = self.variables.as_ref()?;
        for (i, &value) in vars.iter().enumerate() {
            visitor(i, value);
        }
        Some(())
    }

    /// Get dual values (shadow prices)
    ///
    /// Computed once and cached for subsequent calls.
    /// Uses interior mutability so it can be called with `&self`.
    pub fn dual_values(&self) -> Option<&[f64]> {
        if !self.is_feasible() {
            return None;
        }

        // Compute and cache if not already done
        self.duals.get_or_init(|| {
            let rows = self.problem.num_rows() as usize;
            let cols = self.problem.num_cols() as usize;
            let mut duals = vec![0.0; rows + cols];
            if self.problem.get_dual_solution(&mut duals).is_some() {
                Some(duals)
            } else {
                None
            }
        }).as_ref().map(|v| v.as_slice())
    }

    /// Get constraint values
    ///
    /// Computed once and cached for subsequent calls.
    /// Uses interior mutability so it can be called with `&self`.
    pub fn constraint_values(&self) -> Option<&[f64]> {
        if !self.is_feasible() {
            return None;
        }

        // Compute and cache if not already done
        self.constraints.get_or_init(|| {
            let rows = (self.problem.num_rows() + 1) as usize;
            let mut constraints = vec![0.0; rows];
            if self.problem.get_constraints(&mut constraints).is_some() {
                Some(constraints)
            } else {
                None
            }
        }).as_ref().map(|v| v.as_slice())
    }

    /// Get total iterations performed
    pub fn iterations(&self) -> i64 {
        self.problem.get_total_iter()
    }

    /// Get total nodes evaluated (for MIP problems)
    pub fn nodes(&self) -> i64 {
        self.problem.get_total_nodes()
    }

    /// Get time elapsed during solve
    pub fn time_elapsed(&self) -> f64 {
        self.problem.time_elapsed()
    }

    /// Get the underlying problem for further operations
    pub fn problem(&self) -> &Problem {
        &self.problem
    }

    /// Get the underlying problem mutably
    pub fn problem_mut(&mut self) -> &mut Problem {
        &mut self.problem
    }

    /// Consume the solution and return the underlying problem
    pub fn into_problem(self) -> Problem {
        self.problem
    }
}