From 352417cb13a104115eba204582d1d9806231ffe2 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 11:46:47 +0100 Subject: [PATCH 1/9] first untested impl of Sphere2 --- src/constraints/sphere2.rs | 42 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 src/constraints/sphere2.rs diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs new file mode 100644 index 00000000..aeb3d72c --- /dev/null +++ b/src/constraints/sphere2.rs @@ -0,0 +1,42 @@ +use super::Constraint; + +#[derive(Copy, Clone)] +/// A Euclidean sphere, that is, a set given by $B_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert = r\\}$ +/// or a Euclidean sphere centered at a point $x_c$, that is, $B_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> { + fn project(&self, x: &mut [f64]) { + if let Some(center) = &self.center { + let mut norm_difference = 0.0; + x.iter().zip(center.iter()).for_each(|(a, b)| { + let diff_ = *a - *b; + norm_difference += diff_ * diff_ + }); + norm_difference = norm_difference.sqrt(); + x.iter_mut().zip(center.iter()).for_each(|(x, c)| { + *x = *c + (*x - *c) / norm_difference; + }); + } else { + let norm_x = crate::matrix_operations::norm2(x); + let norm_over_radius = norm_x / self.radius; + x.iter_mut().for_each(|x_| *x_ /= norm_over_radius); + } + } + + fn is_convex(&self) -> bool { + true + } +} From d227dfd93f2917bbb8cd35606ec4c5c7185679df Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 11:56:41 +0100 Subject: [PATCH 2/9] sphere2: test with no centre --- src/constraints/mod.rs | 2 ++ src/constraints/tests.rs | 14 ++++++++++++++ 2 files changed, 16 insertions(+) diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index 7939e288..ea8895da 100644 --- a/src/constraints/mod.rs +++ b/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 diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index eecb1a6c..4d770515 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -763,6 +763,20 @@ fn t_ball1_random_optimality_conditions_centered() { } } +#[test] +fn t_sphere2_no_center() { + let radius = 1.0; + let mut x_out = [1.0, 1.0]; + let mut x_in = [1.0, 1.0]; + 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(1.0, norm_out, 1e-10, 1e-12, "norm_out is not 1.0"); + unit_test_utils::assert_nearly_equal(1.0, norm_in, 1e-10, 1e-12, "norm_in is not 1.0"); +} + #[test] #[should_panic] fn t_ball1_alpha_negative() { From 757d72ef10480a523ac06141d3500dcdf1c3764c Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 12:09:46 +0100 Subject: [PATCH 3/9] sphere2: first tests --- src/constraints/tests.rs | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 4d770515..d85f9a9e 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -773,8 +773,28 @@ fn t_sphere2_no_center() { 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(1.0, norm_out, 1e-10, 1e-12, "norm_out is not 1.0"); - unit_test_utils::assert_nearly_equal(1.0, norm_in, 1e-10, 1e-12, "norm_in is not 1.0"); + 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_center() { + let radius = 1.0; + let center = [-3.0, 5.0]; + let mut x = [1.0, 1.0]; + let unit_sphere = Sphere2::new(Some(¢er), 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] From 125444d0343921b12c5629f4cae80972442e0d9c Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 14:13:54 +0100 Subject: [PATCH 4/9] sphere: fix bug scaling properly by radius more tests --- src/constraints/sphere2.rs | 8 ++++---- src/constraints/tests.rs | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs index aeb3d72c..d0110ea4 100644 --- a/src/constraints/sphere2.rs +++ b/src/constraints/sphere2.rs @@ -21,13 +21,13 @@ impl<'a> Constraint for Sphere2<'a> { fn project(&self, x: &mut [f64]) { if let Some(center) = &self.center { let mut norm_difference = 0.0; - x.iter().zip(center.iter()).for_each(|(a, b)| { - let diff_ = *a - *b; + x.iter().zip(center.iter()).for_each(|(xi, ci)| { + let diff_ = *xi - *ci; norm_difference += diff_ * diff_ }); norm_difference = norm_difference.sqrt(); x.iter_mut().zip(center.iter()).for_each(|(x, c)| { - *x = *c + (*x - *c) / norm_difference; + *x = *c + self.radius * (*x - *c) / norm_difference; }); } else { let norm_x = crate::matrix_operations::norm2(x); @@ -37,6 +37,6 @@ impl<'a> Constraint for Sphere2<'a> { } fn is_convex(&self) -> bool { - true + false } } diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index d85f9a9e..09a660da 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -765,9 +765,9 @@ fn t_ball1_random_optimality_conditions_centered() { #[test] fn t_sphere2_no_center() { - let radius = 1.0; + let radius = 0.9; let mut x_out = [1.0, 1.0]; - let mut x_in = [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); @@ -779,7 +779,7 @@ fn t_sphere2_no_center() { #[test] fn t_sphere2_center() { - let radius = 1.0; + let radius = 1.3; let center = [-3.0, 5.0]; let mut x = [1.0, 1.0]; let unit_sphere = Sphere2::new(Some(¢er), radius); From 5d8b7ab4b6a675849b0a11e0d9a9d3bd2aca512e Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 15:47:25 +0100 Subject: [PATCH 5/9] update Chargo.toml and CHANGELOG --- CHANGELOG.md | 10 ++++++++++ Cargo.toml | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index aaafe05b..503b3e19 100644 --- a/CHANGELOG.md +++ b/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.7.7] - 2023-10-27 + +### Added + +- New constraint: sphere of Euclidean norm + @@ -247,6 +256,7 @@ This is a breaking API change. --------------------- --> +[v0.7.8]: https://github.com/alphaville/optimization-engine/compare/v0.7.7...v0.7.8 [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 diff --git a/Cargo.toml b/Cargo.toml index 77f18b03..231a14ac 100644 --- a/Cargo.toml +++ b/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.7.8" edition = "2018" From a0646eeb525b1d34f887289f3bca449fb3756d24 Mon Sep 17 00:00:00 2001 From: Ruairi Moran Date: Fri, 27 Oct 2023 10:17:37 -0500 Subject: [PATCH 6/9] typo --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 503b3e19..a1bd7ed9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,7 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo -## [v0.7.7] - 2023-10-27 +## [v0.7.8] - 2023-10-27 ### Added From cdcf995456b18cc3497fad58280c032367864b87 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 17:25:55 +0100 Subject: [PATCH 7/9] sphere: handling projection of centre --- src/constraints/sphere2.rs | 33 +++++++++++++++++++++++++-------- src/constraints/tests.rs | 30 ++++++++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 8 deletions(-) diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs index d0110ea4..060b5d31 100644 --- a/src/constraints/sphere2.rs +++ b/src/constraints/sphere2.rs @@ -18,21 +18,38 @@ impl<'a> Sphere2<'a> { } impl<'a> Constraint for Sphere2<'a> { + /// Projection onto the set, that is, + /// + /// $$ + /// \Pi_C(v) = \mathrm{argmin}_{z\in C}\Vert{}z-v{}\Vert + /// $$ + /// + /// + /// + /// ## 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 mut norm_difference = 0.0; - x.iter().zip(center.iter()).for_each(|(xi, ci)| { - let diff_ = *xi - *ci; - norm_difference += diff_ * diff_ - }); - norm_difference = norm_difference.sqrt(); + let norm_difference = crate::matrix_operations::norm2_squared_diff(x, center).sqrt(); + if norm_difference <= epsilon { + x.copy_from_slice(¢er); + 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); - let norm_over_radius = norm_x / self.radius; - x.iter_mut().for_each(|x_| *x_ /= norm_over_radius); + 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); } } diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 09a660da..f73a3d7b 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -777,6 +777,16 @@ fn t_sphere2_no_center() { 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; @@ -797,6 +807,26 @@ fn t_sphere2_center() { 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(¢er), 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() { From d0834161ea79803c4178c6e96cd1988fb89ce217 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 17:29:07 +0100 Subject: [PATCH 8/9] Update version number again According to semver we bump the MINOR version number whenever we add functionality in a backwards-comp. manner New version: 0.8.0 --- CHANGELOG.md | 6 +++--- Cargo.toml | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a1bd7ed9..aba8337c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,9 +9,9 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo -## [v0.7.8] - 2023-10-27 +## [v0.8.1] - 2023-10-27 ### Added @@ -256,7 +256,7 @@ This is a breaking API change. --------------------- --> -[v0.7.8]: https://github.com/alphaville/optimization-engine/compare/v0.7.7...v0.7.8 +[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 diff --git a/Cargo.toml b/Cargo.toml index 231a14ac..13adc99f 100644 --- a/Cargo.toml +++ b/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.8" +version = "0.8.0" edition = "2018" From 493c5eb9a4dbfb6bbdf678a51802d3746dcafee3 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Fri, 27 Oct 2023 17:42:20 +0100 Subject: [PATCH 9/9] sphere2: API documentation --- src/constraints/sphere2.rs | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs index 060b5d31..2befc82b 100644 --- a/src/constraints/sphere2.rs +++ b/src/constraints/sphere2.rs @@ -1,8 +1,8 @@ use super::Constraint; #[derive(Copy, Clone)] -/// A Euclidean sphere, that is, a set given by $B_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert = r\\}$ -/// or a Euclidean sphere centered at a point $x_c$, that is, $B_2^{x_c, r} = \\{x \in \mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert = r\\}$ +/// 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, @@ -18,13 +18,16 @@ impl<'a> Sphere2<'a> { } impl<'a> Constraint for Sphere2<'a> { - /// Projection onto the set, that is, + /// Projection onto the sphere, $S_{r, c}$ with radius $r$ and center $c$. + /// If $x\neq c$, the projection is uniquely defined by /// /// $$ - /// \Pi_C(v) = \mathrm{argmin}_{z\in C}\Vert{}z-v{}\Vert + /// 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 /// @@ -53,6 +56,8 @@ impl<'a> Constraint for Sphere2<'a> { } } + /// Returns false (the sphere is not a convex set) + /// fn is_convex(&self) -> bool { false }