problemreductions/models/optimization/
qubo.rs

1//! QUBO (Quadratic Unconstrained Binary Optimization) problem implementation.
2//!
3//! QUBO minimizes a quadratic function over binary variables.
4
5use crate::graph_types::SimpleGraph;
6use crate::traits::Problem;
7use crate::types::{EnergyMode, ProblemSize, SolutionSize};
8use serde::{Deserialize, Serialize};
9
10/// The QUBO (Quadratic Unconstrained Binary Optimization) problem.
11///
12/// Given n binary variables x_i ∈ {0, 1} and a matrix Q,
13/// minimize the quadratic form:
14///
15/// f(x) = Σ_i Σ_j Q_ij * x_i * x_j = x^T Q x
16///
17/// The matrix Q is typically upper triangular, with diagonal elements
18/// representing linear terms and off-diagonal elements representing
19/// quadratic interactions.
20///
21/// # Example
22///
23/// ```
24/// use problemreductions::models::optimization::QUBO;
25/// use problemreductions::{Problem, Solver, BruteForce};
26///
27/// // Q matrix: minimize x0 - 2*x1 + x0*x1
28/// // Q = [[1, 1], [0, -2]]
29/// let problem = QUBO::from_matrix(vec![
30///     vec![1.0, 1.0],
31///     vec![0.0, -2.0],
32/// ]);
33///
34/// let solver = BruteForce::new();
35/// let solutions = solver.find_best(&problem);
36///
37/// // Optimal is x = [0, 1] with value -2
38/// assert!(solutions.contains(&vec![0, 1]));
39/// ```
40#[derive(Debug, Clone, Serialize, Deserialize)]
41pub struct QUBO<W = f64> {
42    /// Number of variables.
43    num_vars: usize,
44    /// Q matrix stored as upper triangular (row-major).
45    /// `Q[i][j]` for i <= j represents the coefficient of x_i * x_j
46    matrix: Vec<Vec<W>>,
47}
48
49impl<W: Clone + Default> QUBO<W> {
50    /// Create a QUBO problem from a full matrix.
51    ///
52    /// The matrix should be square. Only the upper triangular part
53    /// (including diagonal) is used.
54    pub fn from_matrix(matrix: Vec<Vec<W>>) -> Self {
55        let num_vars = matrix.len();
56        Self { num_vars, matrix }
57    }
58
59    /// Create a QUBO from linear and quadratic terms.
60    ///
61    /// # Arguments
62    /// * `linear` - Linear coefficients (diagonal of Q)
63    /// * `quadratic` - Quadratic coefficients as ((i, j), value) for i < j
64    pub fn new(linear: Vec<W>, quadratic: Vec<((usize, usize), W)>) -> Self
65    where
66        W: num_traits::Zero,
67    {
68        let num_vars = linear.len();
69        let mut matrix = vec![vec![W::zero(); num_vars]; num_vars];
70
71        // Set diagonal (linear terms)
72        for (i, val) in linear.into_iter().enumerate() {
73            matrix[i][i] = val;
74        }
75
76        // Set off-diagonal (quadratic terms)
77        for ((i, j), val) in quadratic {
78            if i < j {
79                matrix[i][j] = val;
80            } else {
81                matrix[j][i] = val;
82            }
83        }
84
85        Self { num_vars, matrix }
86    }
87
88    /// Get the number of variables.
89    pub fn num_vars(&self) -> usize {
90        self.num_vars
91    }
92
93    /// Get the Q matrix.
94    pub fn matrix(&self) -> &[Vec<W>] {
95        &self.matrix
96    }
97
98    /// Get a specific matrix element `Q[i][j]`.
99    pub fn get(&self, i: usize, j: usize) -> Option<&W> {
100        self.matrix.get(i).and_then(|row| row.get(j))
101    }
102}
103
104impl<W> Problem for QUBO<W>
105where
106    W: Clone
107        + Default
108        + PartialOrd
109        + num_traits::Num
110        + num_traits::Zero
111        + std::ops::AddAssign
112        + std::ops::Mul<Output = W>
113        + 'static,
114{
115    const NAME: &'static str = "QUBO";
116    type GraphType = SimpleGraph;
117    type Weight = W;
118    type Size = W;
119
120    fn num_variables(&self) -> usize {
121        self.num_vars
122    }
123
124    fn num_flavors(&self) -> usize {
125        2 // Binary
126    }
127
128    fn problem_size(&self) -> ProblemSize {
129        ProblemSize::new(vec![("num_vars", self.num_vars)])
130    }
131
132    fn energy_mode(&self) -> EnergyMode {
133        EnergyMode::SmallerSizeIsBetter // Minimize
134    }
135
136    fn solution_size(&self, config: &[usize]) -> SolutionSize<Self::Size> {
137        let value = self.evaluate(config);
138        SolutionSize::valid(value)
139    }
140}
141
142impl<W> QUBO<W>
143where
144    W: Clone + num_traits::Zero + std::ops::AddAssign + std::ops::Mul<Output = W>,
145{
146    /// Evaluate the QUBO objective for a configuration.
147    pub fn evaluate(&self, config: &[usize]) -> W {
148        let mut value = W::zero();
149
150        for i in 0..self.num_vars {
151            let x_i = config.get(i).copied().unwrap_or(0);
152            if x_i == 0 {
153                continue;
154            }
155
156            for j in i..self.num_vars {
157                let x_j = config.get(j).copied().unwrap_or(0);
158                if x_j == 0 {
159                    continue;
160                }
161
162                if let Some(q_ij) = self.matrix.get(i).and_then(|row| row.get(j)) {
163                    value += q_ij.clone();
164                }
165            }
166        }
167
168        value
169    }
170}
171
172#[cfg(test)]
173mod tests {
174    use super::*;
175    use crate::solvers::{BruteForce, Solver};
176
177    #[test]
178    fn test_qubo_from_matrix() {
179        let problem = QUBO::from_matrix(vec![vec![1.0, 2.0], vec![0.0, 3.0]]);
180        assert_eq!(problem.num_vars(), 2);
181        assert_eq!(problem.get(0, 0), Some(&1.0));
182        assert_eq!(problem.get(0, 1), Some(&2.0));
183        assert_eq!(problem.get(1, 1), Some(&3.0));
184    }
185
186    #[test]
187    fn test_qubo_new() {
188        let problem = QUBO::new(vec![1.0, 2.0], vec![((0, 1), 3.0)]);
189        assert_eq!(problem.get(0, 0), Some(&1.0));
190        assert_eq!(problem.get(1, 1), Some(&2.0));
191        assert_eq!(problem.get(0, 1), Some(&3.0));
192    }
193
194    #[test]
195    fn test_evaluate() {
196        // Q = [[1, 2], [0, 3]]
197        // f(x) = x0 + 3*x1 + 2*x0*x1
198        let problem = QUBO::from_matrix(vec![vec![1.0, 2.0], vec![0.0, 3.0]]);
199
200        assert_eq!(problem.evaluate(&[0, 0]), 0.0);
201        assert_eq!(problem.evaluate(&[1, 0]), 1.0);
202        assert_eq!(problem.evaluate(&[0, 1]), 3.0);
203        assert_eq!(problem.evaluate(&[1, 1]), 6.0); // 1 + 3 + 2 = 6
204    }
205
206    #[test]
207    fn test_solution_size() {
208        let problem = QUBO::from_matrix(vec![vec![1.0, 2.0], vec![0.0, 3.0]]);
209
210        let sol = problem.solution_size(&[0, 0]);
211        assert!(sol.is_valid);
212        assert_eq!(sol.size, 0.0);
213
214        let sol = problem.solution_size(&[1, 1]);
215        assert_eq!(sol.size, 6.0);
216    }
217
218    #[test]
219    fn test_brute_force_minimize() {
220        // Q = [[1, 0], [0, -2]]
221        // f(x) = x0 - 2*x1
222        // Minimum at x = [0, 1] with value -2
223        let problem = QUBO::from_matrix(vec![vec![1.0, 0.0], vec![0.0, -2.0]]);
224        let solver = BruteForce::new();
225
226        let solutions = solver.find_best(&problem);
227        assert_eq!(solutions.len(), 1);
228        assert_eq!(solutions[0], vec![0, 1]);
229        assert_eq!(problem.solution_size(&solutions[0]).size, -2.0);
230    }
231
232    #[test]
233    fn test_brute_force_with_interaction() {
234        // Q = [[-1, 2], [0, -1]]
235        // f(x) = -x0 - x1 + 2*x0*x1
236        // x=[0,0] -> 0, x=[1,0] -> -1, x=[0,1] -> -1, x=[1,1] -> 0
237        let problem = QUBO::from_matrix(vec![vec![-1.0, 2.0], vec![0.0, -1.0]]);
238        let solver = BruteForce::new();
239
240        let solutions = solver.find_best(&problem);
241        // Minimum is -1 at [1,0] or [0,1]
242        assert_eq!(solutions.len(), 2);
243        for sol in &solutions {
244            assert_eq!(problem.solution_size(sol).size, -1.0);
245        }
246    }
247
248    #[test]
249    fn test_energy_mode() {
250        let problem = QUBO::<f64>::from_matrix(vec![vec![1.0]]);
251        assert!(problem.energy_mode().is_minimization());
252    }
253
254    #[test]
255    fn test_num_variables_flavors() {
256        let problem = QUBO::<f64>::from_matrix(vec![vec![0.0; 5]; 5]);
257        assert_eq!(problem.num_variables(), 5);
258        assert_eq!(problem.num_flavors(), 2);
259    }
260
261    #[test]
262    fn test_problem_size() {
263        let problem = QUBO::<f64>::from_matrix(vec![vec![0.0; 3]; 3]);
264        let size = problem.problem_size();
265        assert_eq!(size.get("num_vars"), Some(3));
266    }
267
268    #[test]
269    fn test_matrix_access() {
270        let problem = QUBO::from_matrix(vec![
271            vec![1.0, 2.0, 3.0],
272            vec![0.0, 4.0, 5.0],
273            vec![0.0, 0.0, 6.0],
274        ]);
275        let matrix = problem.matrix();
276        assert_eq!(matrix.len(), 3);
277        assert_eq!(matrix[0], vec![1.0, 2.0, 3.0]);
278    }
279
280    #[test]
281    fn test_empty_qubo() {
282        let problem = QUBO::<f64>::from_matrix(vec![]);
283        assert_eq!(problem.num_vars(), 0);
284        assert_eq!(problem.evaluate(&[]), 0.0);
285    }
286
287    #[test]
288    fn test_single_variable() {
289        let problem = QUBO::from_matrix(vec![vec![-5.0]]);
290        let solver = BruteForce::new();
291
292        let solutions = solver.find_best(&problem);
293        assert_eq!(solutions.len(), 1);
294        assert_eq!(solutions[0], vec![1]); // x=1 gives -5, x=0 gives 0
295    }
296
297    #[test]
298    fn test_qubo_new_reverse_indices() {
299        // Test the case where (j, i) is provided with i < j
300        let problem = QUBO::new(vec![1.0, 2.0], vec![((1, 0), 3.0)]); // j > i
301        assert_eq!(problem.get(0, 1), Some(&3.0)); // Should be stored at (0, 1)
302    }
303
304    #[test]
305    fn test_get_out_of_bounds() {
306        let problem = QUBO::from_matrix(vec![vec![1.0, 2.0], vec![0.0, 3.0]]);
307        assert_eq!(problem.get(5, 5), None);
308        assert_eq!(problem.get(0, 5), None);
309    }
310}
problemreductions/models/optimization/qubo.rs

problemreductions/models/optimization/
qubo.rs
