problemreductions/solvers/ilp/
solver.rs

1//! ILP solver implementation using HiGHS.
2
3use crate::models::algebraic::{Comparison, ObjectiveSense, VariableDomain, ILP};
4use crate::models::misc::TimetableDesign;
5use crate::rules::{ReduceTo, ReductionMode, ReductionResult};
6#[cfg(not(feature = "ilp-highs"))]
7use good_lp::default_solver;
8#[cfg(feature = "ilp-highs")]
9use good_lp::highs;
10#[cfg(feature = "ilp-highs")]
11use good_lp::solvers::highs::HighsParallelType;
12use good_lp::{variable, ProblemVariables, Solution, SolverModel, Variable};
13
14/// An ILP solver using the HiGHS backend.
15///
16/// This solver solves Integer Linear Programming problems directly using the HiGHS solver.
17///
18/// # Example
19///
20/// ```rust,ignore
21/// use problemreductions::models::algebraic::{ILP, LinearConstraint, ObjectiveSense};
22/// use problemreductions::solvers::ILPSolver;
23///
24/// // Create a simple binary ILP: maximize x0 + 2*x1 subject to x0 + x1 <= 1
25/// let ilp = ILP::<bool>::new(
26///     2,
27///     vec![LinearConstraint::le(vec![(0, 1.0), (1, 1.0)], 1.0)],
28///     vec![(0, 1.0), (1, 2.0)],
29///     ObjectiveSense::Maximize,
30/// );
31///
32/// let solver = ILPSolver::new();
33/// if let Some(solution) = solver.solve(&ilp) {
34///     println!("Solution: {:?}", solution);
35/// }
36/// ```
37#[derive(Debug, Clone, Default)]
38pub struct ILPSolver {
39    /// Time limit in seconds (None = no limit).
40    pub time_limit: Option<f64>,
41}
42
43#[derive(Debug, Clone, PartialEq, Eq)]
44pub enum SolveViaReductionError {
45    WitnessPathRequired { name: String },
46    NoReductionPath { name: String },
47    NoSolution { name: String },
48}
49
50impl std::fmt::Display for SolveViaReductionError {
51    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
52        match self {
53            SolveViaReductionError::WitnessPathRequired { name } => write!(
54                f,
55                "ILP solving requires a witness-capable source problem and reduction path; only aggregate-value solving is available for {}.",
56                name
57            ),
58            SolveViaReductionError::NoReductionPath { name } => {
59                write!(f, "No reduction path from {} to ILP", name)
60            }
61            SolveViaReductionError::NoSolution { name } => {
62                write!(f, "ILP solver found no solution for {}", name)
63            }
64        }
65    }
66}
67
68impl std::error::Error for SolveViaReductionError {}
69
70impl ILPSolver {
71    /// Create a new ILP solver with default settings.
72    pub fn new() -> Self {
73        Self::default()
74    }
75
76    /// Create an ILP solver with a time limit.
77    pub fn with_time_limit(seconds: f64) -> Self {
78        Self {
79            time_limit: Some(seconds),
80        }
81    }
82
83    /// Solve an ILP problem directly.
84    ///
85    /// Returns `None` if the problem is infeasible or the solver fails.
86    /// The returned solution is a configuration vector where each element
87    /// is the variable value (config index = value).
88    pub fn solve<V: VariableDomain>(&self, problem: &ILP<V>) -> Option<Vec<usize>> {
89        let n = problem.num_vars;
90        if n == 0 {
91            return problem.is_feasible(&[]).then_some(vec![]);
92        }
93
94        // Derive tighter per-variable upper bounds from single-variable ≤ constraints.
95        // This avoids giving HiGHS the full domain (e.g. 2^31 for i32), which can
96        // cause severe performance degradation even when constraints already bound
97        // the variable to a small range.
98        let default_ub = (V::DIMS_PER_VAR - 1) as f64;
99        let mut upper_bounds = vec![default_ub; n];
100        for constraint in &problem.constraints {
101            if constraint.cmp == crate::models::algebraic::Comparison::Le
102                && constraint.terms.len() == 1
103            {
104                let (var_idx, coef) = constraint.terms[0];
105                if coef > 0.0 && var_idx < n {
106                    let ub = constraint.rhs / coef;
107                    if ub < upper_bounds[var_idx] {
108                        upper_bounds[var_idx] = ub;
109                    }
110                }
111            }
112        }
113
114        // Create integer variables with tightened bounds
115        let mut vars_builder = ProblemVariables::new();
116        let vars: Vec<Variable> = (0..n)
117            .map(|i| {
118                let mut v = variable().integer();
119                v = v.min(0.0);
120                v = v.max(upper_bounds[i]);
121                vars_builder.add(v)
122            })
123            .collect();
124
125        // Build objective expression
126        let objective: good_lp::Expression = problem
127            .objective
128            .iter()
129            .map(|&(var_idx, coef)| coef * vars[var_idx])
130            .sum();
131
132        // Build the model with objective
133        let unsolved = match problem.sense {
134            ObjectiveSense::Maximize => vars_builder.maximise(&objective),
135            ObjectiveSense::Minimize => vars_builder.minimise(&objective),
136        };
137
138        // Create the solver model
139        #[cfg(feature = "ilp-highs")]
140        let mut model = {
141            let mut model = unsolved
142                .using(highs)
143                .set_option("random_seed", 0i32)
144                .set_option("presolve", "off")
145                .set_parallel(HighsParallelType::Off)
146                .set_threads(1);
147            if let Some(seconds) = self.time_limit {
148                model = model.set_time_limit(seconds);
149            }
150            model
151        };
152
153        #[cfg(not(feature = "ilp-highs"))]
154        let mut model = unsolved.using(default_solver);
155
156        // Add constraints
157        for constraint in &problem.constraints {
158            // Build left-hand side expression
159            let lhs: good_lp::Expression = constraint
160                .terms
161                .iter()
162                .map(|&(var_idx, coef)| coef * vars[var_idx])
163                .sum();
164
165            // Create the constraint based on comparison type
166            let good_lp_constraint = match constraint.cmp {
167                Comparison::Le => lhs.leq(constraint.rhs),
168                Comparison::Ge => lhs.geq(constraint.rhs),
169                Comparison::Eq => lhs.eq(constraint.rhs),
170            };
171
172            model = model.with(good_lp_constraint);
173        }
174
175        // Solve
176        let solution = model.solve().ok()?;
177
178        // Extract solution: config index = value (no lower bound offset)
179        let result: Vec<usize> = vars
180            .iter()
181            .map(|v| {
182                let val = solution.value(*v);
183                val.round().max(0.0) as usize
184            })
185            .collect();
186
187        Some(result)
188    }
189
190    /// Solve any problem that reduces to `ILP<bool>`.
191    ///
192    /// This method first reduces the problem to a binary ILP, solves the ILP,
193    /// and then extracts the solution back to the original problem space.
194    ///
195    /// # Example
196    ///
197    /// ```no_run
198    /// use problemreductions::prelude::*;
199    /// use problemreductions::solvers::ILPSolver;
200    ///
201    /// // Create a problem that reduces directly to ILP.
202    /// let problem = MaximumSetPacking::<i32>::new(vec![
203    ///     vec![0, 1],
204    ///     vec![1, 2],
205    ///     vec![3, 4],
206    /// ]);
207    ///
208    /// // Solve using ILP solver
209    /// let solver = ILPSolver::new();
210    /// if let Some(solution) = solver.solve_reduced(&problem) {
211    ///     println!("Solution: {:?}", solution);
212    /// }
213    /// ```
214    pub fn solve_reduced<P>(&self, problem: &P) -> Option<Vec<usize>>
215    where
216        P: ReduceTo<ILP<bool>>,
217    {
218        let reduction = problem.reduce_to();
219        let ilp_solution = self.solve(reduction.target_problem())?;
220        Some(reduction.extract_solution(&ilp_solution))
221    }
222
223    /// Solve a type-erased problem directly when a native solver hook exists.
224    ///
225    /// Returns `None` if the input type has no direct solver or the solver finds no solution.
226    pub fn solve_dyn(&self, any: &dyn std::any::Any) -> Option<Vec<usize>> {
227        if let Some(ilp) = any.downcast_ref::<ILP<bool>>() {
228            return self.solve(ilp);
229        }
230        if let Some(ilp) = any.downcast_ref::<ILP<i32>>() {
231            return self.solve(ilp);
232        }
233        if let Some(problem) = any.downcast_ref::<TimetableDesign>() {
234            return problem.solve_via_required_assignments();
235        }
236        None
237    }
238
239    fn supports_direct_dyn(&self, any: &dyn std::any::Any) -> bool {
240        any.is::<ILP<bool>>() || any.is::<ILP<i32>>() || any.is::<TimetableDesign>()
241    }
242
243    /// Two-level path selection:
244    /// 1. Dijkstra finds the cheapest path to each ILP variant using
245    ///    `MinimizeStepsThenOverhead` (additive edge costs: step count + log overhead).
246    /// 2. Across ILP variants, we pick the path whose composed final output size
247    ///    is smallest — this is the actual ILP problem size the solver will face.
248    fn best_path_to_ilp(
249        &self,
250        graph: &crate::rules::ReductionGraph,
251        name: &str,
252        variant: &std::collections::BTreeMap<String, String>,
253        mode: ReductionMode,
254        instance: &dyn std::any::Any,
255    ) -> Option<crate::rules::ReductionPath> {
256        let ilp_variants = graph.variants_for("ILP");
257        let input_size = crate::rules::ReductionGraph::compute_source_size(name, instance);
258        let mut best_path: Option<crate::rules::ReductionPath> = None;
259        let mut best_cost = f64::INFINITY;
260
261        for dv in &ilp_variants {
262            if let Some(path) = graph.find_cheapest_path_mode(
263                name,
264                variant,
265                "ILP",
266                dv,
267                mode,
268                &input_size,
269                &crate::rules::MinimizeStepsThenOverhead,
270            ) {
271                // Use composed final output size for cross-variant comparison,
272                // since this determines the actual ILP problem size.
273                let final_size = graph
274                    .evaluate_path_overhead(&path, &input_size)
275                    .unwrap_or_default();
276                let cost = final_size.total() as f64;
277                if cost < best_cost {
278                    best_cost = cost;
279                    best_path = Some(path);
280                }
281            }
282        }
283
284        best_path
285    }
286
287    pub fn try_solve_via_reduction(
288        &self,
289        name: &str,
290        variant: &std::collections::BTreeMap<String, String>,
291        instance: &dyn std::any::Any,
292    ) -> Result<Vec<usize>, SolveViaReductionError> {
293        if self.supports_direct_dyn(instance) {
294            return self
295                .solve_dyn(instance)
296                .ok_or_else(|| SolveViaReductionError::NoSolution {
297                    name: name.to_string(),
298                });
299        }
300
301        let graph = crate::rules::ReductionGraph::new();
302
303        let Some(path) =
304            self.best_path_to_ilp(&graph, name, variant, ReductionMode::Witness, instance)
305        else {
306            if self
307                .best_path_to_ilp(&graph, name, variant, ReductionMode::Aggregate, instance)
308                .is_some()
309            {
310                return Err(SolveViaReductionError::WitnessPathRequired {
311                    name: name.to_string(),
312                });
313            }
314
315            return Err(SolveViaReductionError::NoReductionPath {
316                name: name.to_string(),
317            });
318        };
319
320        let chain = graph.reduce_along_path(&path, instance).ok_or_else(|| {
321            SolveViaReductionError::WitnessPathRequired {
322                name: name.to_string(),
323            }
324        })?;
325        let ilp_solution = self.solve_dyn(chain.target_problem_any()).ok_or_else(|| {
326            SolveViaReductionError::NoSolution {
327                name: name.to_string(),
328            }
329        })?;
330        Ok(chain.extract_solution(&ilp_solution))
331    }
332
333    /// Solve a type-erased problem by finding a reduction path to ILP.
334    ///
335    /// Tries all ILP variants, picks the cheapest path, reduces, solves,
336    /// and extracts the solution back. Falls back to direct ILP solve if
337    /// the problem is already an ILP type.
338    ///
339    /// Returns `None` if no path to ILP exists or the solver finds no solution.
340    pub fn solve_via_reduction(
341        &self,
342        name: &str,
343        variant: &std::collections::BTreeMap<String, String>,
344        instance: &dyn std::any::Any,
345    ) -> Option<Vec<usize>> {
346        self.try_solve_via_reduction(name, variant, instance).ok()
347    }
348}
349
350#[cfg(test)]
351#[path = "../../unit_tests/solvers/ilp/solver.rs"]
352mod tests;
problemreductions/solvers/ilp/solver.rs

problemreductions/solvers/ilp/
solver.rs
