problemreductions/rules/
minimumdiscreteplanarinversekinematics_qubo.rs

1//! Reduction from MinimumDiscretePlanarInverseKinematics to QUBO.
2//!
3//! Each sampled orientation becomes a binary one-hot variable. The planar
4//! end-effector coordinates are linear in those selectors, so the squared
5//! distance to the target expands directly into a quadratic objective. One-hot
6//! penalties enforce exactly one orientation per link, and forbidden-pair
7//! penalties enforce the consecutive-joint feasibility relations.
8//!
9//! Reference: Salloum et al. (2025).
10
11use crate::models::algebraic::QUBO;
12use crate::models::misc::MinimumDiscretePlanarInverseKinematics;
13use crate::reduction;
14use crate::rules::traits::{ReduceTo, ReductionResult};
15
16fn block_offsets(block_sizes: &[usize]) -> Vec<usize> {
17    let mut offsets = Vec::with_capacity(block_sizes.len());
18    let mut offset = 0;
19    for &size in block_sizes {
20        offsets.push(offset);
21        offset += size;
22    }
23    offsets
24}
25
26/// Result of reducing MinimumDiscretePlanarInverseKinematics to QUBO.
27#[derive(Debug, Clone)]
28pub struct ReductionMinimumDiscretePlanarInverseKinematicsToQUBO {
29    target: QUBO<f64>,
30    block_offsets: Vec<usize>,
31    block_sizes: Vec<usize>,
32    allowed_pairs: Vec<Vec<(usize, usize)>>,
33}
34
35impl ReductionResult for ReductionMinimumDiscretePlanarInverseKinematicsToQUBO {
36    type Source = MinimumDiscretePlanarInverseKinematics;
37    type Target = QUBO<f64>;
38
39    fn target_problem(&self) -> &Self::Target {
40        &self.target
41    }
42
43    fn extract_solution(
44        &self,
45        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
46    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
47        crate::rules::traits::validate_target_solution(self.target_problem(), target_solution)?;
48
49        let config: Vec<usize> = self
50            .block_offsets
51            .iter()
52            .zip(&self.block_sizes)
53            .enumerate()
54            .map(|(link, (&start, &size))| {
55                let mut selected = target_solution[start..start + size]
56                    .iter()
57                    .enumerate()
58                    .filter_map(|(orientation, &bit)| bit.then_some(orientation));
59                match (selected.next(), selected.next()) {
60                    (Some(orientation), None) => Ok(orientation),
61                    (None, _) => Err(crate::rules::ExtractionError::invalid(format!(
62                        "link {link} has no selected orientation"
63                    ))),
64                    (Some(_), Some(_)) => Err(crate::rules::ExtractionError::invalid(format!(
65                        "link {link} has multiple selected orientations"
66                    ))),
67                }
68            })
69            .collect::<crate::rules::ExtractionResult<_>>()?;
70        for (junction, (pair, allowed)) in config.windows(2).zip(&self.allowed_pairs).enumerate() {
71            if !allowed.contains(&(pair[0], pair[1])) {
72                return Err(crate::rules::ExtractionError::invalid(format!(
73                    "junction {junction} has a forbidden orientation pair"
74                )));
75            }
76        }
77        Ok(config)
78    }
79}
80
81#[reduction(transform = {
82    exact {
83        num_vars = "num_orientation_samples",
84    },
85    upper_bound {
86        num_quadratic_terms = "num_orientation_samples * (num_orientation_samples - 1) / 2",
87    },
88})]
89impl ReduceTo<QUBO<f64>> for MinimumDiscretePlanarInverseKinematics {
90    type Result = ReductionMinimumDiscretePlanarInverseKinematicsToQUBO;
91
92    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
93        let block_sizes: Vec<usize> = self.orientation_samples().iter().map(Vec::len).collect();
94        let block_offsets = block_offsets(&block_sizes);
95        let total_vars: usize = block_sizes.iter().sum();
96        let (gx, gy) = self.target_point();
97
98        let mut x_coeffs = Vec::with_capacity(total_vars);
99        let mut y_coeffs = Vec::with_capacity(total_vars);
100        for (&length, samples) in self.link_lengths().iter().zip(self.orientation_samples()) {
101            for &angle in samples {
102                x_coeffs.push(length * angle.cos());
103                y_coeffs.push(length * angle.sin());
104            }
105        }
106
107        // A violation contributes at least one full penalty unit. This bound
108        // exceeds the largest possible squared distance in exact arithmetic.
109        // A proportional gap avoids rounding B + 1 back to B at large scales;
110        // the floating-point objective still has finite-precision limitations.
111        let sum_abs_x: f64 = x_coeffs.iter().map(|coeff| coeff.abs()).sum();
112        let sum_abs_y: f64 = y_coeffs.iter().map(|coeff| coeff.abs()).sum();
113        let distance_bound = (sum_abs_x + gx.abs()).powi(2) + (sum_abs_y + gy.abs()).powi(2);
114        let penalty = 2.0 * (1.0 + distance_bound);
115        if !penalty.is_finite() {
116            return Err(crate::rules::ReductionError::non_finite_result::<
117                Self,
118                QUBO<f64>,
119            >("computing the inverse-kinematics penalty"));
120        }
121
122        let mut matrix = vec![vec![0.0; total_vars]; total_vars];
123        let mut add_upper = |i: usize, j: usize, value: f64| {
124            let (lo, hi) = if i <= j { (i, j) } else { (j, i) };
125            matrix[lo][hi] += value;
126        };
127
128        // Position objective: (X - g_x)^2 + (Y - g_y)^2, dropping the
129        // additive constant g_x^2 + g_y^2.
130        for (idx, (&x_coeff, &y_coeff)) in x_coeffs.iter().zip(&y_coeffs).enumerate() {
131            add_upper(
132                idx,
133                idx,
134                x_coeff * x_coeff - 2.0 * gx * x_coeff + y_coeff * y_coeff - 2.0 * gy * y_coeff,
135            );
136        }
137        for i in 0..total_vars {
138            for j in (i + 1)..total_vars {
139                add_upper(
140                    i,
141                    j,
142                    2.0 * (x_coeffs[i] * x_coeffs[j] + y_coeffs[i] * y_coeffs[j]),
143                );
144            }
145        }
146
147        // One-hot penalty per link: P * (sum_a y_{j,a} - 1)^2.
148        for (&start, &size) in block_offsets.iter().zip(&block_sizes) {
149            for a in 0..size {
150                add_upper(start + a, start + a, -penalty);
151            }
152            for a in 0..size {
153                for b in (a + 1)..size {
154                    add_upper(start + a, start + b, 2.0 * penalty);
155                }
156            }
157        }
158
159        // Forbidden-pair penalties between consecutive links.
160        for (junction, pairs) in self.allowed_pairs().iter().enumerate() {
161            let prev_size = block_sizes[junction];
162            let curr_size = block_sizes[junction + 1];
163            let prev_start = block_offsets[junction];
164            let curr_start = block_offsets[junction + 1];
165
166            let mut allowed = vec![vec![false; curr_size]; prev_size];
167            for &(a_prev, a_curr) in pairs {
168                allowed[a_prev][a_curr] = true;
169            }
170
171            for (a_prev, row) in allowed.iter().enumerate() {
172                for (a_curr, &is_allowed) in row.iter().enumerate() {
173                    if !is_allowed {
174                        add_upper(prev_start + a_prev, curr_start + a_curr, penalty);
175                    }
176                }
177            }
178        }
179
180        Ok(ReductionMinimumDiscretePlanarInverseKinematicsToQUBO {
181            target: QUBO::from_matrix(matrix).map_err(|message| {
182                crate::rules::ReductionError::construction::<
183                    MinimumDiscretePlanarInverseKinematics,
184                    QUBO<f64>,
185                >(message)
186            })?,
187            block_offsets,
188            block_sizes,
189            allowed_pairs: self.allowed_pairs().to_vec(),
190        })
191    }
192}
193
194#[cfg(feature = "example-db")]
195pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
196    use crate::export::SolutionPair;
197    use std::f64::consts::FRAC_PI_2;
198
199    vec![crate::example_db::specs::RuleExampleSpec {
200        id: "minimumdiscreteplanarinversekinematics_to_qubo",
201        build: || {
202            crate::example_db::specs::rule_example_with_witness::<_, QUBO<f64>>(
203                MinimumDiscretePlanarInverseKinematics::new(
204                    vec![2.0, 1.0],
205                    (2.0, 1.0),
206                    vec![vec![0.0, FRAC_PI_2], vec![0.0, FRAC_PI_2]],
207                    vec![vec![(0, 0), (0, 1), (1, 1)]],
208                )
209                .unwrap(),
210                SolutionPair {
211                    source_config: serde_json::json!(vec![0, 1]),
212                    target_config: serde_json::json!(vec![true, false, false, true]),
213                },
214            )
215        },
216    }]
217}
218
219#[cfg(test)]
220#[path = "../unit_tests/rules/minimumdiscreteplanarinversekinematics_qubo.rs"]
221mod tests;
problemreductions/rules/minimumdiscreteplanarinversekinematics_qubo.rs

problemreductions/rules/
minimumdiscreteplanarinversekinematics_qubo.rs
