Code generation support for Sphere2

docs/python-interface.md

@@ -80,6 +80,7 @@ following types of constraints:
 | `Ball2`            | Euclidean ball: `Ball2(None, r)` creates a Euclidean ball of radius `r` centered at the origin, and `Ball2(xc, r)` is a ball centered at point `xc` (list/np.array) |
 | `BallInf`          | Ball of infinity norm:`BallInf(None, r)` creates an infinity-norm ball of radius `r` centered at the origin, and `BallInf(xc, r)` is an infinity ball centered at point `xc` (list/np.array) |
 | `Ball1`            | L1 ball: `Ball(None, r)` creates an ell1-ball of radius `r` centered at the origin, and `BallInf(xc, r)` is an ell1-ball centered at point `xc` (list/np.array)|
+| `Sphere2`            | Euclidean sphere: `Sphere2(None, r)` creates a Euclidean sphere of radius `r` centered at the origin, and `Sphere2(xc, r)` is a sphere centered at point `xc` (list/np.array) |
 | `Simplex`            | A simplex of <em>size</em> $\alpha$ is a set of the form $\Delta_\alpha = \\{x \in \mathbb{R}^n {}:{} x_i \geq 0, \sum_i x_i = \alpha\\}$. Create one with `Simplex(alpha)`. Projections are computed using Condat's [fast projection method](https://link.springer.com/article/10.1007/s10107-015-0946-6).  |
 | `Halfspace`        | A halfspace is a set of the form $\\{u \in \mathbb{R}^{n_u} {}:{} \langle c, u\rangle \leq b \\}$, for a vector $c$ and a scalar $b$. The syntax is straightforwarrd: `Halfspace(c, b)`. |
 | `FiniteSet`        | Finite set, $\\{u^{(1)},\ldots,u^{(m)}\\}$; the set of point is provided as a list of lists, for example, `FiniteSet([[1,2],[2,3],[4,5]])`. The commonly used set of binary numbers, $\\{0, 1\\}$, is created with `FiniteSet([[0], [1]])`. |

open-codegen/CHANGELOG.md

@@ -7,7 +7,11 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 Note: This is the Changelog file of `opengen` - the Python interface of OpEn
 
+## [0.8.0] - 2024-03-20
 
+### Added
+
+* Code generation support for Sphere2
 
 ## [0.7.1] - 2022-02-14
 
@@ -181,6 +185,7 @@ Note: This is the Changelog file of `opengen` - the Python interface of OpEn
 * Project-specific `tcp_iface` TCP interface
 * Fixed `lbfgs` typo
 
+[0.8.0]: https://github.com/alphaville/optimization-engine/compare/opengen-0.7.1...opengen-0.8.0
 [0.7.1]: https://github.com/alphaville/optimization-engine/compare/opengen-0.7.0...opengen-0.7.1
 [0.7.0]: https://github.com/alphaville/optimization-engine/compare/opengen-0.6.13...opengen-0.7.0
 [0.6.13]: https://github.com/alphaville/optimization-engine/compare/opengen-0.6.12...opengen-0.6.13

open-codegen/VERSION

@@ -1 +1 @@
-0.7.1
+0.8.0

open-codegen/opengen/constraints/__init__.py

@@ -1,5 +1,6 @@
 from .ball1 import *
 from .ball2 import *
+from .sphere2 import *
 from .rectangle import *
 from .constraint import *
 from .ball_inf import *

open-codegen/opengen/constraints/sphere2.py

@@ -0,0 +1,77 @@
+import casadi.casadi as cs
+import numpy as np
+from .constraint import Constraint
+import opengen.functions as fn
+
+
+class Sphere2(Constraint):
+    """A Euclidean sphere constraint
+
+    A constraint of the form :math:`\|u-u_0\| = r`, where :math:`u_0` is the center
+    of the ball and `r` is its radius
+
+    """
+
+    def __init__(self, center=None, radius: float = 1.0):
+        """Constructor for a Euclidean sphere constraint
+
+        :param center: center of the sphere; if this is equal to Null, the
+            sphere is centered at the origin
+
+        :param radius: radius of the sphere
+
+        :return: New instance of Sphere2 with given center and radius
+        """
+        if radius <= 0:
+            raise Exception("The radius must be a positive number")
+
+        if center is not None and not isinstance(center, (list, np.ndarray)):
+            raise Exception("center is neither None nor a list nor np.ndarray")
+
+        self.__center = None if center is None else np.array(
+            [float(i) for i in center])
+        self.__radius = float(radius)
+
+    @property
+    def center(self):
+        """Returns the center of the ball"""
+        return self.__center
+
+    @property
+    def radius(self):
+        """Returns the radius of the sphere"""
+        return self.__radius
+
+    def distance_squared(self, u):
+        """Computes the squared distance between a given point `u` and this sphere
+
+            :param u: given point; can be a list of float, a numpy
+                n-dim array (`ndarray`) or a CasADi SX/MX symbol
+
+            :return: distance from set as a float or a CasADi symbol
+        """
+        if fn.is_symbolic(u):
+            # Case I: `u` is a CasADi SX symbol
+            v = u if self.__center is None else u - self.__center
+        elif (isinstance(u, list) and all(isinstance(x, (int, float)) for x in u))\
+                or isinstance(u, np.ndarray):
+            if self.__center is None:
+                v = u
+            else:
+                # Note: self.__center is np.ndarray (`u` might be a list)
+                z = self.__center.reshape(len(u))
+                u = np.array(u).reshape(len(u))
+                v = np.subtract(u, z)
+        else:
+            raise Exception("u is of invalid type")
+
+        return (self.__radius - fn.norm2(v))**2
+
+    def project(self, u):
+        raise NotImplementedError()
+
+    def is_convex(self):
+        return False
+
+    def is_compact(self):
+        return True

