diff --git a/docs/python-interface.md b/docs/python-interface.md index 679ef0ba..a332ea29 100644 --- a/docs/python-interface.md +++ b/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 size $\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]])`. | diff --git a/open-codegen/CHANGELOG.md b/open-codegen/CHANGELOG.md index 0a0a22f2..81a2df78 100644 --- a/open-codegen/CHANGELOG.md +++ b/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 diff --git a/open-codegen/VERSION b/open-codegen/VERSION index 7deb86fe..8adc70fd 100644 --- a/open-codegen/VERSION +++ b/open-codegen/VERSION @@ -1 +1 @@ -0.7.1 \ No newline at end of file +0.8.0 \ No newline at end of file diff --git a/open-codegen/opengen/constraints/__init__.py b/open-codegen/opengen/constraints/__init__.py index a0aa9477..4fbf168b 100644 --- a/open-codegen/opengen/constraints/__init__.py +++ b/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 * diff --git a/open-codegen/opengen/constraints/sphere2.py b/open-codegen/opengen/constraints/sphere2.py new file mode 100644 index 00000000..f8fca5f2 --- /dev/null +++ b/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 diff --git a/open-codegen/opengen/templates/optimizer.rs b/open-codegen/opengen/templates/optimizer.rs index b6019a61..912784db 100644 --- a/open-codegen/opengen/templates/optimizer.rs +++ b/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}}); diff --git a/open-codegen/test/test_constraints.py b/open-codegen/test/test_constraints.py index a2022513..4421cbd8 100644 --- a/open-codegen/test/test_constraints.py +++ b/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()