Merge pull request #329 from alphaville/feature/sphere.rs

Spherical constraint

CHANGELOG.md

@@ -8,6 +8,15 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 Note: This is the main Changelog file for the Rust solver. The Changelog file for the Python interface (`opengen`) can be found in [/open-codegen/CHANGELOG.md](open-codegen/CHANGELOG.md)
 
 
+<!-- ---------------------
+      v0.8.0
+     --------------------- -->
+## [v0.8.1] - 2023-10-27
+
+### Added
+
+- New constraint: sphere of Euclidean norm
+
 <!-- ---------------------
       v0.7.7
      --------------------- -->
@@ -247,6 +256,7 @@ This is a breaking API change.
      --------------------- -->
 
 <!-- Releases -->
+[v0.7.8]: https://github.com/alphaville/optimization-engine/compare/v0.7.7...v0.8.0
 [v0.7.7]: https://github.com/alphaville/optimization-engine/compare/v0.7.6...v0.7.7 
 [v0.7.6]: https://github.com/alphaville/optimization-engine/compare/v0.7.5...v0.7.6 
 [v0.7.5]: https://github.com/alphaville/optimization-engine/compare/v0.7.4...v0.7.5 

Cargo.toml

@@ -42,7 +42,7 @@ homepage = "https://alphaville.github.io/optimization-engine/"
 repository = "https://github.com/alphaville/optimization-engine"
 
 # Version of this crate (SemVer)
-version = "0.7.7"
+version = "0.8.0"
 
 edition = "2018"
 

src/constraints/mod.rs

@@ -19,6 +19,7 @@ mod no_constraints;
 mod rectangle;
 mod simplex;
 mod soc;
+mod sphere2;
 mod zero;
 
 pub use ball1::Ball1;
@@ -32,6 +33,7 @@ pub use no_constraints::NoConstraints;
 pub use rectangle::Rectangle;
 pub use simplex::Simplex;
 pub use soc::SecondOrderCone;
+pub use sphere2::Sphere2;
 pub use zero::Zero;
 
 /// A set which can be used as a constraint

src/constraints/sphere2.rs

@@ -0,0 +1,64 @@
+use super::Constraint;
+
+#[derive(Copy, Clone)]
+/// A Euclidean sphere, that is, a set given by $S_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert = r\\}$
+/// 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\\}$
+pub struct Sphere2<'a> {
+    center: Option<&'a [f64]>,
+    radius: f64,
+}
+
+impl<'a> Sphere2<'a> {
+    /// Construct a new Euclidean sphere with given center and radius
+    /// If no `center` is given, then it is assumed to be in the origin
+    pub fn new(center: Option<&'a [f64]>, radius: f64) -> Self {
+        assert!(radius > 0.0);
+        Sphere2 { center, radius }
+    }
+}
+
+impl<'a> Constraint for Sphere2<'a> {
+    /// Projection onto the sphere, $S_{r, c}$ with radius $r$ and center $c$.
+    /// If $x\neq c$, the projection is uniquely defined by
+    ///
+    /// $$
+    /// P_{S_{r, c}}(x) = c + r\frac{x-c}{\Vert{}x-c\Vert_2},
+    /// $$
+    ///
+    /// but for $x=c$, the projection is multi-valued. In particular, let
+    /// $y = P_{S_{r, c}}(c)$. Then $y_1 = c_1 + r$ and $y_i = c_i$ for
+    /// $i=2,\ldots, n$.
+    ///
+    /// ## Arguments
+    ///
+    /// - `x`: The given vector $x$ is updated with the projection on the set
+    ///
+    fn project(&self, x: &mut [f64]) {
+        let epsilon = 1e-12;
+        if let Some(center) = &self.center {
+            let norm_difference = crate::matrix_operations::norm2_squared_diff(x, center).sqrt();
+            if norm_difference <= epsilon {
+                x.copy_from_slice(&center);
+                x[0] += self.radius;
+                return;
+            }
+            x.iter_mut().zip(center.iter()).for_each(|(x, c)| {
+                *x = *c + self.radius * (*x - *c) / norm_difference;
+            });
+        } else {
+            let norm_x = crate::matrix_operations::norm2(x);
+            if norm_x <= epsilon {
+                x[0] += self.radius;
+                return;
+            }
+            let norm_over_radius = self.radius / norm_x;
+            x.iter_mut().for_each(|x_| *x_ *= norm_over_radius);
+        }
+    }
+
+    /// Returns false (the sphere is not a convex set)
+    ///
+    fn is_convex(&self) -> bool {
+        false
+    }
+}

src/constraints/tests.rs

@@ -763,6 +763,70 @@ fn t_ball1_random_optimality_conditions_centered() {
     }
 }
 
+#[test]
+fn t_sphere2_no_center() {
+    let radius = 0.9;
+    let mut x_out = [1.0, 1.0];
+    let mut x_in = [-0.3, -0.2];
+    let unit_sphere = Sphere2::new(None, radius);
+    unit_sphere.project(&mut x_out);
+    unit_sphere.project(&mut x_in);
+    let norm_out = crate::matrix_operations::norm2(&x_out);
+    let norm_in = crate::matrix_operations::norm2(&x_in);
+    unit_test_utils::assert_nearly_equal(radius, norm_out, 1e-10, 1e-12, "norm_out is not 1.0");
+    unit_test_utils::assert_nearly_equal(radius, norm_in, 1e-10, 1e-12, "norm_in is not 1.0");
+}
+
+#[test]
+fn t_sphere2_no_center_projection_of_zero() {
+    let radius = 0.9;
+    let mut x = [0.0, 0.0];
+    let unit_sphere = Sphere2::new(None, radius);
+    unit_sphere.project(&mut x);
+    let norm_result = crate::matrix_operations::norm2(&x);
+    unit_test_utils::assert_nearly_equal(radius, norm_result, 1e-10, 1e-12, "norm_out is not 1.0");
+}
+
+#[test]
+fn t_sphere2_center() {
+    let radius = 1.3;
+    let center = [-3.0, 5.0];
+    let mut x = [1.0, 1.0];
+    let unit_sphere = Sphere2::new(Some(&center), radius);
+
+    unit_sphere.project(&mut x);
+    let mut x_minus_c = [0.0; 2];
+    x.iter()
+        .zip(center.iter())
+        .zip(x_minus_c.iter_mut())
+        .for_each(|((a, b), c)| {
+            *c = a - b;
+        });
+
+    let norm_out = crate::matrix_operations::norm2(&x_minus_c);
+    unit_test_utils::assert_nearly_equal(radius, norm_out, 1e-10, 1e-12, "norm_out is not 1.0");
+}
+
+#[test]
+fn t_sphere2_center_projection_of_center() {
+    let radius = 1.3;
+    let center = [-3.0, 5.0];
+    let mut x = [-3.0, 5.0];
+    let unit_sphere = Sphere2::new(Some(&center), radius);
+
+    unit_sphere.project(&mut x);
+    let mut x_minus_c = [0.0; 2];
+    x.iter()
+        .zip(center.iter())
+        .zip(x_minus_c.iter_mut())
+        .for_each(|((a, b), c)| {
+            *c = a - b;
+        });
+
+    let norm_out = crate::matrix_operations::norm2(&x_minus_c);
+    unit_test_utils::assert_nearly_equal(radius, norm_out, 1e-10, 1e-12, "norm_out is not 1.0");
+}
+
 #[test]
 #[should_panic]
 fn t_ball1_alpha_negative() {