diff --git a/doc/manual/manual/contractors/index.rst b/doc/manual/manual/contractors/index.rst index 2f3dea786..e86d51bdc 100644 --- a/doc/manual/manual/contractors/index.rst +++ b/doc/manual/manual/contractors/index.rst @@ -1,3 +1,4 @@ +.. _sec-ctc: Contractors, separators ======================= diff --git a/doc/manual/manual/geometry/zonotope.png b/doc/manual/manual/geometry/zonotope.png new file mode 100644 index 000000000..fb74b6718 Binary files /dev/null and b/doc/manual/manual/geometry/zonotope.png differ diff --git a/doc/manual/manual/geometry/zonotope.rst b/doc/manual/manual/geometry/zonotope.rst index 8e6c54d24..2fe1f1354 100644 --- a/doc/manual/manual/geometry/zonotope.rst +++ b/doc/manual/manual/geometry/zonotope.rst @@ -12,15 +12,69 @@ The following classes represent zonotopes. Zonotope -------- -.. doxygenclass:: codac2::Zonotope - :project: codac - :members: +A zonotope is a convex and symmetric polytope. +It can be represented as the Minkowski sum of a finite number of line segments, which are called its generators. + +In Codac, zonotopes are represented by the class :class:`Zonotope`. A zonotope is defined by its center and its generators. +The center is a :class:`Vector`, noted :math:`z`, and the generators are stored in a :class:`Matrix`, noted :math:`A`. +Each column of the matrix corresponds to a generator. + +The resulting zonotope is: + +.. math:: + Z = z + A \cdot [-1,1]^m + +Where :math:`m` is the number of generators (i.e., the number of columns of the matrix :math:`A`). + +For example, let us consider the following zonotope in :math:`\mathbb{R}^2`: + +.. math:: + Z = \left(\begin{array}{c} 1\\ 2\\ \end{array}\right) + \left(\begin{array}{ccc} 1 & 0 & 2 \\ 0 & 1 & 1 \end{array}\right) \cdot [-1,1]^3 + +It can be constructed in Codac as follows: + +.. tabs:: + + .. code-tab:: py + + z = Vector([1,2]) + A = Matrix([[1,0,2],[0,1,1]]) + + Z = Zonotope(z, A) + + .. code-tab:: c++ + + Vector z({1,2}); + Matrix A({{1,0,2},{0,1,1}}); + + Zonotope Z(z, A); + + .. code-tab:: matlab + + z = Vector({1,2}); + A = Matrix({{1,0,2},{0,1,1}}); + + Z = Zonotope(z,A); + +The resulting zonotope is represented in the figure below: + +.. figure:: zonotope.png + :width: 400px + +Additional methods are provided to handle zonotopes. + +- `box()`: gives the bounding box of the zonotope. Returns an :class:`IntervalVector`. +- `proj(v)`: Projects the zonotope on a given subspace, defined by the vector of indices `v` (:class:`std::vector\`). Returns a :class:`Zonotope`. .. _subsec-zonotope-parallelepiped: Parallelepiped -------------- -.. doxygenclass:: codac2::Parallelepiped - :project: codac - :members: \ No newline at end of file +A parallelepiped is a special case of zonotope, where the number of generators is equal or less than the dimension of the space. + +It inherits from the `Zonotope` class, and thus has the same properties and methods. In addition, it defines the following methods: + +- `vertices()`: gives the vertices of the parallelepiped. Returns a :class:`std::vector\`. +- `contains(v)`: checks if the point `v` is contained in the parallelepiped. Returns a :class:`BoolInterval`. **The matrix A must be invertible** +- `is_superset(x)`: checks if the parallelepiped is a superset of the :class:`IntervalVector` `x`. Returns a :class:`BoolInterval`. **The matrix A must be invertible** \ No newline at end of file diff --git a/doc/manual/manual/geometry/zonotope.xml b/doc/manual/manual/geometry/zonotope.xml new file mode 100644 index 000000000..c64635b6f --- /dev/null +++ b/doc/manual/manual/geometry/zonotope.xml @@ -0,0 +1,422 @@ + + + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + +0.6 0 0 0.6 0 0 e + + + + + +0.5 0 0 0.5 0 0 e + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h + + + + + +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + + + +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h + + +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +-1 0.333 m +0 0 l +-1 -0.333 l + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.3 0 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.3 0 l +-0.5 -0.333 l +h + + + + +1 0 m +0 0.333 l +0 -0.333 l +h +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +1 0 m +0 0.333 l +0 -0.333 l +h +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +317.961 419.094 m +75.2427 257.282 l +75.2427 95.4693 l +196.602 95.4693 l +439.32 257.282 l +439.32 419.094 l +h + + +257.282 257.282 m +317.961 257.282 l + + +257.282 257.282 m +257.282 338.188 l + + +257.282 257.282 m +378.641 338.188 l + + +14.5631 14.5631 m +500 14.5631 l + + +14.5631 14.5631 m +14.5631 500 l + + +75.2427 4.85437 m +75.2427 14.5631 l + +-2 + +196.602 4.85437 m +196.602 14.5631 l + +0 + +317.961 4.85437 m +317.961 14.5631 l + +2 + +439.32 4.85437 m +439.32 14.5631 l + +4 + +4.85437 14.5631 m +14.5631 14.5631 l + +-1 + +4.85437 95.4693 m +14.5631 95.4693 l + +0 + +4.85437 176.375 m +14.5631 176.375 l + +1 + +4.85437 257.282 m +14.5631 257.282 l + +2 + +4.85437 338.188 m +14.5631 338.188 l + +3 + +4.85437 419.094 m +14.5631 419.094 l + +4 + +4.85437 500 m +14.5631 500 l + +5 +a_1 +a_2 +a_3 +Z + +z + + diff --git a/doc/manual/manual/installation/cpp.rst b/doc/manual/manual/installation/cpp.rst index c8f020144..1b97da510 100644 --- a/doc/manual/manual/installation/cpp.rst +++ b/doc/manual/manual/installation/cpp.rst @@ -15,39 +15,39 @@ Linux Installation ------------------ -Install from packages (latest release) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +.. Install from packages (latest release) +.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -This installation procedure is valid for Ubuntu (amd64, arm64), Debian (arm64, armhf) and possibly others. -A Debian package is available for the last release |version| of the library: +.. This installation procedure is valid for Ubuntu (amd64, arm64), Debian (arm64, armhf) and possibly others. +.. A Debian package is available for the last release |version| of the library: -.. code-block:: bash +.. .. code-block:: bash - sudo sh -c 'echo "deb [trusted=yes] https://webperso.ensta.fr/packages/$(if [ -z "$(. /etc/os-release && echo $UBUNTU_CODENAME)" ]; then echo debian/$(. /etc/os-release && echo $VERSION_CODENAME); else echo ubuntu/$(. /etc/os-release && echo $UBUNTU_CODENAME); fi) ./" > /etc/apt/sources.list.d/ensta-bretagne.list' - sudo apt update - sudo apt install libcodac-dev +.. sudo sh -c 'echo "deb [trusted=yes] https://webperso.ensta.fr/packages/$(if [ -z "$(. /etc/os-release && echo $UBUNTU_CODENAME)" ]; then echo debian/$(. /etc/os-release && echo $VERSION_CODENAME); else echo ubuntu/$(. /etc/os-release && echo $UBUNTU_CODENAME); fi) ./" > /etc/apt/sources.list.d/ensta-bretagne.list' +.. sudo apt update +.. sudo apt install libcodac-dev -Then, check your installation :ref:`with the instructions of this page `. +.. Then, check your installation :ref:`with the instructions of this page `. -.. admonition:: Uninstall Codac +.. .. admonition:: Uninstall Codac - To uninstall Codac, you might want to do the following: +.. To uninstall Codac, you might want to do the following: - .. code-block:: bash +.. .. code-block:: bash - sudo apt remove libcodac-dev libibex-dev - sudo rm -f /etc/apt/sources.list.d/ensta-bretagne.list - sudo apt update +.. sudo apt remove libcodac-dev libibex-dev +.. sudo rm -f /etc/apt/sources.list.d/ensta-bretagne.list +.. sudo apt update -.. admonition:: Standalone archives +.. .. admonition:: Standalone archives - Standalone archives exist also for all the supported configurations, *e.g.* for a Raspberry Pi running Raspberry Pi OS Bookworm 32 bit, download and extract ``codac_standalone_armhf_bookworm.zip`` from ``_, then in the ``example`` folder run: +.. Standalone archives exist also for all the supported configurations, *e.g.* for a Raspberry Pi running Raspberry Pi OS Bookworm 32 bit, download and extract ``codac_standalone_armhf_bookworm.zip`` from ``_, then in the ``example`` folder run: - .. code-block:: bash +.. .. code-block:: bash - cmake . ; cmake --build . ; ./my_project +.. cmake . ; cmake --build . ; ./my_project - and check that the graphical output appears. +.. and check that the graphical output appears. Install from the sources (latest developments) diff --git a/doc/manual/manual/tools/octasym.rst b/doc/manual/manual/tools/octasym.rst index d9d33b8a9..374a50f65 100644 --- a/doc/manual/manual/tools/octasym.rst +++ b/doc/manual/manual/tools/octasym.rst @@ -3,4 +3,149 @@ Octahedral symmetries ===================== -Further documentation upcoming. \ No newline at end of file +In Codac, the class :class:`OctaSym` is used to represent the symmetries of the hypercube, which are elements of the hyperoctahedral group. This class provides methods to create, manipulate, and apply these symmetries to objects. + +Hyperoctahedral group +--------------------- + +The hyperoctahedral :math:`B_n` group is the group of symmetries of the hypercube :math:`[-1,1]^n`. It contains :math:`2^n \cdot n!` elements. + +It is a signed symmetric group of permutations. It means that each symmetry of the hypercube can be represented by a +permutation of the coordinates, and a sign for each coordinate (i.e., a reflection or not). + +Two main representations of the hyperoctahedral group exist: + +- The **Cauchy's representation**: It is a two-row representation where the first row contains the indices of the coordinates, and the second row contains the new indices after the permutation, with a sign indicating whether there is a reflection or not. + +- The **Permutation representation**: It is a :math:`n\times n` matrix representation where each row and column has exactly one non-zero entry, which is either +1 or -1, indicating the permutation and reflection. + +For example, in dimension three the same symmetry :math:`\sigma` can be written in Cauchy's representation as: + +.. math:: + + \sigma = \begin{pmatrix} + 1 & 2 & 3 \\ + 3 & -1 & 2 + \end{pmatrix} + +And in permutation representation as: + +.. math:: + + P_{\sigma} = \begin{pmatrix} + 0 & 0 & 1 \\ + -1 & 0 & 0 \\ + 0 & 1 & 0 + \end{pmatrix} + +If we apply the symmetry :math:`\sigma` to a point :math:`\mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}`, we get: + +.. math:: + + \sigma(\mathbf{x}) = P_{\sigma}\cdot \mathbf{x} = \begin{pmatrix} x_3 \\ -x_1 \\ x_2 \end{pmatrix} + + +Codac implementation +-------------------- + +The main representation used in Codac is the Cauchy's notation. **We suppose that the first row of the cauchy's representation is in increasing order, i.e. (1,2,...n)**. + +The constructors of the :class:`OctaSym` class then require a single vector reprensenting the second row of the Cauchy's representation. + +For example, the symmetry :math:`\sigma` mentionned earlier can be created in Codac as follows: + +.. tabs:: + + .. code-tab:: py + + sigma = OctaSym([3,-1,2]) + + .. code-tab:: c++ + + OctaSym sigma ({3,-1,2}); + + .. code-tab:: matlab + + sigma = OctaSym({3,-1,2}); + +The permutation matrix can be recovered using the method `permutation_matrix()`. It returns the matrix :math:`P_{\sigma}` as a :class:`Matrix` object + +.. tabs:: + + .. code-tab:: py + + P_sigma = sigma.permutation_matrix() + + .. code-tab:: c++ + + Matrix P_sigma = sigma.permutation_matrix(); + + .. code-tab:: matlab + + P_sigma = sigma.permutation_matrix(); + +The inverse of a symmetry can be computed using the method `invert()`. These symmetries can also be composed using the multiplication operator `*`. + +.. tabs:: + + .. code-tab:: py + + sigma_inv = sigma.invert() + sigma_comp = sigma * sigma_inv # identity symmetry + + .. code-tab:: c++ + + OctaSym sigma_inv = sigma.invert(); + OctaSym sigma_comp = sigma * sigma_inv; // identity symmetry + + .. code-tab:: matlab + + sigma_inv = sigma.invert(); + sigma_comp = sigma * sigma_inv; % identity symmetry + +To apply the symmetry to an object, the classical operator `()` is used. Objects which support the application of a symmetry are listed below. +**Objects in bold are not supported yet but will be in the future**. + +- Vectors + + - :class:`Vector`. :ref:`see more ` + - :class:`IntervalVector`. :ref:`see more ` + - :class:`VectorVar` and :class:`VectorExpr` (for functions). :ref:`see more ` + +- Contractors and separators. :ref:`see more ` + +- **Analytic** and sampled trajectories. :class:`SampledTraj`. + +- **Sliced tubes**. + + +For example, wtih the symmetry :math:`\sigma` defined earlier: + +.. tabs:: + + .. code-tab:: py + + X = Vector([4,5,6]) + X_s = sigma(X) # X_s is the (poncual) interval vector (6,-4,5) + + x = VectorVar(3) + f = AnalyticFunction([x], sigma(2*x)) + Y = f.eval(X) # Y is the (poncual) interval vector (12,-8,10) + + .. code-tab:: c++ + + Vector X ({4,5,6}); + Vector X_s = sigma(X); // X_s is (poncual) interval the vector (6,-4,5) + + VectorVar x(3); + AnalyticFunction f({x}, sigma(2*x)); + auto Y = f.eval(X); // Y is the (poncual) interval vector (12,-8,10) + + .. code-tab:: matlab + + X = Vector({4,5,6}); + X_s = sigma(X); % X_s is the (poncual) interval vector (6,-4,5) + + x = VectorVar(3); + f = AnalyticFunction({x}, sigma(2*x)); + Y = f.eval(X); % Y is the (poncual) interval vector (12,-8,10) \ No newline at end of file diff --git a/src/core/actions/codac2_OctaSym.h b/src/core/actions/codac2_OctaSym.h index 46a9f6037..81cdb1d7c 100644 --- a/src/core/actions/codac2_OctaSym.h +++ b/src/core/actions/codac2_OctaSym.h @@ -33,20 +33,56 @@ namespace codac2 /** * \class OctaSym + * \brief Represents an hyperoctahedral symmetry. */ class OctaSym : public std::vector, public Action { public: + /** + * \brief Constructs an hyperoctahedral symmetry from a list of integers. The list represents the second line in the Cauchy' representation. + * It is supposed that the first line of the representation is (1 2 3 ... n). + * + * \param s The list of integers representing the hyperoctahedral symmetry. + */ OctaSym(std::initializer_list s); + + /** + * \brief Constructs an hyperoctahedral symmetry from a vector of integers. The vector represents the second line in the Cauchy' representation. + * It is supposed that the first line of the representation is (1 2 3 ... n). + * + * \param s The vector of integers representing the hyperoctahedral symmetry. + */ OctaSym(const std::vector& s); + /** + * \brief Inverts of the hyperoctahedral symmetry. + * + * \return The inverse of the hyperoctahedral symmetry. + */ OctaSym invert() const; + /** + * \brief Composes the hyperoctahedral symmetry with another one. + * + * \param s The hyperoctahedral symmetry to compose with. + * \return The composition of the two hyperoctahedral symmetries. + */ OctaSym operator*(const OctaSym& s) const; + /** + * \brief Computes the permutation matrix associated to the hyperoctahedral symmetry. + * + * \return The permutation matrix associated to the hyperoctahedral symmetry. + */ Matrix permutation_matrix() const; + /** + * \brief Applies the hyperoctahedral symmetry to a vector. + * + * \param x The Vector or IntervalVector to which the hyperoctahedral symmetry is applied. + * \return The result of the application of the hyperoctahedral symmetry to the vector. + */ template requires (Derived::ColsAtCompileTime == 1) Mat operator()(const Eigen::MatrixBase& x) const @@ -58,16 +94,34 @@ namespace codac2 return x_; } + /** + * \brief Applies the hyperoctahedral symmetry to a Contractor. + * + * \param c The Contractor to which the hyperoctahedral symmetry is applied. + * \return The result of the application of the hyperoctahedral symmetry to the Contractor. + */ template requires IsCtcBaseOrPtr CtcAction operator()(const C& c) const; // -> is defined in CtcAction class + /** + * \brief Applies the hyperoctahedral symmetry to a Separator. + * + * \param s The Separator to which the hyperoctahedral symmetry is applied. + * \return The result of the application of the hyperoctahedral symmetry to the Separator. + */ template requires is_sep_v SepAction operator()(const S& s) const; // -> is defined in SepAction class + /** + * \brief Applies the hyperoctahedral symmetry to a SetExpr. + * + * \param x1 The SetExpr to which the hyperoctahedral symmetry is applied. + * \return The result of the application of the hyperoctahedral symmetry to the SetExpr. + */ std::shared_ptr operator()(const std::shared_ptr& x1) const; // -> is defined in set operations file @@ -80,6 +134,12 @@ namespace codac2 return y; } + /** + * \brief Applies the hyperoctahedral symmetry to a VectorExpr or a VectorVar. + * + * \param x1 The VectorExpr or VectorVar to which the hyperoctahedral symmetry is applied. + * \return The result of the application of the hyperoctahedral symmetry to the VectorExpr or VectorVar. + */ template // To avoid ambiguity with operator()(const Eigen::MatrixBase& x): requires (std::is_same_v || std::is_same_v) @@ -92,6 +152,13 @@ namespace codac2 return { std::make_shared>(*this, x1) }; } + /** + * \brief Overloads the stream insertion operator to print the hyperoctahedral symmetry in a human-readable format. + * + * \param str The output stream to which the hyperoctahedral symmetry is printed. + * \param s The hyperoctahedral symmetry to print. + * \return The output stream to which the hyperoctahedral symmetry is printed. + */ friend std::ostream& operator<<(std::ostream& str, const OctaSym& s) { str << "("; diff --git a/src/core/contractors/codac2_CtcLohner.cpp b/src/core/contractors/codac2_CtcLohner.cpp new file mode 100644 index 000000000..1d53fb1dc --- /dev/null +++ b/src/core/contractors/codac2_CtcLohner.cpp @@ -0,0 +1,22 @@ +/** + * codac2_CtcLohner.cpp + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Simon Rohou + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_CtcLohner.h" +#include "codac2_Slice.h" +#include "codac2_Interval.h" + +using namespace std; +using namespace codac2; + +CtcLohner::CtcLohner(const Function &f, int contractions, double eps) + : Ctc(), + m_f(f), + contractions(contractions), + dim(f.nb_var()), + eps(eps) {} diff --git a/src/core/contractors/codac2_CtcLohner.h b/src/core/contractors/codac2_CtcLohner.h new file mode 100644 index 000000000..3037bddbb --- /dev/null +++ b/src/core/contractors/codac2_CtcLohner.h @@ -0,0 +1,41 @@ +/** + * \file codac2_CtcLohner.h + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Simon Rohou + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include "codac2_Ctc.h" +#include "codac2_AnalyticFunction.h" +#include "codac2_TimePropag.h" +#include "codac2_TDomain.h" + +namespace codac2 +{ + template + class Slice; + + template + class SlicedTube; + + class CtcLohner + { + public: + + CtcLohner(const AnalyticFunction& f, int contractions = 5, double eps = 0.1); + + protected: + + AnalyticFunction m_f; //!< forward function + int contractions; //!< number of contractions of the global enclosure by the estimated local enclosure + int dim; //!< dimension of the state vector + double eps; //!< inflation parameter for the global enclosure + + static const std::string m_ctc_name; //!< class name (mainly used for CN Exceptions) + static std::vector m_str_expected_doms; + }; +} \ No newline at end of file