diff --git a/CHANGELOG.md b/CHANGELOG.md index aaafe05b..aba8337c 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.8.1] - 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.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 77f18b03..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.7" +version = "0.8.0" edition = "2018" 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/sphere2.rs b/src/constraints/sphere2.rs new file mode 100644 index 00000000..2befc82b --- /dev/null +++ b/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(¢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); + 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 + } +} diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index eecb1a6c..f73a3d7b 100644 --- a/src/constraints/tests.rs +++ b/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(¢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] +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() {