optimization_engine/constraints/
sphere2.rs

1use super::Constraint;
2
3#[derive(Copy, Clone)]
4/// A Euclidean sphere, that is, a set given by $S_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert = r\\}$
5/// or a Euclidean sphere centered at a point $x_c$, that is, $S_2^{x_c, r} = \\{x \in \mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert = r\\}$
6pub struct Sphere2<'a> {
7    center: Option<&'a [f64]>,
8    radius: f64,
9}
10
11impl<'a> Sphere2<'a> {
12    /// Construct a new Euclidean sphere with given center and radius
13    /// If no `center` is given, then it is assumed to be in the origin
14    pub fn new(center: Option<&'a [f64]>, radius: f64) -> Self {
15        assert!(radius > 0.0);
16        Sphere2 { center, radius }
17    }
18}
19
20impl<'a> Constraint for Sphere2<'a> {
21    /// Projection onto the sphere, $S_{r, c}$ with radius $r$ and center $c$.
22    /// If $x\neq c$, the projection is uniquely defined by
23    ///
24    /// $$
25    /// P_{S_{r, c}}(x) = c + r\frac{x-c}{\Vert{}x-c\Vert_2},
26    /// $$
27    ///
28    /// but for $x=c$, the projection is multi-valued. In particular, let
29    /// $y = P_{S_{r, c}}(c)$. Then $y_1 = c_1 + r$ and $y_i = c_i$ for
30    /// $i=2,\ldots, n$.
31    ///
32    /// ## Arguments
33    ///
34    /// - `x`: The given vector $x$ is updated with the projection on the set
35    ///
36    /// ## Panics
37    ///
38    /// Panics if `x` is empty or, when a center is provided, if `x` and
39    /// `center` have incompatible dimensions.
40    ///
41    fn project(&self, x: &mut [f64]) {
42        let epsilon = 1e-12;
43        assert!(!x.is_empty(), "x must be nonempty");
44        if let Some(center) = &self.center {
45            assert_eq!(
46                x.len(),
47                center.len(),
48                "x and center have incompatible dimensions"
49            );
50            let norm_difference = crate::matrix_operations::norm2_squared_diff(x, center).sqrt();
51            if norm_difference <= epsilon {
52                x.copy_from_slice(center);
53                x[0] += self.radius;
54                return;
55            }
56            x.iter_mut().zip(center.iter()).for_each(|(x, c)| {
57                *x = *c + self.radius * (*x - *c) / norm_difference;
58            });
59        } else {
60            let norm_x = crate::matrix_operations::norm2(x);
61            if norm_x <= epsilon {
62                x[0] += self.radius;
63                return;
64            }
65            let norm_over_radius = self.radius / norm_x;
66            x.iter_mut().for_each(|x_| *x_ *= norm_over_radius);
67        }
68    }
69
70    /// Returns false (the sphere is not a convex set)
71    ///
72    fn is_convex(&self) -> bool {
73        false
74    }
75}
optimization_engine/constraints/sphere2.rs

optimization_engine/constraints/
sphere2.rs
