Merge pull request #344 from alphaville/feature/333-affine-space-codegen

Support for affine spaces in code generation

Cargo.toml

@@ -98,7 +98,8 @@ rpmalloc = { version = "0.2", features = [
 # epigraph of the squared Euclidean norm
 roots = "0.0.8"
 
-# Least squares solver
+# Least squares solver (NOTE: ndarray must be version 0.15 - not 0.16)
+# Bug report: https://github.com/argmin-rs/modcholesky/issues/34
 ndarray = { version = "0.15", features = ["approx"] }
 modcholesky = "0.1"
 

docs/python-interface.md

@@ -77,18 +77,21 @@ following types of constraints:
 
 | Constraint         | Explanation                                    |
 |--------------------|------------------------------------------------|
-| `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]])`. |
-| `NoConstraints`    | No constraints - the whole $\mathbb{R}^{n}$|
-| `Rectangle`        | Rectangle, $$R = \\{u \in \mathbb{R}^{n_u} {}:{} f_{\min} \leq u \leq f_{\max}\\},$$ for example, `Rectangle(fmin, fmax)` |
-| `SecondOrderCone`  | Second-order aka "ice cream" aka "Lorenz" cone |
-| `EpigraphSquaredNorm`| The epigraph of the squared Eucliden norm is a set of the form $X = \\{(z, t) \in \mathbb{R}^{n+1}:  \Vert z \Vert  \leq t\\}$. | 
-| `CartesianProduct` | Cartesian product of any of the above. See more information below. |
+| `AffineSpace`      | An affine space is a set of the form $\\{x\in\mathbb{R}^n {}:{} Ax = b\\}$ for a matrix $A\in \mathbb{R}^p$ and vector $b$.  Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.AffineSpace.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.affine_space)) |
+| `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) Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Ball2.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.ball2)) |
+| `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) Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.BallInf.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.ball_inf)) |
+| `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) Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Ball1.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.ball1)) |
+| `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) Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Sphere2.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.sphere2)) |
+| `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). Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Simplex.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.simplex)) |
+| `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 straightforward: `Halfspace(c, b)`. Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Halfspace.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.halfspace)) |
+| `Hyperplane` | A hyperplane is a set given by $H=\\{x \in \mathbb{R}^n {}:{} c^\intercal x =b \\}$. Docs: [Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Hyperplane.html) |
+| `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]])`. Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.FiniteSet.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.finite_set)) |
+| `NoConstraints`    | No constraints - the whole $\mathbb{R}^{n}$ Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.NoConstraints.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.no_constraints)) |
+| `Rectangle`        | Rectangle, $$R = \\{u \in \mathbb{R}^{n_u} {}:{} f_{\min} \leq u \leq f_{\max}\\},$$ for example, `Rectangle(fmin, fmax)` Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Rectangle.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.rectangle)) |
+| `SecondOrderCone`  | Second-order aka "ice cream" aka "Lorenz" cone. Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.SecondOrderCone.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.soc)) |
+| `EpigraphSquaredNorm`| The epigraph of the squared Euclidean norm is a set of the form $X = \\{(z, t) \in \mathbb{R}^{n+1}:  \Vert z \Vert  \leq t\\}$. Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.EpigraphSquaredNorm.html), Python: to be implemented) |
+| `CartesianProduct` | Cartesian product of any of the above. See more information below. Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.CartesianProduct.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.cartesian)) |
+| `Zero` | The set $\\{0\\}$. Docs: ([Rust](https://docs.rs/optimization_engine/latest/optimization_engine/constraints/struct.Zero.html), [Python](https://alphaville.github.io/optimization-engine/api-dox/html/opengen.constraints.html#module-opengen.constraints.zero)) |
 
 
 
@@ -424,7 +427,7 @@ Note here that the solver performed 6 outer
 [penalty-type iterations](https://en.wikipedia.org/wiki/Penalty_method)
 and 35 inner iterations overall. The infinity norm of the constraint
 violations is approximately $7.66\cdot 10^{-5}$ (which is below the
-default tolerance of $10^{-4}$. This means that the solution $u^\star$
+default tolerance of $10^{-4}$.) This means that the solution $u^\star$
 satisfies
 
 <div class="math">

open-codegen/CHANGELOG.md

@@ -7,6 +7,15 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 Note: This is the Changelog file of `opengen` - the Python interface of OpEn
 
+## [0.9.0] - 2024-08-15
+
+### Added 
+
+* Support for affine spaces in code generation (see `opengen.constraints.AffineSpace`)
+* TCP solver status can be printed
+
+
+
 ## [0.8.1] - 2024-07-17
 
 ### Added
@@ -193,6 +202,8 @@ Note: This is the Changelog file of `opengen` - the Python interface of OpEn
 * Project-specific `tcp_iface` TCP interface
 * Fixed `lbfgs` typo
 
+[0.9.0]: https://github.com/alphaville/optimization-engine/compare/opengen-0.8.1...opengen-0.9.0
+[0.8.1]: https://github.com/alphaville/optimization-engine/compare/v0.9.0...opengen-0.8.1
 [0.8.1]: https://github.com/alphaville/optimization-engine/compare/v0.9.0...opengen-0.8.1
 [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

open-codegen/VERSION

@@ -1 +1 @@
-0.8.1
+0.9.0

open-codegen/opengen/builder/optimizer_builder.py

@@ -26,6 +26,11 @@
 _ICASADI_PREFIX = 'icasadi_'
 _ROS_PREFIX = 'ros_node_'
 
+# Template files
+_OPTIMIZER_RS = "optimizer.rs.jinja"
+_OPTIMIZER_CARGO_TOML = "optimizer_cargo.toml.jinja"
+_OPTIMIZER_BUILD_RS = "optimizer_build.rs.jinja"
+
 
 def make_dir_if_not_exists(directory):
     if not os.path.exists(directory):
@@ -243,7 +248,7 @@ def __generate_cargo_toml(self):
         self.__logger.info("Generating Cargo.toml for target optimizer")
         target_dir = self.__target_dir()
         cargo_template = OpEnOptimizerBuilder.__get_template(
-            'optimizer_cargo.toml')
+            _OPTIMIZER_CARGO_TOML)
         cargo_output_template = cargo_template.render(
             meta=self.__meta,
             build_config=self.__build_config,
@@ -540,7 +545,7 @@ def __generate_main_project_code(self):
             "Generating main code for target optimizer (lib.rs)")
         target_dir = self.__target_dir()
         optimizer_rs_template = OpEnOptimizerBuilder.__get_template(
-            'optimizer.rs')
+            _OPTIMIZER_RS)
         optimizer_rs_output_template = optimizer_rs_template.render(
             solver_config=self.__solver_config,
             meta=self.__meta,
@@ -557,7 +562,7 @@ def __generate_build_rs(self):
         self.__logger.info("Generating build.rs for target optimizer")
         target_dir = self.__target_dir()
         build_rs_template = OpEnOptimizerBuilder.__get_template(
-            'optimizer_build.rs')
+            _OPTIMIZER_BUILD_RS)
         build_rs_output_template = build_rs_template.render(
             meta=self.__meta,
             activate_clib_generation=self.__build_config.build_c_bindings)

open-codegen/opengen/constraints/__init__.py

@@ -11,3 +11,4 @@
 from .finite_set import *
 from .halfspace import *
 from .simplex import *
+from .affine_space import *

open-codegen/opengen/constraints/affine_space.py

@@ -0,0 +1,56 @@
+import casadi.casadi as cs
+import numpy as np
+from .constraint import Constraint
+import opengen.functions as fn
+
+
+class AffineSpace(Constraint):
+    """An affine constraint
+
+    A constraint of the form :math:`Ax = b`, where :math:`A` and :math:`b` are 
+    a matrix and a vector of appropriate dimensions
+    """
+
+    def __init__(self, A, b):
+        """Constructor for an affine space
+
+        :return: new instance of AffineSpace
+        """
+        self.__A = A.flatten('C')
+        self.__b = b
+
+    @property
+    def matrix_a(self):
+        """Matrix A
+        """
+        return self.__A
+
+    @property
+    def vector_b(self):
+        """Vector b        
+        """
+        return self.__b
+
+    def distance_squared(self, u):
+        """Squared distance to affine space
+
+        Not implemented yet
+        """
+        raise NotImplementedError()
+
+    def project(self, u):
+        """Projection on affine space
+
+        Not implemented yet
+        """
+        raise NotImplementedError()
+
+    def is_convex(self):
+        """Affine spaces are convex sets
+        """
+        return True
+
+    def is_compact(self):
+        """Affine spaces are not compact sets
+        """
+        return False

open-codegen/opengen/tcp/solver_status.py

@@ -115,3 +115,15 @@ def cost(self):
         :return: Value of cost function at the solution
         """
         return self.__dict__["__cost"]
+
+    def __repr__(self):
+        return "Solver Status Report:\n" + \
+            f"Exit status....... {self.exit_status}\n" + \
+            f"Num Outer Iters... {self.num_outer_iterations}\n" + \
+            f"Num Inner Iters... {self.num_inner_iterations}\n" + \
+            f"FPR............... {self.last_problem_norm_fpr}\n" + \
+            "Infeasibility\n" + \
+            f"  L  F1............ {self.f1_infeasibility}\n" + \
+            f"  L  Fw............ {self.f2_norm}\n" + \
+            f"Penalty............ {self.penalty}\n" + \
+            f"Time............... {self.solve_time_ms} ms\n"

open-codegen/opengen/templates/icasadi/icasadi_lib.rs

@@ -268,7 +268,7 @@ pub fn precondition(
 
 /// Computes the initial penalty
 ///
-/// Make sure that you have called init_
+/// Make sure that you have called init_{{ meta.optimizer_name }}
 pub fn initial_penalty(
     u: &[f64],
     static_params: &[f64],

open-codegen/opengen/templates/optimizer.rs → open-codegen/opengen/templates/optimizer.rs.jinja

@@ -61,17 +61,18 @@ pub const {{meta.optimizer_name|upper}}_N1: usize = {{problem.dim_constraints_au
 /// Number of penalty constraints
 pub const {{meta.optimizer_name|upper}}_N2: usize = {{problem.dim_constraints_penalty() or 0}};
 
-{% include "c/optimizer_cinterface.rs" %}
+{% include "c/optimizer_cinterface.rs.jinja" %}
 
 // ---Parameters of the constraints----------------------------------------------------------------------
-
+{# CASE I: Ball* or Sphere #}
 {% 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 %};
-
 /// Constraints: Radius of Ball
 const CONSTRAINTS_BALL_RADIUS : f64 = {{problem.constraints.radius}};
-{% elif 'Rectangle' == problem.constraints.__class__.__name__ -%}
+{% endif %}
+{# CASE II: Rectangle #}
+{% if 'Rectangle' == problem.constraints.__class__.__name__ -%}
 const CONSTRAINTS_XMIN :Option<&[f64]> = {% if problem.constraints.xmin is not none %}Some(&[
 {%- for xmini in problem.constraints.xmin -%}
 {%- if float('-inf') == xmini -%}std::f64::NEG_INFINITY{%- else -%}{{xmini}}{%- endif -%},
@@ -85,7 +86,6 @@ const CONSTRAINTS_XMAX :Option<&[f64]> = {% if problem.constraints.xmax is not n
 {% endif %}
 
 
-
 {% if problem.alm_set_c is not none %}
 // ---Parameters of ALM-type constraints (Set C)---------------------------------------------------------
 {% if 'Ball2' == problem.alm_set_c.__class__.__name__ or 'BallInf' == problem.alm_set_c.__class__.__name__ -%}
@@ -150,6 +150,10 @@ fn make_constraints() -> impl Constraint {
     {% elif 'Rectangle' == problem.constraints.__class__.__name__ -%}
     // - Rectangle:
     Rectangle::new(CONSTRAINTS_XMIN, CONSTRAINTS_XMAX)
+    {% elif 'AffineSpace' == problem.constraints.__class__.__name__ -%}
+    let constraints_affine_a = vec![{{problem.constraints.matrix_a | join(', ')}}];
+    let constraints_affine_b  = vec![{{problem.constraints.vector_b | join(', ')}}];
+    AffineSpace::new(constraints_affine_a, constraints_affine_b)
     {% elif 'FiniteSet' == problem.constraints.__class__.__name__ -%}
     // - Finite Set:
     let data: &[&[f64]] = &[
@@ -192,6 +196,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}} = Sphere2::new(center_{{loop.index}}, radius_{{loop.index}});
         let bounds = bounds.add_constraint(idx_{{loop.index}}, set_{{loop.index}});
+        {% elif 'AffineSpace' == set_i.__class__.__name__ -%}
+        let constraints_affine_a = vec![{{problem.constraints.matrix_a | join(', ')}}];
+        let constraints_affine_b  = vec![{{problem.constraints.vector_b | join(', ')}}];
+        let set_{{loop.index}} = AffineSpace::new(constraints_affine_a, constraints_affine_b)
+        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}});