diff --git a/doc/manual/development/link_v1.rst b/doc/manual/development/link_v1.rst new file mode 100644 index 000000000..7f2eaf30d --- /dev/null +++ b/doc/manual/development/link_v1.rst @@ -0,0 +1,8 @@ +.. _sec-dev-v1: + +Old version (Codac V1) +====================== + +.. raw:: html + + diff --git a/doc/manual/index.rst b/doc/manual/index.rst index dc2e192c7..e4f545767 100644 --- a/doc/manual/index.rst +++ b/doc/manual/index.rst @@ -180,6 +180,8 @@ Overview of Codac User manual ----------- +* :ref:`sec-intro` + * :ref:`sec-install` * :ref:`sec-install-py` * :ref:`sec-install-cpp` @@ -334,6 +336,7 @@ User manual * :ref:`sec-tools` * :ref:`sec-tools-serialization` * :ref:`sec-tools-registration` + * :ref:`sec-tools-octasym` * Codac extensions * :ref:`sec-extensions-capd` @@ -374,6 +377,7 @@ Development * :ref:`sec-dev-common-issues` * :ref:`sec-dev-changelog` * :ref:`sec-dev-api` +* :ref:`sec-dev-v1` @@ -392,13 +396,13 @@ Development :caption: User manual :maxdepth: 2 + manual/introduction/index.rst manual/installation/index.rst manual/intervals/index.rst manual/linear/index.rst manual/functions/index.rst manual/contractors/index.rst manual/geometry/index.rst - manual/actions/index.rst manual/ellipsoids/index.rst manual/visualization/index.rst manual/tools/index.rst @@ -453,6 +457,7 @@ Development development/common_issues.rst development/changelog.rst development/api_redirect.rst + development/link_v1.rst How to cite Codac diff --git a/doc/manual/manual/actions/index.rst b/doc/manual/manual/actions/index.rst deleted file mode 100644 index eb7bde9d3..000000000 --- a/doc/manual/manual/actions/index.rst +++ /dev/null @@ -1,11 +0,0 @@ -.. _sec-actions: - -Actions -======= - -Actions are in codac - -.. toctree:: - :maxdepth: 1 - - octasym.rst \ No newline at end of file diff --git a/doc/manual/manual/contractors/analytic/ctcinverse.rst b/doc/manual/manual/contractors/analytic/ctcinverse.rst index b893b19cb..600c3bf2a 100644 --- a/doc/manual/manual/contractors/analytic/ctcinverse.rst +++ b/doc/manual/manual/contractors/analytic/ctcinverse.rst @@ -5,7 +5,188 @@ The CtcInverse contractor Main author: `Simon Rohou `_ -Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p`. The ``CtcInverse`` contractor allows to deal with constraints under the form :math:`\mathbf{f}(\mathbf{x})\in[\mathbf{y}]` by contracting boxes :math:`[\mathbf{x}]\in\mathbb{IR}^n`. +Definition +---------- + +Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p`. The ``CtcInverse`` contractor allows one to handle constraints of the form :math:`\mathbf{f}(\mathbf{x})\in\mathbb{Y}` by contracting input boxes :math:`[\mathbf{x}]\in\mathbb{IR}^n` while preserving all solutions consistent with the constraint set :math:`\mathbb{Y}`. + +:math:`\mathbb{Y}` is a constraint set on the codomain of :math:`\mathbf{f}`. Note that contractors are typically used when :math:`\mathbb{Y}` is thin (*e.g.*, equality constraints). For thick sets (inequalities), separators are often preferable to obtain both outer and inner approximations, and faster computations. + +In Codac, :math:`\mathbb{Y}` can be represented in two main ways: + +.. important:: + + **(i) Box representation** + + When :math:`\mathbb{Y}` is given as a box :math:`[\mathbf{y}]`, the constraint reads: :math:`\mathbf{f}(\mathbf{x}) \in [\mathbf{y}]`. + + **(ii) Contractor representation** + + More generally, :math:`\mathbb{Y}` may be any output set representable by a contractor + :math:`\mathcal{C}_{\mathbb{Y}}`: :math:`\mathbf{f}(\mathbf{x}) \in \mathcal{C}_{\mathbb{Y}}`. + + In this case, :math:`\mathcal{C}_{\mathbb{Y}}` is any contractor acting on boxes in :math:`\mathbb{R}^p`. + +In the current version of Codac (v2), ``CtcInverse`` is the successor of the v1 ``CtcFunction`` contractor (renamed for clarity). It is also closely related to the ``CtcFwdBwd`` of the `IBEX library `_. + + +Construction and basic usage +---------------------------- + +The typical workflow is: + +1. Define analytic variables (scalar, vector, matrix) associated with the domain of the function. +2. Build an :class:`~codac2.AnalyticFunction`. +3. Instantiate ``CtcInverse`` with a codomain constraint (a box :math:`[\mathbf{y}]` or a contractor :math:`\mathcal{C}_{\mathbb{Y}}`). +4. Contract an input box :math:`[\mathbf{x}]` or pave the contractor in the domain space. + + +Dealing with equalities and inequalities +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Equalities are expressed with a degenerate interval (or box) in the codomain, *e.g.*: + +- :math:`f(x)=0` is encoded as :math:`f(x)\in[0,0]`. + +Inequalities are encoded by widening the codomain interval(s). For instance: + +- :math:`f(x)\leqslant 0` is encoded as :math:`f(x)\in[-\infty,0]`. + +In the case of inequalities, the use of the equivalent ``SepInverse`` separator should be preferred for efficiency reasons or in order to obtain an inner approximation. + + +Example +^^^^^^^ + +Consider for instance Himmelblau's function: + +.. math:: + + f(\mathbf{x}) = (x_1^2 + x_2 - 11)^2 + (x_1 + x_2^2 - 7)^2. + +To represent the vectors :math:`\mathbf{x}\in\mathbb{R}^n` consistent with the constraint :math:`f(\mathbf{x})=50`, a ``CtcInverse`` contractor can be built as follows: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcinv-1-beg] + :end-before: [ctcinv-1-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcinv-1-beg] + :end-before: [ctcinv-1-end] + :dedent: 4 + +The contractor can be used as an operator to contract a 2d box :math:`[\mathbf{x}]`. It can also be involved in a paver in order to reveal the constraint: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcinv-2-beg] + :end-before: [ctcinv-2-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcinv-2-beg] + :end-before: [ctcinv-2-end] + :dedent: 4 + +Which produces the following output: + +.. figure:: ./himmelblau_50.png + :width: 400px + + Outer approximation of the solution set for :math:`f(\mathbf{x})\in[50,50]`. This paving result reveals three connected subsets approximated with outer boxes. Blue parts are guaranteed not to contain solutions. + +We recall that for thick solution sets, one should prefer the use of the ``SepInverse`` separator in order to get both outer and inner approximations. For instance, the solution set associated with :math:`f(\mathbf{x})\leqslant 50` corresponds to: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcinv-3-beg] + :end-before: [ctcinv-3-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcinv-3-beg] + :end-before: [ctcinv-3-end] + :dedent: 4 + +.. figure:: ./himmelblau_50_inner.png + :width: 400px + + Inner and outer approximations of the solution set for :math:`f(\mathbf{x})\leqslant 50`, equivalent to :math:`f(\mathbf{x})\in[0,50]`, obtained by paving a separator ``SepInverse(f, [0,50])``. Green parts are boxes containing only vectors consistent with the constraint. + +Finally, it is easy to consider disjunction of constraints by making a union of contractors. For instance, the solution set + +.. math:: + :label: himmelblau_union + + \mathbb{X}=\left\{\mathbf{x}, \left(f(\mathbf{x}) = 50\right) \lor \left(f(\mathbf{x}) = 150\right) \lor \left(f(\mathbf{x}) = 250\right)\right\} + +can be easily approximated by the following union of contractors: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcinv-4-beg] + :end-before: [ctcinv-4-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcinv-4-beg] + :end-before: [ctcinv-4-end] + :dedent: 4 + +.. figure:: ./himmelblau_50_150_250.png + :width: 400px + + Outer approximation of the solution set :math:`\mathbb{X}` for :eq:`himmelblau_union`. + + +.. from codac import * +.. +.. # Example of Himmelblau's function +.. a = 11; b = 7 +.. x = VectorVar(2) +.. f = AnalyticFunction([x], sqr(sqr(x[0])+x[1]-a)+sqr(x[0]+sqr(x[1])-b)) +.. +.. c1 = CtcInverse(f, 50) +.. c2 = CtcInverse(f, 150) +.. c3 = CtcInverse(f, 250) +.. +.. g = Figure2D("test", GraphicOutput.VIBES | GraphicOutput.IPE) +.. g.set_axes(axis(0,[-6,6]), axis(1,[-6,6])) +.. g.pave([[-6,6],[-6,6]], c1|c2|c3, 1e-2) + + +Possible usages of ``CtcInverse`` with respect to domain and codomain types +----------------------------------------------------------------------------- :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p` is the most classical case, but Codac provides generic solutions to define a wide variety of expressions such as functions defined as :math:`\mathbf{f}:X\to Y`, where: @@ -14,7 +195,7 @@ Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p`. The ``CtcIn * :math:`X=\mathbb{R}`: one scalar input ; * :math:`X=\mathbb{R}^n`: one :math:`n`-d vector input ; * :math:`X=\mathbb{R}^{n\times m}`: one :math:`n\times m` matrix input ; - * :math:`X=\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}`: several types of mixed input, with no restrictions on order or number. + * :math:`X=\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}`: several types of mixed inputs, with no restrictions on order or number. - :math:`Y` denotes the codomain of :math:`\mathbf{f}`, that can be either :math:`\mathbb{R}`, :math:`\mathbb{R}^p` or :math:`\mathbb{R}^{p\times q}`. @@ -36,7 +217,7 @@ Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p`. The ``CtcIn The following table lists the combinations implemented in Codac. The :bg-ok:`✓` cases correspond to efficient contractions involving centered form expressions, while :bg-bok:`✓` refer to classical forward/backward propagations. -Note that since Codac's design is based on template programming, the most generic use-cases of the ``CtcInverse`` are not available in Python or Matlab. +Note that since Codac's design is based on C++ template programming, the most generic use-cases of the ``CtcInverse`` are not available in Python or Matlab. .. tabs:: @@ -95,4 +276,75 @@ Note that since Codac's design is based on template programming, the most generi | | :math:`\mathbb{R}^{r\times c}` | |noa| | |noa| | |noa| | | +--------------------------------+-------------------------+---------------------------+-------------------------------------+ | | any mixed types | |noa| | |noa| | |noa| | - +--------------------+--------------------------------+-------------------------+---------------------------+-------------------------------------+ \ No newline at end of file + +--------------------+--------------------------------+-------------------------+---------------------------+-------------------------------------+ + + +Not-in constraints: ``CtcInverseNotIn`` +--------------------------------------- + +When the constraint is a complement constraint :math:`\mathbf{f}(\mathbf{x})\notin\mathbb{Y}`, Codac provides ``CtcInverseNotIn`` which internally builds a union of inverse contractors over the complementary parts of :math:`\mathbb{Y}`. + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcinv-5-beg] + :end-before: [ctcinv-5-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcinv-5-beg] + :end-before: [ctcinv-5-end] + :dedent: 4 + + +Miscellaneous +------------- + +Access to the underlying function +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The underlying analytic function can be accessed through ``.function()`` (useful for dimension checks or meta-programming). + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcinv-6-beg] + :end-before: [ctcinv-6-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcinv-6-beg] + :end-before: [ctcinv-6-end] + :dedent: 4 + +Centered form option +^^^^^^^^^^^^^^^^^^^^ + +For some signatures (see :bg-ok:`✓` in the table), ``CtcInverse`` can leverage centered-form +expressions to improve contraction. +The constructor exposes: + +- ``with_centered_form`` (default: ``true``): enable/disable centered-form based steps. + +Disabling it can be useful to reduce computation cost when forward/backward propagation is +sufficient for your application, or for comparisons with the centered form. + + +.. admonition:: Technical documentation + + See the C++ API documentation: + + - ``CtcInverse``: `codac2_CtcInverse.h <../../../api/html/codac2___ctc_inverse_8h.html>`_ + - ``CtcInverseNotIn``: `codac2_CtcInverseNotIn.h <../../../api/html/codac2___ctc_inverse_not_in_8h.html>`_ \ No newline at end of file diff --git a/doc/manual/manual/contractors/analytic/himmelblau_50.png b/doc/manual/manual/contractors/analytic/himmelblau_50.png new file mode 100644 index 000000000..239f44eaa Binary files /dev/null and b/doc/manual/manual/contractors/analytic/himmelblau_50.png differ diff --git a/doc/manual/manual/contractors/analytic/himmelblau_50_150_250.png b/doc/manual/manual/contractors/analytic/himmelblau_50_150_250.png new file mode 100644 index 000000000..aacf4851b Binary files /dev/null and b/doc/manual/manual/contractors/analytic/himmelblau_50_150_250.png differ diff --git a/doc/manual/manual/contractors/analytic/himmelblau_50_inner.png b/doc/manual/manual/contractors/analytic/himmelblau_50_inner.png new file mode 100644 index 000000000..471bdf1f0 Binary files /dev/null and b/doc/manual/manual/contractors/analytic/himmelblau_50_inner.png differ diff --git a/doc/manual/manual/contractors/analytic/src.cpp b/doc/manual/manual/contractors/analytic/src.cpp new file mode 100644 index 000000000..1d9d4cf02 --- /dev/null +++ b/doc/manual/manual/contractors/analytic/src.cpp @@ -0,0 +1,74 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("CtcInverse - manual") +{ + { + // [ctcinv-1-beg] + // Example of Himmelblau's function + double a = 11, b = 7; + VectorVar x(2); + AnalyticFunction f({x}, sqr(sqr(x[0])+x[1]-a)+sqr(x[0]+sqr(x[1])-b)); + CtcInverse c(f, 50); + // [ctcinv-1-end] + + // [ctcinv-2-beg] + DefaultFigure::pave({{-6,6},{-6,6}}, c, 1e-2); + // [ctcinv-2-end] + + // [ctcinv-3-beg] + SepInverse s(f, {0,50}); + DefaultFigure::pave({{-6,6},{-6,6}}, s, 1e-2); + // [ctcinv-3-end] + + // [ctcinv-4-beg] + auto cu = CtcInverse(f,50) | CtcInverse(f,150) | CtcInverse(f,250); + DefaultFigure::pave({{-6,6},{-6,6}}, cu, 1e-2); + // [ctcinv-4-end] + } + + { + // [ctcinv-5-beg] + VectorVar x(2); + AnalyticFunction f({x}, x[0]); + + // Enforce the first component not in [0,1] + CtcInverseNotIn c(f, {0,1}); + + IntervalVector y({{0.5,3},{-1,1}}); + c.contract(y); // {{1,3},{-1,1}} + // Only the first component is constrained by the not-in condition + // [ctcinv-5-end] + + CHECK(y == IntervalVector({{1,3},{-1,1}})); + } + + { + // [ctcinv-6-beg] + VectorVar x(2); + AnalyticFunction f({x}, x[0]-x[1]); + CtcInverse c(f, 0); + // c.function().input_size() == 2 + // c.function().output_size() == 1 + // [ctcinv-6-end] + + CHECK(c.function().input_size() == 2); + CHECK(c.function().output_size() == 1); + } +} \ No newline at end of file diff --git a/doc/manual/manual/contractors/analytic/src.py b/doc/manual/manual/contractors/analytic/src.py new file mode 100644 index 000000000..34fe2214a --- /dev/null +++ b/doc/manual/manual/contractors/analytic/src.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2026 +# \author Simon Rohou +# \copyright Copyright 2026 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import sys, os +import unittest +import math +from codac import * + +class TestCtcAnalyticManual(unittest.TestCase): + + def tests_CtcInverse_manual(test): + + # [ctcinv-1-beg] + # Example of Himmelblau's function + a = 11; b = 7 + x = VectorVar(2) + f = AnalyticFunction([x], sqr(sqr(x[0])+x[1]-a)+sqr(x[0]+sqr(x[1])-b)) + c = CtcInverse(f, 50) + # [ctcinv-1-end] + + # [ctcinv-2-beg] + DefaultFigure.pave([[-6,6],[-6,6]], c, 1e-2) + # [ctcinv-2-end] + + # [ctcinv-3-beg] + s = SepInverse(f, [0,50]) + DefaultFigure.pave([[-6,6],[-6,6]], s, 1e-2) + # [ctcinv-3-end] + + # [ctcinv-4-beg] + cu = CtcInverse(f,50) | CtcInverse(f,150) | CtcInverse(f,250) + DefaultFigure.pave([[-6,6],[-6,6]], cu, 1e-2) + # [ctcinv-4-end] + + # [ctcinv-5-beg] + x = VectorVar(2) + f = AnalyticFunction([x], x[0]) + + # Enforce the first component not in [0,1] + c = CtcInverseNotIn(f, [0,1]) + + y = IntervalVector([[0.5,3],[-1,1]]) + c.contract(y) # [[1,3],[-1,1]] + # Only the first component is constrained by the not-in condition + # [ctcinv-5-end] + + test.assertTrue(y == IntervalVector([[1,3],[-1,1]])) + + # [ctcinv-6-beg] + x = VectorVar(2) + f = AnalyticFunction([x], x[0]-x[1]) + c = CtcInverse(f, 0) + assert c.function().input_size() == 2 + assert c.function().output_size() == 1 + # [ctcinv-6-end] + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/doc/manual/manual/functions/analytic/analytic_functions.rst b/doc/manual/manual/functions/analytic/analytic_functions.rst index 9dd4db277..1752148f7 100644 --- a/doc/manual/manual/functions/analytic/analytic_functions.rst +++ b/doc/manual/manual/functions/analytic/analytic_functions.rst @@ -19,7 +19,7 @@ This page provides an overview of the ``AnalyticFunction`` class, its key featur .. note:: - .. Figure:: CtcInverse_small.png + .. figure:: CtcInverse_small.png :align: right For defining a contractor based on an ``AnalyticFunction``, the ``CtcInverse`` class is available. diff --git a/doc/manual/manual/intervals/IntervalVector_class.rst b/doc/manual/manual/intervals/IntervalVector_class.rst index fff8fe17d..2f441c74a 100644 --- a/doc/manual/manual/intervals/IntervalVector_class.rst +++ b/doc/manual/manual/intervals/IntervalVector_class.rst @@ -3,26 +3,248 @@ The IntervalVector class ======================== -Some specific commands in Python are provided below: + Main author: `Simon Rohou `_ + +An :class:`~codac.IntervalVector` represents an axis-aligned *box* (Cartesian product of intervals) in :math:`\mathbb{R}^n`: + +.. math:: + + [\mathbf{x}] = [x_1]\times [x_2]\times \dots \times [x_n]. + +In C++, ``IntervalVector`` is an alias for an Eigen dynamic column vector whose scalar type is :ref:`Interval ` (i.e., ``Eigen::Matrix``). This means you can use the usual Eigen vector API (size/resize/segment/concatenation patterns) while benefiting from Codac’s interval arithmetic and utilities. In Python and Matlab, some of these functions are available. + +Creating IntervalVectors +------------------------ + +A box can be created from: + +* a dimension :math:`n` (components initialized with the default interval :math:`[-\infty,\infty]`), +* a dimension :math:`n` and a common interval (a “cube”), +* a list of intervals or bounds pairs, +* a real vector (degenerate box wrapping a point). + +.. tabs:: + + .. group-tab:: Python + + .. code-block:: py + + # Default box: [-oo,oo]^n + x = IntervalVector(3) + + # Cube [-1,3]^2 + y = IntervalVector(2, Interval(-1,3)) + + # From a list of bounds (each entry is [lb,ub]) + z = IntervalVector([[3,4],[4,6]]) # [3,4]×[4,6] + + # From a list of components (Intervals and/or bounds pairs) + q = IntervalVector([y[1], z[0], [0,oo]]) # [-1,3]×[3,4]×[0,oo] + + # From a point (degenerate intervals) + p = Vector([0.42,0.42,0.42]) + bp = IntervalVector(p) # [0.42,0.42]^3 + + .. group-tab:: C++ + + .. code-block:: c++ + + // Default box: [-oo,oo]^n (Interval default constructor) + IntervalVector x(3); + + // Cube [-1,3]^2 + IntervalVector y(2, Interval(-1,3)); + + // From a list of bounds (initializer-list style) + IntervalVector z{{3,4},{4,6}}; // [3,4]×[4,6] + + // From a point (degenerate intervals) + Vector p({0.42,0.42,0.42}); + IntervalVector bp(p); // [0.42,0.42]^3 + +.. note:: + + Codac intentionally does **not** enable Eigen’s “initialize matrices by zero” option for interval-based vectors and matrices, because default-initializing an interval should produce :math:`[-\infty,\infty]` rather than ``0``. + As a result, resizing an ``IntervalVector`` typically creates new components initialized as :math:`[-\infty,\infty]`. + +Accessing and modifying components +---------------------------------- + +Components are intervals; indexing is 0-based in Python/C++ and 1-based in Matlab. .. tabs:: - - .. code-tab:: py - x = IntervalVector([[1,2],[2,3],[3,4]]) - y = IntervalVector(3) + .. group-tab:: Python + + .. code-block:: py + + x = IntervalVector(2, [-1,3]) # [-1,3]^2 + x[1] = Interval(0,10) # [-1,3]×[0,10] + + # Iterating/accessing over components + y = IntervalVector(2) + for i, xi in enumerate(x): + y[i] = xi + + # Unpacking (Python convenience) + a,b = x + assert a == x[0] and b == x[1] + + # Building a new box from existing components + v = IntervalVector([*x, [3,6]]) # concatenation in Python + # v == [[-1,3]×[0,10]×[3,6]] + + # Resize: new components are default-initialized ([-oo,oo]) + v.resize(4) # v == [[-1,3]×[0,10]×[3,6]×[-oo,oo]] + s = v.subvector(1,2) # [0,10]×[3,6] + + .. group-tab:: C++ + + .. code-block:: c++ + + IntervalVector x(2, Interval(-1,3)); // [-1,3]^2 + x[1] = Interval(0, 10); // [-1,3]×[0,10] + + // Accessing components + const Interval& x0 = x[0]; + + // Resize: new components are default-initialized ([-oo,oo]) + x.resize(4); // x == [-1,3]×[0,10]×[-oo,oo]×[-oo,oo] + + // Subvector / segment extraction + IntervalVector s = x.subvector(1,2); // [0,10]×[-oo,oo] + +Box properties +-------------- + +Most scalar properties of :ref:`Interval ` have natural extensions to boxes: + +* component-wise bounds: :math:`\mathbf{x}^-` and :math:`\mathbf{x}^+`, +* component-wise widths (diameters/radii), +* geometric metrics such as volume or max diameter. + +Typical accessors you will use in practice: + +.. tabs:: + + .. group-tab:: Python + + .. code-block:: py + + x = IntervalVector([[0,2],[-1,3]]) + + n = x.size() # dimension + # Common box information (component-wise): + lo = x.lb() # Vector of lower bounds + hi = x.ub() # Vector of upper bounds + m = x.mid() # Vector of midpoints + d = x.diam() # Vector of diameters + + .. group-tab:: C++ + + .. code-block:: c++ + + IntervalVector x{{0,2},{-1,3}}; + + Index n = x.size(); + Vector lo = x.lb(); + Vector hi = x.ub(); + Vector m = x.mid(); + Vector d = x.diam(); + +.. note:: + + The exact set of box-wide metrics available (*e.g.*, ``max_diam()``, ``volume()``, ``argmax_diam()``) is provided by the interval-based Eigen extensions used in Codac. See the technical documentation section at the end of this page for the comprehensive list. + +Testing boxes +------------- + +Common predicates include: + +- ``is_empty()``: tests whether the box represents the empty set (at least one empty component), +- ``is_degenerated()``: all components are degenerate, +- ``is_degenerated()``: at least one component is degenerate, +- ``is_unbounded()``: at least one bound is infinite, +- ``contains(p)``: inclusion test of a point/vector, +- ``intersects(y)``: non-empty intersection test with another box, +- ``is_subset(y)``, ``is_strict_subset(y)``, *etc.* + +.. tabs:: + + .. group-tab:: Python + + .. code-block:: py + + x = IntervalVector([[0,1],[2,3]]) + y = IntervalVector([[0.5,2],[1,4]]) + + assert x.intersects(y) + assert x.is_subset(y) + + .. group-tab:: C++ + + .. code-block:: c++ + + IntervalVector x{{0,1},{2,3}}; + IntervalVector y{{0.5,2},{1,4}}; + + assert(x.intersects(y)); + assert(x.is_subset(y)); + +Advanced operations +------------------- + +As for intervals, typical advanced operations on boxes include: + +.. list-table:: Common advanced methods for a given box :math:`[\mathbf{x}]` + :widths: 30 70 + :header-rows: 1 + + * - Method + - Description + * - ``inflate(rad)`` + - Expands the box by ``±rad`` (scalar radius for common inflation, or vector radius for specifying different inflation on each component). + * - ``bisect([ratio])`` + - Splits the box into two sub-boxes along the widest dimension. + * - ``rand()`` + - Returns a random sample ``Vector`` inside the box (when non-empty). + * - ``init()`` + - Re-initializes to :math:`[-\infty,\infty]^n`. + +Arithmetic and linear algebra +----------------------------- + +Because ``IntervalVector`` is a vector of ``Interval``, standard arithmetic is naturally extended component-wise: + +* box–box and box–scalar operations (``+``, ``-``, ``*``, ``/``), +* interactions with real vectors/matrices, +* matrix–vector products involving :class:`~codac.IntervalMatrix`. + +.. tabs:: + + .. group-tab:: Python + + .. code-block:: py + + x = IntervalVector([[0,1],[2,3]]) + y = IntervalVector([[1,2],[0,1]]) + + z1 = x+y + z2 = 2*x + z3 = x-1 + + .. group-tab:: C++ + + .. code-block:: c++ + + IntervalVector x{{0,1},{2,3}}; + IntervalVector y{{1,2},{0,1}}; - i = 0 - for xi in x: - y[i] = xi - i = i+1 + IntervalVector z1 = x+y; + IntervalVector z2 = 2.*x; + IntervalVector z3 = x-1.; - # x == y - a,b,c = x - # a == x[0] - # b == x[1] - # c == x[2] +.. admonition:: Technical documentation - v = IntervalVector([*x, [3,6]]) - # v == [[1,2],[2,3],[3,4],[3,6]] \ No newline at end of file + See the `C++ API documentation of the IntervalVector class <../../api/html/codac2___interval_vector_8h.html>`_. \ No newline at end of file diff --git a/doc/manual/manual/intervals/Interval_class.rst b/doc/manual/manual/intervals/Interval_class.rst index 25aafb0c3..48b24989a 100644 --- a/doc/manual/manual/intervals/Interval_class.rst +++ b/doc/manual/manual/intervals/Interval_class.rst @@ -193,7 +193,27 @@ Advanced operations Interval arithmetic ------------------- -All standard arithmetic operations are supported, both element-wise and with real numbers. +Interval analysis is based on the extension of all classical real arithmetic operators. +Consider two intervals :math:`[x]` and :math:`[y]` and an operator :math:`\diamond\in\left\{+,-,\cdot,/\right\}`. We define :math:`[x]\diamond[y]` as the smallest interval containing all feasible values for :math:`x\diamond y`, assuming that :math:`x\in[x]` and :math:`y\in[y]`: + +.. math:: + + [x]\diamond[y] = \big[\left\{x\diamond y \mid x\in[x],y\in[y]\right\}\big]. + +Dealing with closed intervals, most of the operations can rely on their bounds. It is for instance the case of addition, difference, union, *etc.*: + +.. math:: + + \begin{eqnarray} + [x]+[y]&=&\left[x^-+y^-,x^++y^+\right],\\ + \left[x\right]-\left[y\right]& = &\left[x^--y^+,x^+-y^-\right],\\ + \left[x\right]\sqcup\left[y\right]& = &\left[\min\left(x^-,y^-\right),\max\left(x^+,y^+\right)\right],\\ + \left[x\right]\cap\left[y\right]& = &\left[\max\left(x^-,y^-\right),\min\left(x^+,y^+\right)\right] \\ & & \textrm{if} \max\left\{x^-,y^-\right\}\leqslant\min\left\{x^+,y^+\right\}, \varnothing \textrm{ otherwise}. + \end{eqnarray} + +Note that Codac is built upon a low-level interval library, `GAOL `_, which has been built to provide functionalities for computing arithmetic on intervals, involving basic operations as well as non-linear functions. + +All standard arithmetic operations are supported in Codac, both element-wise and with real numbers. .. tabs:: @@ -304,4 +324,4 @@ These functions are useful in the context of interval arithmetic to tightly cont .. admonition:: Technical documentation - See the `C++ API documentation of this class <../../api/html/classcodac2_1_1_interval.html>`_. \ No newline at end of file + See the `C++ API documentation of the Interval class <../../api/html/classcodac2_1_1_interval.html>`_. \ No newline at end of file diff --git a/doc/manual/manual/intervals/index.rst b/doc/manual/manual/intervals/index.rst index 403845ee5..cdbc3828f 100644 --- a/doc/manual/manual/intervals/index.rst +++ b/doc/manual/manual/intervals/index.rst @@ -6,7 +6,7 @@ Intervals Codac provides data structures for handling basic interval sets. These structures represent the building blocks of most domain structures provided in Codac. The elementary interval structures are: - :ref:`The Interval class `: represents a real bounded interval :math:`[x^{-},x^{+}]`. -- ``IntervalVector`` (or ``IntervalRow``): represents a vector (or row) where each component is an interval. +- :ref:`The IntervalVector class ` (or ``IntervalRow``): represents a vector (or row) where each component is an interval. - ``IntervalMatrix``: represents a matrix where each element is an interval. diff --git a/doc/manual/manual/introduction/cst.png b/doc/manual/manual/introduction/cst.png new file mode 100644 index 000000000..e426c34ae Binary files /dev/null and b/doc/manual/manual/introduction/cst.png differ diff --git a/doc/manual/manual/introduction/ctc.png b/doc/manual/manual/introduction/ctc.png new file mode 100644 index 000000000..b92761b6b Binary files /dev/null and b/doc/manual/manual/introduction/ctc.png differ diff --git a/doc/manual/manual/introduction/index.rst b/doc/manual/manual/introduction/index.rst new file mode 100644 index 000000000..fc5cef30f --- /dev/null +++ b/doc/manual/manual/introduction/index.rst @@ -0,0 +1,93 @@ +.. _sec-intro: + +Introduction +============ + +Codac is designed to deal with static and dynamical systems that are usually encountered in robotics. +It stands on constraint programming: a paradigm in which the properties of a solution to be found (*e.g.* the pose of a robot, the location of a landmark) are stated by constraints coming from the equations of the problem. + +Then, a solver performs **constraint propagation** on the variables and provides a reliable set of feasible solutions corresponding to the considered problem. In this approach, the user concentrates on *what* is the problem instead of *how* to solve it, thus leaving the computer dealing with the *how*. The strength of this declarative paradigm lies in its simpleness, as it allows one to describe a complex problem without requiring the knowledge of resolution tools coming with specific parameters to choose. + +The following figure presents a conceptual view of a constraint propagation. + +.. figure:: ./cst.png + + In this theoretical view, a domain :math:`\mathbb{X}` depicted in yellow is known to enclose a solution :math:`\mathbf{x}^*` pictured by a red dot and consistent with two constraints :math:`\mathcal{L}_1` and :math:`\mathcal{L}_2`. The estimation of :math:`\mathbf{x}^*` consists in reducing :math:`\mathbb{X}` while satisfying :math:`\mathcal{L}_1` and :math:`\mathcal{L}_2`. After the propagation process, the set of feasible solutions has been reduced by removing vectors that were not consistent with the constraints: we obtain a finer approximation of :math:`\mathbf{x}^*`. + +.. When working with finite domains, a propagation technique can be used to simplify a problem. The process is run several times up to a fixed point reached when the domains cannot be reduced anymore. Interval analysis can be efficiently used for this purpose, taking advantage of interval arithmetic and its capacity to preserve any feasible solution. + +Variables +--------- + +The unknown solutions of a system are reals :math:`x`, vectors :math:`\mathbf{x}`, or in the dynamic case: trajectories denoted by :math:`x(\cdot)`, or :math:`\mathbf{x}(\cdot)` in the vector case. + +.. admonition:: Dot notation :math:`(\cdot)` + + Note that in the literature, the dot notation :math:`(\cdot)` may be used to represent the independent system variable (time). + This notation may be chosen in order to clearly distinguish a whole trajectory :math:`\mathbf{x}(\cdot):\mathbb{R}\to\mathbb{R}^n` from a local evaluation: :math:`\mathbf{x}(t)\in\mathbb{R}^n`. Indeed, the time :math:`t` may also be a variable to be estimated. + +In the constraint programming approach, the estimation of a variable consists in computing its reliable enclosure set, defined here as an interval of values. +The obtained set is said **reliable**: the resolution guarantees that no solution can be lost during the solving process, according to the equations defining the problem. In other words, a variable :math:`x` to be estimated will surely be enclosed in the resulting set :math:`[x]`. + + +Domains +------- + +As depicted in the figure, a variable is known to be enclosed in some domain :math:`\mathbb{X}` on which we will apply constraints. +In Codac, the domains are represented by the following sets: + +- for reals :math:`x`, domains are intervals :math:`[x]`; +- for vectors :math:`\mathbf{x}`, we define boxes :math:`[\mathbf{x}]`; +- for trajectories :math:`x(\cdot)`, we introduce tubes :math:`[x](\cdot)`. + + +Constraints +----------- + +The approach consists in defining the problem as a set of constraints. They usually come from state equations, numerical models, or measurements. +In mobile robotics, we usually have to deal with: + +- non-linear functions: :math:`x_3=\cos\big(x_1+\exp(x_2)\big)` +- differential equations: :math:`\dot{x}(\cdot)=f\big(x(\cdot),u(\cdot)\big)` +- temporal delays: :math:`x(t)=y(t-\tau)` +- time uncertainties: :math:`x(t)=y`, with :math:`t\in[t]` +- *etc.* + +The aim of Codac is to easily deal with these constraints in order to eventually characterize sets of variables compliant with the defined rules. +Each constraint comes from equations or measurements. And each constraint will be applied by means of an operator called **contractor**. + + +Contractors +----------- + +Contractors are operators designed to reduce domains without losing any solution consistent with the related constraint. In Codac, the contractors act on intervals, boxes and tubes. + +.. figure:: ./ctc.png + + Implementation of the constraints :math:`\mathcal{L}_1` and :math:`\mathcal{L}_2` of the previous figure by means of respective contractors :math:`\mathcal{C}_1` and :math:`\mathcal{C}_2`. On this theoretical example, the domain :math:`\mathbb{X}` is now a subset of a box :math:`[\mathbf{x}]`, easily representable and contractible. + +Codac already provides a catalog a contractors that one can use to deal with many applications. New contractors can also be designed and embedded in this framework without difficulty. + +.. admonition:: Constraint programming *vs* probabilistic methods + + The approach does differ from probabilistic methods in regards of the nature of the estimated solution. Probabilistic methods compute punctual potential solutions – for instance a vector – while in set-membership ones, it is the **set of all feasible solutions** that is evaluated, and thus an infinity of potential solutions. + + Another main distinction lies in the way things are computed: with set-membership methods, estimations are not randomly performed. **Computations are deterministic**: given a set of parameters or inputs, algorithms will always output the same result. + + +Reliable outputs +---------------- + +One of the advantages of this set-membership approach is the reliable outputs that are obtained. +By *reliable*, we mean that all sources of uncertainties are taken into account, including: + +- model parameter uncertainties +- measurement noise +- uncertainties related to time discretization +- linearization truncations +- approximation of real numbers by floating-point values + +The outcomes are intervals and tubes that are guaranteed to contain the solutions of the system. +This is well suited for proof purposes as we always consider worst-case possibilities when delineating the boundaries of the solution sets. + +The main drawback however, is that we may obtain large sets that may not be useful to characterize the solutions of the problem. We call this *pessimism*. This can be overcome by reformulating some constraints or by using bisections on sets. \ No newline at end of file diff --git a/doc/manual/manual/tools/index.rst b/doc/manual/manual/tools/index.rst index 065f21701..88bc066a1 100644 --- a/doc/manual/manual/tools/index.rst +++ b/doc/manual/manual/tools/index.rst @@ -7,4 +7,5 @@ Tools :maxdepth: 1 serialization.rst - registration.rst \ No newline at end of file + registration.rst + octasym.rst \ No newline at end of file diff --git a/doc/manual/manual/actions/octasym.rst b/doc/manual/manual/tools/octasym.rst similarity index 75% rename from doc/manual/manual/actions/octasym.rst rename to doc/manual/manual/tools/octasym.rst index c8dcc4383..d9d33b8a9 100644 --- a/doc/manual/manual/actions/octasym.rst +++ b/doc/manual/manual/tools/octasym.rst @@ -1,4 +1,4 @@ -.. _sec-actions-octasym: +.. _sec-tools-octasym: Octahedral symmetries ===================== diff --git a/src/core/contractors/codac2_CtcInverseNotIn.h b/src/core/contractors/codac2_CtcInverseNotIn.h index 4f715b90e..3f18a7209 100644 --- a/src/core/contractors/codac2_CtcInverseNotIn.h +++ b/src/core/contractors/codac2_CtcInverseNotIn.h @@ -59,4 +59,69 @@ namespace codac2 CtcUnion::contract(x); } }; + + // Template deduction guides + // Same as CtcInverse + + // ScalarType + + CtcInverseNotIn(const AnalyticFunction&, std::initializer_list, bool = true, bool = false) -> + CtcInverseNotIn; + + template + CtcInverseNotIn(const AnalyticFunction&, std::initializer_list, bool = true, bool = false) -> + CtcInverseNotIn; + + CtcInverseNotIn(const AnalyticFunction&, const Interval&, bool = true, bool = false) -> + CtcInverseNotIn; + + CtcInverseNotIn(const AnalyticFunction&, double, bool = true, bool = false) -> + CtcInverseNotIn; + + template + requires IsCtcBaseOrPtr + CtcInverseNotIn(const AnalyticFunction&, const C&, bool = true, bool = false) -> + CtcInverseNotIn; + + // VectorType + + CtcInverseNotIn(const AnalyticFunction&, std::initializer_list, bool = true, bool = false) -> + CtcInverseNotIn; + + CtcInverseNotIn(const AnalyticFunction&, std::initializer_list>, bool = true, bool = false) -> + CtcInverseNotIn; + + CtcInverseNotIn(const AnalyticFunction&, const Vector&, bool = true, bool = false) -> + CtcInverseNotIn; + + CtcInverseNotIn(const AnalyticFunction&, const IntervalVector&, bool = true, bool = false) -> + CtcInverseNotIn; + + template + requires (OtherDerived::RowsAtCompileTime == -1 && OtherDerived::ColsAtCompileTime == 1) + CtcInverseNotIn(const AnalyticFunction&, const Eigen::MatrixBase&, bool = true, bool = false) -> + CtcInverseNotIn; + + template + requires IsCtcBaseOrPtr + CtcInverseNotIn(const AnalyticFunction&, const C&, bool = true, bool = false) -> + CtcInverseNotIn; + + // MatrixType + + CtcInverseNotIn(const AnalyticFunction&, std::initializer_list>, bool = true, bool = false) -> + CtcInverseNotIn; + + CtcInverseNotIn(const AnalyticFunction&, std::initializer_list>>, bool = true, bool = false) -> + CtcInverseNotIn; + + template + requires (OtherDerived::RowsAtCompileTime == -1 && OtherDerived::ColsAtCompileTime == -1) + CtcInverseNotIn(const AnalyticFunction&, const Eigen::MatrixBase&, bool = true, bool = false) -> + CtcInverseNotIn; + + template + requires IsCtcBaseOrPtr + CtcInverseNotIn(const AnalyticFunction&, const C&, bool = true, bool = false) -> + CtcInverseNotIn; } \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5185c4be8..38b8f5210 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -42,6 +42,7 @@ list(APPEND SRC_TESTS # listing files without extension core/contractors/codac2_tests_CtcSegment core/contractors/codac2_tests_linear_ctc ../doc/manual/manual/contractors/geometric/src + ../doc/manual/manual/contractors/analytic/src core/domains/codac2_tests_BoolInterval core/domains/ellipsoid/codac2_tests_Ellipsoid diff --git a/tests/core/matrices/codac2_tests_cart_prod.cpp b/tests/core/matrices/codac2_tests_cart_prod.cpp index 28d5b1235..b152aaffd 100644 --- a/tests/core/matrices/codac2_tests_cart_prod.cpp +++ b/tests/core/matrices/codac2_tests_cart_prod.cpp @@ -49,4 +49,5 @@ TEST_CASE("cart_prod IntervalVector") CHECK(cart_prod(IntervalVector::empty(3)) == IntervalVector::empty(3)); CHECK(cart_prod(IntervalVector({{0,1},{2,3},{4,5}}),IntervalVector({{8,9}})) == IntervalVector({{0,1},{2,3},{4,5},{8,9}})); CHECK(cart_prod(25.,IntervalVector({{0,1},{2,3},{4,5}}),IntervalVector({{8,9}}),Vector::ones(3)) == IntervalVector({{25},{0,1},{2,3},{4,5},{8,9},{1},{1},{1}})); + CHECK(cart_prod(Interval(),Interval(),42.) == IntervalVector({{-oo,oo},{-oo,oo},{42.,42.}})); // cart_prod({-oo,oo},{-oo,oo},42.) not supported } \ No newline at end of file diff --git a/tests/core/matrices/codac2_tests_cart_prod.py b/tests/core/matrices/codac2_tests_cart_prod.py index 8e6f136c0..902e312d4 100644 --- a/tests/core/matrices/codac2_tests_cart_prod.py +++ b/tests/core/matrices/codac2_tests_cart_prod.py @@ -32,6 +32,7 @@ def test_cart_prod_intervalvector(self): self.assertTrue(cart_prod(IntervalVector.empty(3)) == IntervalVector.empty(3)) self.assertTrue(cart_prod([[0,1],[2,3],[4,5]],[[8,9]]) == IntervalVector([[0,1],[2,3],[4,5],[8,9]])) self.assertTrue(cart_prod(25.,[[0,1],[2,3],[4,5]],[[8,9]],Vector.ones(3)) == IntervalVector([[25],[0,1],[2,3],[4,5],[8,9],[1],[1],[1]])) + self.assertTrue(cart_prod([-oo,oo],[-oo,oo],42.) == IntervalVector([[-oo,oo],[-oo,oo],[42.,42.]])) if __name__ == '__main__': unittest.main() \ No newline at end of file