Merge pull request #193 from alphaville/feature/hyperplane

Projection on hyperplanes

CHANGELOG.md

@@ -9,17 +9,22 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo
 
 
 <!-- ---------------------
-      Unreleased; 
-      Next version: 0.7.1
+      Unreleased
      --------------------- -->
 
-## [Unreleased]
 
-### Added
+<!-- ---------------------
+      v0.7.1-alpha.1
+     --------------------- -->
 
-* Introduced `Halfspace` (implemented and tested)
+## [v0.7.1-alpha.1]
 
+### Added
 
+* Introduced `Halfspace` (implemented and tested)
+* Introduced `Hyperplane` (implemented and tested)
+* New types: `FunctionCallResult`, `MappingType` and `JacobianMappingType`
+* Various clippy-related code improvements
 
 
 
@@ -175,6 +180,7 @@ This is a breaking API change.
 
 <!-- Releases -->
 [Unreleased]: https://github.com/alphaville/optimization-engine/compare/v0.7.0...master
+[v0.7.1]: https://github.com/alphaville/optimization-engine/compare/v0.7.0...v0.7.1
 [v0.7.0]: https://github.com/alphaville/optimization-engine/compare/v0.6.2...v0.7.0
 [v0.6.2]: https://github.com/alphaville/optimization-engine/compare/v0.6.1-alpha.2...v0.6.2
 [v0.6.1-alpha.2]: https://github.com/alphaville/optimization-engine/compare/v0.5.0...v0.6.1-alpha.2

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.0"
+version = "0.7.1-alpha.1"
 
 edition = "2018"
 

src/constraints/halfspace.rs

@@ -56,13 +56,36 @@ impl<'a> Halfspace<'a> {
 }
 
 impl<'a> Constraint for Halfspace<'a> {
+    /// Projects on halfspace using the following formula:
+    ///
+    /// $$\begin{aligned}
+    /// \mathrm{proj}_{H}(x) = \begin{cases}
+    /// x,& \text{ if } \langle c, x\rangle \leq b
+    /// \\\\
+    /// x - \frac{\langle c, x\rangle - b}
+    ///          {\\|c\\|}c,& \text{else}
+    /// \end{cases}
+    /// \end{aligned}$$
+    ///
+    /// where $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle \leq b\\}$
+    ///
+    /// # Arguments
+    ///
+    /// - `x`: (in) vector to be projected on the current instance of a halfspace,
+    ///    (out) projection on the second-order cone
+    ///
+    /// # Panics
+    ///
+    /// This method panics if the length of `x` is not equal to the dimension
+    /// of the halfspace.
+    ///
     fn project(&self, x: &mut [f64]) {
         let inner_product = matrix_operations::inner_product(x, self.normal_vector);
         if inner_product > self.offset {
             let factor = (inner_product - self.offset) / self.normal_vector_squared_norm;
             x.iter_mut()
                 .zip(self.normal_vector.iter())
-                .for_each(|(x, nrm_vct)| *x -= factor * nrm_vct);
+                .for_each(|(x, normal_vector_i)| *x -= factor * normal_vector_i);
         }
     }
 

src/constraints/hyperplane.rs

@@ -0,0 +1,95 @@
+use super::Constraint;
+use crate::matrix_operations;
+
+#[derive(Clone)]
+/// A hyperplane is a set given by $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle = b\\}$.
+pub struct Hyperplane<'a> {
+    /// normal vector
+    normal_vector: &'a [f64],
+    /// offset
+    offset: f64,
+    /// squared Euclidean norm of the normal vector (computed once upon construction)
+    normal_vector_squared_norm: f64,
+}
+
+impl<'a> Hyperplane<'a> {
+    /// A hyperplane is a set given by $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle = b\\}$,
+    /// where $c$ is the normal vector of the hyperplane and $b$ is an offset.
+    ///
+    /// This method constructs a new instance of `Hyperplane` with a given normal
+    /// vector and bias
+    ///
+    /// # Arguments
+    ///
+    /// - `normal_vector`: the normal vector, $c$, as a slice
+    /// - `offset`: the offset parameter, $b$
+    ///
+    /// # Returns
+    ///
+    /// New instance of `Hyperplane`
+    ///
+    /// # Panics
+    ///
+    /// Does not panic. Note: it does not panic if you provide an empty slice as `normal_vector`,
+    /// but you should avoid doing that.
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use optimization_engine::constraints::{Constraint, Hyperplane};
+    ///
+    /// let normal_vector = [1., 2.];
+    /// let offset = 1.0;
+    /// let hyperplane = Hyperplane::new(&normal_vector, offset);
+    /// let mut x = [-1., 3.];
+    /// hyperplane.project(&mut x);
+    /// ```
+    ///
+    pub fn new(normal_vector: &'a [f64], offset: f64) -> Self {
+        let normal_vector_squared_norm = matrix_operations::norm2_squared(normal_vector);
+        Hyperplane {
+            normal_vector,
+            offset,
+            normal_vector_squared_norm,
+        }
+    }
+}
+
+impl<'a> Constraint for Hyperplane<'a> {
+    /// Projects on the hyperplane using the formula:
+    ///
+    /// $$\begin{aligned}
+    /// \mathrm{proj}_{H}(x) =
+    /// x - \frac{\langle c, x\rangle - b}
+    ///          {\\|c\\|}c.
+    /// \end{aligned}$$
+    ///
+    /// where $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle = b\\}$
+    ///
+    /// # Arguments
+    ///
+    /// - `x`: (in) vector to be projected on the current instance of a hyperplane,
+    ///    (out) projection on the second-order cone
+    ///
+    /// # Panics
+    ///
+    /// This method panics if the length of `x` is not equal to the dimension
+    /// of the hyperplane.
+    ///
+    fn project(&self, x: &mut [f64]) {
+        let inner_product = matrix_operations::inner_product(x, self.normal_vector);
+        let factor = (inner_product - self.offset) / self.normal_vector_squared_norm;
+        x.iter_mut()
+            .zip(self.normal_vector.iter())
+            .for_each(|(x, nrm_vct)| *x -= factor * nrm_vct);
+    }
+
+    /// Hyperplanes are convex sets
+    ///
+    /// # Returns
+    ///
+    /// Returns `true`
+    fn is_convex(&self) -> bool {
+        true
+    }
+}