open-codegen/opengen/templates/optimizer.rs

@@ -65,7 +65,7 @@ pub const {{meta.optimizer_name|upper}}_N2: usize = {{problem.dim_constraints_pe
 
 // ---Parameters of the constraints----------------------------------------------------------------------
 
-{% if 'Ball1' == problem.constraints.__class__.__name__ or 'Ball2' == problem.constraints.__class__.__name__ or 'BallInf' == problem.constraints.__class__.__name__ -%}
+{% if 'Ball1' == problem.constraints.__class__.__name__ or 'Ball2' == problem.constraints.__class__.__name__ or 'BallInf' == problem.constraints.__class__.__name__ or 'Sphere2' == problem.constraints.__class__.__name__ -%}
 /// Constraints: Centre of Ball
 const CONSTRAINTS_BALL_XC: Option<&[f64]> = {% if problem.constraints.center is not none %}Some(&[{{problem.constraints.center | join(', ')}}]){% else %}None{% endif %};
 
@@ -140,6 +140,9 @@ fn make_constraints() -> impl Constraint {
     {% elif 'Ball1' == problem.constraints.__class__.__name__ -%}
     // - Ball1:
     Ball1::new(CONSTRAINTS_BALL_XC, CONSTRAINTS_BALL_RADIUS)
+    {% elif 'Sphere2' == problem.constraints.__class__.__name__ -%}
+    // - Sphere2:
+    Sphere2::new(CONSTRAINTS_BALL_XC, CONSTRAINTS_BALL_RADIUS)
     {% elif 'Simplex' == problem.constraints.__class__.__name__ -%}
     // - Simplex:
     let alpha_simplex : f64 = {{problem.constraints.alpha}};
@@ -184,6 +187,11 @@ fn make_constraints() -> impl Constraint {
         let center_{{loop.index}}: Option<&[f64]> = {% if set_i.center is not none %}Some(&[{{set_i.center | join(', ')}}]){% else %}None{% endif %};
         let set_{{loop.index}} = Ball1::new(center_{{loop.index}}, radius_{{loop.index}});
         let bounds = bounds.add_constraint(idx_{{loop.index}}, set_{{loop.index}});
+        {% elif 'Sphere2' == set_i.__class__.__name__ -%}
+        let radius_{{loop.index}} = {{set_i.radius}};
+        let center_{{loop.index}}: Option<&[f64]> = {% if set_i.center is not none %}Some(&[{{set_i.center | join(', ')}}]){% else %}None{% endif %};
+        let set_{{loop.index}} = Sphere2::new(center_{{loop.index}}, radius_{{loop.index}});
+        let bounds = bounds.add_constraint(idx_{{loop.index}}, set_{{loop.index}});
         {% elif 'Simplex' == set_i.__class__.__name__ -%}
         let alpha_{{loop.index}} = {{set_i.alpha}};        
         let set_{{loop.index}} = Simplex::new(alpha_{{loop.index}});

open-codegen/test/test_constraints.py

@@ -183,7 +183,8 @@ def test_rectangle_noncompact(self):
         self.assertFalse(rect.is_compact())
 
     def test_rectangle_is_orthant(self):
-        rect = og.constraints.Rectangle([0, float('-inf')], [float('inf'), 0.0])
+        rect = og.constraints.Rectangle(
+            [0, float('-inf')], [float('inf'), 0.0])
         self.assertTrue(rect.is_orthant())
         rect = og.constraints.Rectangle([0, 0], [float('inf'), float('inf')])
         self.assertTrue(rect.is_orthant())
@@ -492,6 +493,26 @@ def test_ball1_project_random_points_center(self):
                 self.assertLessEqual(
                     np.dot(e-x_star, x-x_star), 1e-10, "Ball1 optimality conditions failed (2)")
 
+    # -----------------------------------------------------------------------
+    # Sphere2
+    # -----------------------------------------------------------------------
+
+    def test_sphere2_sq_distance(self):
+        sphere = og.constraints.Sphere2(center=[1, 1, 1, 1], radius=0.5)
+        self.assertFalse(sphere.is_convex())
+        u = [1, 1, 0, 0]
+        dist = sphere.distance_squared(u)
+        self.assertAlmostEqual(0.835786437626905, dist, places=12)
+
+    def test_sphere2_sq_distance_symbolic(self):
+        sphere = og.constraints.Sphere2(center=[1, 1, 1, 1], radius=0.5)
+        self.assertFalse(sphere.is_convex())
+        u = cs.SX.sym("u", 4)
+        dist = sphere.distance_squared(u)
+        fun = cs.Function('fun', [u], [dist])
+        z = fun([1, 1, 0, 0])
+        self.assertAlmostEqual(0.835786437626905, z, places=12)
+
 
 if __name__ == '__main__':
     unittest.main()