problemreductions/rules/
travelingsalesman_qubo.rs

1//! Reduction from TravelingSalesman to QUBO.
2//!
3//! Uses the standard position-based QUBO encoding for TSP:
4//! - Binary variables x_{v,p} = 1 iff vertex v is at position p in the tour
5//! - H_A: each vertex appears exactly once (row constraint)
6//! - H_B: each position has exactly one vertex (column constraint)
7//! - H_C: objective encoding edge costs between consecutive positions
8
9use crate::models::algebraic::QUBO;
10use crate::models::graph::TravelingSalesman;
11use crate::reduction;
12use crate::rules::traits::{ReduceTo, ReductionResult};
13use crate::topology::{Graph, SimpleGraph};
14use std::collections::HashMap;
15
16/// Result of reducing TravelingSalesman to QUBO.
17#[derive(Debug, Clone)]
18pub struct ReductionTravelingSalesmanToQUBO {
19    target: QUBO<i64>,
20    num_vertices: usize,
21    num_edges: usize,
22    edge_index: HashMap<(usize, usize), usize>,
23}
24
25impl ReductionResult for ReductionTravelingSalesmanToQUBO {
26    type Source = TravelingSalesman<SimpleGraph, i64>;
27    type Target = QUBO<i64>;
28
29    fn target_problem(&self) -> &Self::Target {
30        &self.target
31    }
32
33    /// Decode position encoding back to edge-based configuration.
34    ///
35    /// The QUBO solution uses n^2 binary variables x_{v,p} (vertex v at position p).
36    /// We extract the tour order, then map consecutive pairs to edge indices.
37    fn extract_solution(
38        &self,
39        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
40    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
41        crate::rules::traits::validate_target_solution(self.target_problem(), target_solution)?;
42
43        Ok({
44            let n = self.num_vertices;
45
46            let tour: Vec<usize> = (0..n)
47                .map(|position| {
48                    let mut selected =
49                        (0..n).filter(|&vertex| target_solution[vertex * n + position]);
50                    match (selected.next(), selected.next()) {
51                        (Some(vertex), None) => Ok(vertex),
52                        _ => Err(crate::rules::ExtractionError::invalid(format!(
53                            "tour position {position} does not select exactly one vertex"
54                        ))),
55                    }
56                })
57                .collect::<crate::rules::ExtractionResult<_>>()?;
58
59            // Build edge-based config: for each consecutive pair in the tour, mark the edge
60            let mut config = vec![false; self.num_edges];
61            for p in 0..n {
62                let u = tour[p];
63                let v = tour[(p + 1) % n];
64                let key = (u.min(v), u.max(v));
65                let &edge = self.edge_index.get(&key).ok_or_else(|| {
66                    crate::rules::ExtractionError::invalid(format!(
67                        "target tour uses absent source edge ({u}, {v})"
68                    ))
69                })?;
70                config[edge] = true;
71            }
72
73            config
74        })
75    }
76}
77
78#[reduction(
79    transform = exact {
80        num_vars = "num_vertices^2",
81    }
82)]
83impl ReduceTo<QUBO<i64>> for TravelingSalesman<SimpleGraph, i64> {
84    type Result = ReductionTravelingSalesmanToQUBO;
85
86    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
87        let n = self.num_vertices();
88        let edges = self.edges();
89
90        // Build edge weight map (both directions for undirected lookup)
91        let overflow = |operation| {
92            crate::rules::ReductionError::integer_overflow::<
93                TravelingSalesman<SimpleGraph, i64>,
94                QUBO<i64>,
95            >(operation)
96        };
97        let mut edge_weight_map: HashMap<(usize, usize), i64> = HashMap::new();
98        let mut weight_sum = 0i64;
99        for &(u, v, w) in &edges {
100            edge_weight_map.insert((u, v), w);
101            edge_weight_map.insert((v, u), w);
102            let magnitude = w
103                .checked_abs()
104                .ok_or_else(|| overflow("taking the absolute value of a tour weight"))?;
105            weight_sum = weight_sum
106                .checked_add(magnitude)
107                .ok_or_else(|| overflow("summing absolute tour weights"))?;
108        }
109
110        // Build edge index map: canonical (min, max) → edge index
111        let graph_edges = self.graph().edges();
112        let num_edges = graph_edges.len();
113        let mut edge_index: HashMap<(usize, usize), usize> = HashMap::new();
114        for (idx, &(u, v)) in graph_edges.iter().enumerate() {
115            edge_index.insert((u.min(v), u.max(v)), idx);
116        }
117
118        // Penalty weight: must exceed any possible tour cost
119        let a = weight_sum
120            .checked_add(1)
121            .ok_or_else(|| overflow("computing the tour penalty"))?;
122
123        // Build n^2 x n^2 upper-triangular QUBO matrix
124        let dim = n
125            .checked_mul(n)
126            .ok_or_else(|| overflow("computing the number of QUBO variables"))?;
127        let mut matrix = vec![vec![0i64; dim]; dim];
128
129        // Helper: add value to upper-triangular position
130        let mut add_upper = |i: usize, j: usize, val: i64| {
131            let (lo, hi) = if i <= j { (i, j) } else { (j, i) };
132            matrix[lo][hi] = matrix[lo][hi]
133                .checked_add(val)
134                .ok_or_else(|| overflow("adding a tour QUBO coefficient"))?;
135            Ok::<(), crate::rules::ReductionError>(())
136        };
137
138        // H_A: each vertex visited exactly once (row constraint)
139        // For each vertex v: (sum_p x_{v,p} - 1)^2
140        // = sum_p x_{v,p}^2 - 2*sum_p x_{v,p} + 1
141        // = -sum_p x_{v,p} + 2*sum_{p1<p2} x_{v,p1}*x_{v,p2} + const
142        for v in 0..n {
143            for p in 0..n {
144                // Diagonal: -A (from expanding (sum - 1)^2, the -2*x + x^2 = -x for binary)
145                add_upper(
146                    v * n + p,
147                    v * n + p,
148                    a.checked_neg()
149                        .ok_or_else(|| overflow("negating the tour penalty"))?,
150                )?;
151            }
152            for p1 in 0..n {
153                for p2 in (p1 + 1)..n {
154                    // Cross terms: 2*A * x_{v,p1} * x_{v,p2}
155                    add_upper(
156                        v * n + p1,
157                        v * n + p2,
158                        a.checked_mul(2)
159                            .ok_or_else(|| overflow("doubling the tour penalty"))?,
160                    )?;
161                }
162            }
163        }
164
165        // H_B: each position has exactly one vertex (column constraint)
166        // For each position p: (sum_v x_{v,p} - 1)^2
167        for p in 0..n {
168            for v in 0..n {
169                add_upper(
170                    v * n + p,
171                    v * n + p,
172                    a.checked_neg()
173                        .ok_or_else(|| overflow("negating the tour penalty"))?,
174                )?;
175            }
176            for v1 in 0..n {
177                for v2 in (v1 + 1)..n {
178                    add_upper(
179                        v1 * n + p,
180                        v2 * n + p,
181                        a.checked_mul(2)
182                            .ok_or_else(|| overflow("doubling the tour penalty"))?,
183                    )?;
184                }
185            }
186        }
187
188        // H_C: distance objective
189        // For each pair (u, v), add cost for x_{u,p} * x_{v,p_next} and x_{v,p} * x_{u,p_next}
190        for u in 0..n {
191            for v in (u + 1)..n {
192                let cost = edge_weight_map.get(&(u, v)).copied().unwrap_or(a);
193                for p in 0..n {
194                    let p_next = (p + 1) % n;
195                    // x_{u,p} * x_{v,p_next}
196                    add_upper(u * n + p, v * n + p_next, cost)?;
197                    // x_{v,p} * x_{u,p_next}
198                    add_upper(v * n + p, u * n + p_next, cost)?;
199                }
200            }
201        }
202
203        let target = QUBO::from_matrix(matrix).map_err(|message| {
204            crate::rules::ReductionError::construction::<
205                TravelingSalesman<SimpleGraph, i64>,
206                QUBO<i64>,
207            >(message)
208        })?;
209
210        Ok(ReductionTravelingSalesmanToQUBO {
211            target,
212            num_vertices: n,
213            num_edges,
214            edge_index,
215        })
216    }
217}
218
219#[cfg(feature = "example-db")]
220pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
221    use crate::export::SolutionPair;
222    use crate::models::algebraic::QUBO;
223
224    vec![crate::example_db::specs::RuleExampleSpec {
225        id: "travelingsalesman_to_qubo",
226        build: || {
227            let source = TravelingSalesman::new(
228                SimpleGraph::new(3, vec![(0, 1), (0, 2), (1, 2)]),
229                vec![1, 2, 3],
230            );
231            crate::example_db::specs::rule_example_with_witness::<_, QUBO<i64>>(
232                source,
233                SolutionPair {
234                    source_config: serde_json::json!(vec![true, true, true]),
235                    target_config: serde_json::json!(vec![
236                        false, false, true, true, false, false, false, true, false
237                    ]),
238                },
239            )
240        },
241    }]
242}
243
244#[cfg(test)]
245#[path = "../unit_tests/rules/travelingsalesman_qubo.rs"]
246mod tests;
problemreductions/rules/travelingsalesman_qubo.rs

problemreductions/rules/
travelingsalesman_qubo.rs