src/constraints/mod.rs

@@ -13,6 +13,7 @@ mod ballinf;
 mod cartesian_product;
 mod finite;
 mod halfspace;
+mod hyperplane;
 mod no_constraints;
 mod rectangle;
 mod soc;
@@ -23,6 +24,7 @@ pub use ballinf::BallInf;
 pub use cartesian_product::CartesianProduct;
 pub use finite::FiniteSet;
 pub use halfspace::Halfspace;
+pub use hyperplane::Hyperplane;
 pub use no_constraints::NoConstraints;
 pub use rectangle::Rectangle;
 pub use soc::SecondOrderCone;

src/constraints/tests.rs

@@ -15,6 +15,27 @@ fn t_zero_set() {
     );
 }
 
+#[test]
+fn t_hyperplane() {
+    let normal_vector = [1.0, 2.0, 3.0];
+    let offset = 1.0;
+    let hyperplane = Hyperplane::new(&normal_vector, offset);
+    let mut x = [-1., 3., 5.];
+    let x_proj_expected = [
+        -2.357_142_857_142_857,
+        0.285_714_285_714_286,
+        0.928_571_428_571_429,
+    ];
+    hyperplane.project(&mut x);
+    unit_test_utils::assert_nearly_equal_array(
+        &x,
+        &x_proj_expected,
+        1e-8,
+        1e-14,
+        "halfspace projection failed",
+    );
+}
+
 #[test]
 fn t_halfspace_project_inside() {
     let normal_vector = [1., 2.];
@@ -249,7 +270,7 @@ fn t_no_constraints() {
 
     whole_space.project(&mut x);
 
-    assert_eq!([1., 2., 3.], x);
+    unit_test_utils::assert_nearly_equal_array(&[1., 2., 3.], &x, 1e-10, 1e-15, "x is wrong");
 }
 
 #[test]
@@ -577,3 +598,11 @@ fn t_is_convex_cartesian_product() {
     let cartesian_product = cartesian_product.add_constraint(10, finite_noncvx);
     assert!(!cartesian_product.is_convex());
 }
+
+#[test]
+fn t_hyperplane_is_convex() {
+    let normal_vector = [1.0, 2.0, 3.0];
+    let offset = 1.0;
+    let hyperplane = Hyperplane::new(&normal_vector, offset);
+    assert!(hyperplane.is_convex());
+}

src/core/fbs/fbs_optimizer.rs

@@ -3,10 +3,10 @@
 use crate::{
     constraints,
     core::{
-        fbs::fbs_engine::FBSEngine, fbs::FBSCache, AlgorithmEngine, ExitStatus, FunctionCallResult,
-        Optimizer, Problem, SolverStatus,
+        fbs::fbs_engine::FBSEngine, fbs::FBSCache, AlgorithmEngine, ExitStatus, Optimizer, Problem,
+        SolverStatus,
     },
-    matrix_operations, SolverError,
+    matrix_operations, FunctionCallResult, SolverError,
 };
 use std::time;
 
@@ -120,7 +120,7 @@ where
         }
 
         // cost at the solution [propagate error upstream]
-        let mut cost_value = 0.0;
+        let mut cost_value: f64 = 0.0;
         (self.fbs_engine.problem.cost)(u, &mut cost_value)?;
 
         if !matrix_operations::is_finite(&u) || !cost_value.is_finite() {

src/core/fbs/tests.rs

@@ -64,7 +64,6 @@ fn t_fbs_step_no_constraints() {
         let mut u = [1.0, 3.0];
 
         assert!(fbs_engine.step(&mut u).unwrap());
-        assert_eq!([0.5, 2.4], u);
         unit_test_utils::assert_nearly_equal_array(&[0.5, 2.4], &u, 1e-10, 1e-14, "u");
     }
     unit_test_utils::assert_nearly_equal_array(