hercules 0.5.0

A Heuristics toolbox for QUBO in Rust
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

use crate::branch_node::QuboBBNode;
use crate::branch_subproblem::SubProblemSolver;
use crate::branch_subproblem::{SubProblemOptions, SubProblemResult};
use crate::branchbound::BBSolver;
use crate::preprocess::make_sub_problem;
use crate::qubo::Qubo;
use clarabel::algebra::CscMatrix;
use clarabel::solver::{DefaultSettings, DefaultSolver, IPSolver, NonnegativeConeT};
use ndarray::Array1;
use sprs::{CsMat, TriMat};

#[derive(Clone)]
pub struct ClarabelQPSolver {}

impl SubProblemSolver for ClarabelQPSolver {
    fn solve_lower_bound(
        &self,
        bbsolver: &BBSolver,
        node: &QuboBBNode,
        _: Option<SubProblemOptions>,
    ) -> SubProblemResult {
        // solve QP associated with the node
        // generate default settings
        let settings = DefaultSettings {
            verbose: false,
            ..Default::default()
        };

        // if we are fully fixed, then just return the solution
        if node.fixed_variables.len() == bbsolver.qubo.num_x() {
            let mut sol = Array1::zeros(bbsolver.qubo.num_x());
            for (&i, &val) in &node.fixed_variables {
                sol[i] = val as f64;
            }
            let obj = bbsolver.qubo.eval(&sol);
            return (obj, sol);
        }

        // find projected subproblem
        let (sub_qubo, unfixed_map, _constant) =
            make_sub_problem(&bbsolver.qubo, &node.fixed_variables);

        // generate the constraint matrix
        let A_size = 2 * sub_qubo.num_x();
        let mut A = TriMat::new((A_size, sub_qubo.num_x()));
        let mut b = Array1::<f64>::zeros(A_size);

        // add the inequality constraints
        for (index, i) in (0..sub_qubo.num_x()).enumerate() {
            let offset = index * 2;
            A.add_triplet(offset, i, 1.0);
            A.add_triplet(offset + 1, i, -1.0);
            b[offset] = 1.0;
            b[offset + 1] = 0.0;
        }

        // convert the matrix to CSC format and then Clarabel format
        let A_csc = A.to_csc();
        let A_clara = Self::make_cb_form(&A_csc);

        // generate the cones for the solver
        let cones = [NonnegativeConeT(2 * sub_qubo.num_x())];

        // set up the solver with the matrices
        let mut solver = DefaultSolver::new(
            &Self::make_cb_form(&sub_qubo.q),
            sub_qubo.c.as_slice().unwrap(), // unwrap is safe because Array1 is stored in contiguous memory
            &A_clara,
            b.as_slice().unwrap(),
            &cones,
            settings,
        );

        // solve the optimization problem
        solver.solve();

        // convert the solution back to the original space
        let mut x = Array1::<f64>::zeros(bbsolver.qubo.num_x());

        // map out the unfixed variables
        for (&original, &new) in &unfixed_map {
            x[original] = solver.solution.x[new];
        }

        // map out the fixed variables
        for (&i, &val) in &node.fixed_variables {
            x[i] = val as f64;
        }

        let obj = bbsolver.qubo.eval(&x);
        (obj, x)
    }
}

impl ClarabelQPSolver {
    pub fn new(_: &Qubo) -> Self {
        Self {}
    }
    pub fn make_cb_form(p0: &CsMat<f64>) -> CscMatrix {
        let (t, y, u) = p0.to_csc().into_raw_storage();
        CscMatrix::new(p0.rows(), p0.cols(), t, y, u)
    }
}

#[cfg(test)]
mod tests {
    use crate::subproblemsolvers::clarabel_qp::ClarabelQPSolver;
    use crate::tests::make_solver_qubo;

    #[test]
    fn ensure_matrix_equivalence() {
        // the idea of this test is to make sure that the conversion from sprs to clarabel is correct

        let qubo = make_solver_qubo();

        let clarabel_matrix = ClarabelQPSolver::make_cb_form(&qubo.q);

        // check the number of non-zero elements
        assert_eq!(qubo.q.nnz(), clarabel_matrix.nnz());

        // check the values of the non-zero elements
        for (&val, (i, j)) in qubo.q.iter() {
            assert_eq!(val, clarabel_matrix.get_entry((i, j)).unwrap());
        }
    }
}