diff --git a/CMakeLists.txt b/CMakeLists.txt index 7eae26a62..f65251466 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -144,6 +144,11 @@ endif() # Adds Eigen3::Eigen +################################################################################ +# Looking for Threads +################################################################################ + + find_package(Threads REQUIRED) ################################################################################ # Looking for CAPD (if needed) diff --git a/doc/manual/manual/extensions/capd/capd.rst b/doc/manual/manual/extensions/capd/capd.rst new file mode 100644 index 000000000..70db69213 --- /dev/null +++ b/doc/manual/manual/extensions/capd/capd.rst @@ -0,0 +1,203 @@ +.. _sec-extensions-capd-capd: + +CAPD (rigorous numerics in dynamical systems) +============================================= + + Main author: `Maël Godard `_ + +This page describes how to use the CAPD library with Codac. CAPD is a C++ library for rigorous numerics in dynamical systems. + +To use CAPD with Codac, you first need to install the CAPD library. You can find the installation instructions on the `CAPD website `_. + +Note that as CAPD is a C++ only library, the content present in this page is **only available in C++**. + + +.. _subsec-extensions-capd-capd-install: +Installing the ``codac-capd`` extension +--------------------------------------- + +To install the ``codac-capd`` extension, you need to install the Codac library from its sources. This can be done :ref:`by using CMake ` with the option ``WITH_CAPD=ON``. For example: + +.. code-block:: bash + + cmake -DCMAKE_INSTALL_PREFIX=$HOME/ibex-lib/build_install -DCMAKE_BUILD_TYPE=Release -DWITH_CAPD=ON .. + +We highly recommend to test the installation of the library with the provided tests. To do so, you can use the following command: + +.. code-block:: bash + + make test + + +Content +------- + +The ``codac-capd`` extension provides functions to convert CAPD objects to Codac objects and vice versa. + +The functions are ``to_capd`` and ``to_codac``. They can be used to convert the following objects: + +- ``capd::Interval`` :math:`\leftrightarrow` ``codac2::Interval`` +- ``capd::IVector`` :math:`\leftrightarrow` ``codac2::IntervalVector`` +- ``capd::IMatrix`` :math:`\leftrightarrow` ``codac2::IntervalMatrix`` +- ``capd::ITimeMap::SolutionCurve`` :math:`\rightarrow` ``codac2::SlicedTube`` + + +How to use +---------- + +The header of the ``codac-capd`` extension is not included by default. You need to include it manually in your code, together with the CAPD library: + +.. code-block:: c++ + + #include + #include + +You can use the functions ``to_capd`` and ``to_codac`` to convert between CAPD and Codac objects as follows: + +.. tabs:: + + .. code-tab:: c++ + + codac2::Interval codac_interval(0,2); // Codac interval [0, 2] + capd::Interval capd_interval = to_capd(codac_interval); // convert to CAPD interval + codac2::Interval codac_interval2 = to_codac(capd_interval); // convert back to Codac interval + + +Example +------- + +.. image:: img/pendulum.png + :alt: State variables of the pendulum + :align: right + :width: 130px + +For this example we will consider the pendulum with friction. + +The state variables of the pendulum are its angle :math:`\theta` and its angular velocity :math:`\omega`. The pendulum follows the following dynamic: + +.. math:: + + \left(\begin{array}{c} + \dot{\theta}\\ + \dot{\omega} + \end{array}\right)=\left(\begin{array}{c} + \omega\\ + -\sin(\theta)\cdot\frac{g}{l}-0.5\omega + \end{array}\right), + + +where :math:`g` is the gravity constant and :math:`l` is the length of the pendulum. + +This equation can be passed to the CAPD library as follows: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-2-beg] + :end-before: [codac-capd-2-end] + :dedent: 2 + +To solve this ODE, an ``IOdeSolver`` object is necessary. + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-3-beg] + :end-before: [codac-capd-3-end] + :dedent: 2 + +CAPD then uses an ``ITimeMap`` to make the link between a time step and the solution of the ODE at this time. The ``I`` here stands for ``Interval`` as the solution is an interval guaranteed to enclose the solution. Here we will integrate the ODE between :math:`t_0=0s` and :math:`t_f=20s`. + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-4-beg] + :end-before: [codac-capd-4-end] + :dedent: 2 + +To completly define the ODE, we need to define the initial conditions. Here we will set the initial angle to :math:`\theta_0=-\frac{\pi}{2}` and the +initial angular velocity to :math:`\omega_0=0`. For the purpose of this example, we will add a small uncertainty to the initial conditions. The initial conditions are then defined as follows: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-5-beg] + :end-before: [codac-capd-5-end] + :dedent: 2 + +There are then two ways to get the result of the integration depending on the use case. + +If the desired result is the solution of the ODE at a given time (here say :math:`T=1s`), we can do as follows: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-6-beg] + :end-before: [codac-capd-6-end] + :dedent: 2 + +**Be careful, this method modifies the initial set** ``s`` **in place**. + +If the desired result is the solution curve (or tube) of the ODE on the time domain :math:`[t_0,t_f]`, we can do as follows: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-7-beg] + :end-before: [codac-capd-7-end] + :dedent: 2 + +The variable ``solution`` is the desired solution curve (or tube). The operator ``solution(t)`` gives the solution at time :math:`t`. +It can be converted into a Codac ``SlicedTube`` with the function ``to_codac``. This functions takes two arguments: + +- the ``capd::ITimeMap::SolutionCurve`` to convert. +- a ``codac2::TDomain`` object defining the temporal domain of the tube. + +The resulting ``SlicedTube`` will have the same time domain as the one given in argument, completed with the CAPD gates. An example of conversion is : + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-8-beg] + :end-before: [codac-capd-8-end] + :dedent: 2 + +A full display can be done with the following code: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [codac-capd-9-beg] + :end-before: [codac-capd-9-end] + :dedent: 4 + +The result is the following figure, with in green the initial set (:math:`t=0s`) and in red the final set (:math:`t=20s`). The ``SlicedTube`` is displayed in blue with a black edge for better visibility. The orange rectangles correspond to the gates (degenerate slices). + +.. figure:: img/pendulum_result.png + :width: 500px + + Result of the CAPD integration of the pendulum, enclosed in a Codac tube. diff --git a/doc/manual/manual/extensions/capd/img/lorenz.png b/doc/manual/manual/extensions/capd/img/lorenz.png new file mode 100644 index 000000000..cc954a17c Binary files /dev/null and b/doc/manual/manual/extensions/capd/img/lorenz.png differ diff --git a/doc/manual/manual/extensions/capd/img/lorenz_0.png b/doc/manual/manual/extensions/capd/img/lorenz_0.png new file mode 100644 index 000000000..c59205b2e Binary files /dev/null and b/doc/manual/manual/extensions/capd/img/lorenz_0.png differ diff --git a/doc/manual/manual/extensions/capd/img/lorenz_0_05.png b/doc/manual/manual/extensions/capd/img/lorenz_0_05.png new file mode 100644 index 000000000..689927e13 Binary files /dev/null and b/doc/manual/manual/extensions/capd/img/lorenz_0_05.png differ diff --git a/doc/manual/manual/extensions/capd/img/lorenz_0_1.png b/doc/manual/manual/extensions/capd/img/lorenz_0_1.png new file mode 100644 index 000000000..78ff78869 Binary files /dev/null and b/doc/manual/manual/extensions/capd/img/lorenz_0_1.png differ diff --git a/doc/manual/manual/extensions/capd/img/pendulum_peibos.png b/doc/manual/manual/extensions/capd/img/pendulum_peibos.png new file mode 100644 index 000000000..927185ce1 Binary files /dev/null and b/doc/manual/manual/extensions/capd/img/pendulum_peibos.png differ diff --git a/doc/manual/manual/extensions/capd/index.rst b/doc/manual/manual/extensions/capd/index.rst index 3bb0b4a0c..1a683ee7b 100644 --- a/doc/manual/manual/extensions/capd/index.rst +++ b/doc/manual/manual/extensions/capd/index.rst @@ -1,202 +1,9 @@ .. _sec-extensions-capd: -CAPD (rigorous numerics in dynamical systems) -============================================= +CAPD Extension +============== - Main author: `Maël Godard `_ +.. toctree:: -This page describes how to use the CAPD library with Codac. CAPD is a C++ library for rigorous numerics in dynamical systems. - -To use CAPD with Codac, you first need to install the CAPD library. You can find the installation instructions on the `CAPD website `_. - -Note that as CAPD is a C++ only library, the content present in this page is **only available in C++**. - - -Installing the ``codac-capd`` extension ---------------------------------------- - -To install the ``codac-capd`` extension, you need to install the Codac library from its sources. This can be done :ref:`by using CMake ` with the option ``WITH_CAPD=ON``. For example: - -.. code-block:: bash - - cmake -DCMAKE_INSTALL_PREFIX=$HOME/ibex-lib/build_install -DCMAKE_BUILD_TYPE=Release -DWITH_CAPD=ON .. - -We highly recommend to test the installation of the library with the provided tests. To do so, you can use the following command: - -.. code-block:: bash - - make test - - -Content -------- - -The ``codac-capd`` extension provides functions to convert CAPD objects to Codac objects and vice versa. - -The functions are ``to_capd`` and ``to_codac``. They can be used to convert the following objects: - -- ``capd::Interval`` :math:`\leftrightarrow` ``codac2::Interval`` -- ``capd::IVector`` :math:`\leftrightarrow` ``codac2::IntervalVector`` -- ``capd::IMatrix`` :math:`\leftrightarrow` ``codac2::IntervalMatrix`` -- ``capd::ITimeMap::SolutionCurve`` :math:`\rightarrow` ``codac2::SlicedTube`` - - -How to use ----------- - -The header of the ``codac-capd`` extension is not included by default. You need to include it manually in your code, together with the CAPD library: - -.. code-block:: c++ - - #include - #include - -You can use the functions ``to_capd`` and ``to_codac`` to convert between CAPD and Codac objects as follows: - -.. tabs:: - - .. code-tab:: c++ - - codac2::Interval codac_interval(0,2); // Codac interval [0, 2] - capd::Interval capd_interval = to_capd(codac_interval); // convert to CAPD interval - codac2::Interval codac_interval2 = to_codac(capd_interval); // convert back to Codac interval - - -Example -------- - -.. image:: img/pendulum.png - :alt: State variables of the pendulum - :align: right - :width: 130px - -For this example we will consider the pendulum with friction. - -The state variables of the pendulum are its angle :math:`\theta` and its angular velocity :math:`\omega`. The pendulum follows the following dynamic: - -.. math:: - - \left(\begin{array}{c} - \dot{\theta}\\ - \dot{\omega} - \end{array}\right)=\left(\begin{array}{c} - \omega\\ - -\sin(\theta)\cdot\frac{g}{l}-0.5\omega - \end{array}\right), - - -where :math:`g` is the gravity constant and :math:`l` is the length of the pendulum. - -This equation can be passed to the CAPD library as follows: - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-2-beg] - :end-before: [codac-capd-2-end] - :dedent: 2 - -To solve this ODE, an ``IOdeSolver`` object is necessary. - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-3-beg] - :end-before: [codac-capd-3-end] - :dedent: 2 - -CAPD then uses an ``ITimeMap`` to make the link between a time step and the solution of the ODE at this time. The ``I`` here stands for ``Interval`` as the solution is an interval guaranteed to enclose the solution. Here we will integrate the ODE between :math:`t_0=0s` and :math:`t_f=20s`. - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-4-beg] - :end-before: [codac-capd-4-end] - :dedent: 2 - -To completly define the ODE, we need to define the initial conditions. Here we will set the initial angle to :math:`\theta_0=-\frac{\pi}{2}` and the -initial angular velocity to :math:`\omega_0=0`. For the purpose of this example, we will add a small uncertainty to the initial conditions. The initial conditions are then defined as follows: - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-5-beg] - :end-before: [codac-capd-5-end] - :dedent: 2 - -There are then two ways to get the result of the integration depending on the use case. - -If the desired result is the solution of the ODE at a given time (here say :math:`T=1s`), we can do as follows: - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-6-beg] - :end-before: [codac-capd-6-end] - :dedent: 2 - -**Be careful, this method modifies the initial set** ``s`` **in place**. - -If the desired result is the solution curve (or tube) of the ODE on the time domain :math:`[t_0,t_f]`, we can do as follows: - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-7-beg] - :end-before: [codac-capd-7-end] - :dedent: 2 - -The variable ``solution`` is the desired solution curve (or tube). The operator ``solution(t)`` gives the solution at time :math:`t`. -It can be converted into a Codac ``SlicedTube`` with the function ``to_codac``. This functions takes two arguments: - -- the ``capd::ITimeMap::SolutionCurve`` to convert. -- a ``codac2::TDomain`` object defining the temporal domain of the tube. - -The resulting ``SlicedTube`` will have the same time domain as the one given in argument, completed with the CAPD gates. An example of conversion is : - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-8-beg] - :end-before: [codac-capd-8-end] - :dedent: 2 - -A full display can be done with the following code: - -.. tabs:: - - .. group-tab:: C++ - - .. literalinclude:: src.cpp - :language: c++ - :start-after: [codac-capd-9-beg] - :end-before: [codac-capd-9-end] - :dedent: 4 - -The result is the following figure, with in green the initial set (:math:`t=0s`) and in red the final set (:math:`t=20s`). The ``SlicedTube`` is displayed in blue with a black edge for better visibility. The orange rectangles correspond to the gates (degenerate slices). - -.. figure:: img/pendulum_result.png - :width: 500px - - Result of the CAPD integration of the pendulum, enclosed in a Codac tube. \ No newline at end of file + capd.rst + peibos_capd.rst \ No newline at end of file diff --git a/doc/manual/manual/extensions/capd/peibos_capd.rst b/doc/manual/manual/extensions/capd/peibos_capd.rst new file mode 100644 index 000000000..4d5c92216 --- /dev/null +++ b/doc/manual/manual/extensions/capd/peibos_capd.rst @@ -0,0 +1,148 @@ +.. _sec-extensions-capd-peibos: + +PEIBOS-CAPD +=========== + +When compiling CODAC with the codac-capd extension (see :ref:`here `), the CAPD version of the PEIBOS library is also compiled. + +Let us consider an initial set :math:`\mathbb{X}_0 \subset \mathbb{R}^n` with its boundary :math:`\partial \mathbb{X}_0`. +Considering a dynamical system :math:`\dot{\mathbf{x}}=\gamma(\mathbf{x})`, the PEIBOS tool allows to compute the reach set :math:`\mathbf{X}_t=\left\{ \mathbf{x}(t) \mid \mathbf{x} \in \partial \mathbb{X}_0 \right\}`. + +Gnomonic atlas +-------------- + +To handle the boundary of the initial set :math:`\mathbb{X}_0`, the PEIBOS tool relies on a gnomonic atlas. See :ref:`subsec-functions-peibos-gnomonic-atals`. + +Use +--- + +The computation of the reach set is decomposed into two separate functions. + +PEIBOS +~~~~~~ + +The PEIBOS function takes at least six arguments : + +- The CAPD IMap representing :math:`\gamma`. +- The final time for the integration of the ODE. +- A timestep to get intermediate states. +- The inverse chart for the gnomonic atlas (an analytic function). +- The list of symmetries for the gnomonic atlas. Note that each symmetry is represented as a hyperoctahedral symmetry, see :ref:`sec-actions-octasym`. +- A resolution :math:`\epsilon`. The initial box :math:`\left[-1,1\right]^m` will initiallly be splitted in boxes with a diameter smaller than :math:`\epsilon`. +- Eventually an offset vector can be specified if the initial set is not centered around the origin. +- Eventually a flag can be set to True to get the verbose. + +For each of the small box :math:`\left[\mathbf{x}(0)\right]`, this function computes a box containing :math:`\bar{\mathbf{x}}(t)` and an interval matrix enclosing the Jacobian matrix :math:`D\mathbf{\left[x\right]}(t)`. + +It returns a map where the keys are the time steps. For each time step, the value is a vector of tuple, where each tuple contains: + +- A `PEIBOS_CAPD_Key` representing the symmetry and the initial box used (plus an eventual offset), i.e. :math:`\mathbf{x}(0)= \sigma(\psi_0(\mathbf{\text{box}})) + \text{offset}` +- A box enclosing :math:`\bar{\mathbf{x}}(t)` +- An interval matrix enclosing the Jacobian matrix :math:`D\mathbf{\left[x\right]}(t)` + +Each tuple can then be used to build a Parallelepiped enclosing :math:`\mathbf{x}(t)` using the parallelepiped inclusion. This is done in the `reach_set` function. + +The full signature of the function is : + +.. doxygenfunction:: codac2::PEIBOS(const capd::IMap&, double, double, const AnalyticFunction&, const std::vector&, double, const Vector&, bool); + :project: codac + +reach_set +~~~~~~~~~ + +The reach_set function takes only one argument : the map returned by the PEIBOS function. + +It returns a map where the keys are the time steps. For each time step, the value is a vector of :ref:`subsec-zonotope-parallelepiped`. Their union is an outer approximation of :math:`\mathbf{X}_t`. + +This function is then simply : + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [peibos-capd-1-beg] + :end-before: [peibos-capd-1-end] + :dedent: 0 + +The `parallelepiped_inclusion` is the one described in `this article `_. + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [peibos-capd-2-beg] + :end-before: [peibos-capd-2-end] + :dedent: 0 + +Examples +-------- + +2D : Pendulum +~~~~~~~~~~~~~ + +Say that we want to integrate the state of the pendulum starting from the an initial box. It is defined by + +.. math:: + \dot{x}=\gamma(x) = \begin{pmatrix} + x_2 \\ + -5\cdot\sin (x_1) - 0.5\cdot x_2 + \end{pmatrix} + +The corresponding code is: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [peibos-capd-3-beg] + :end-before: [peibos-capd-3-end] + :dedent: 4 + +The result is + +.. image:: img/pendulum_peibos.png + :alt: 20 seconds of integration of the pendulum + :align: center + :width: 400px + +3D : Lorenz system +~~~~~~~~~~~~~~~~~~ + +Say that we want to integrate the state of the pendulum starting from the an initial sphere. It is defined by + +.. math:: + \dot{x}=\gamma(x) = \begin{pmatrix} + \sigma (x_2 - x_1) \\ + x_1 (\rho - x_3) - x_2 \\ + x_1 x_2 - \beta x_3 + \end{pmatrix} + +The corresponding code is: + +.. tabs:: + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [peibos-capd-4-beg] + :end-before: [peibos-capd-4-end] + :dedent: 4 + +The result is + +.. image:: img/lorenz.png + :alt: Integration of the Lorenz system + :align: center + :width: 500px + +Related work +------------ + +This method comes from `this article `_. \ No newline at end of file diff --git a/doc/manual/manual/extensions/capd/src.cpp b/doc/manual/manual/extensions/capd/src.cpp index 2b4e21c9c..57d015666 100644 --- a/doc/manual/manual/extensions/capd/src.cpp +++ b/doc/manual/manual/extensions/capd/src.cpp @@ -8,9 +8,51 @@ #include // [codac-capd-1-end] +using namespace codac2; + +// [peibos-capd-2-beg] +Parallelepiped parallelepiped_inclusion_(const IntervalVector& Y, const IntervalMatrix& Jf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X) +{ + // Computation of the Jacobian of g = f o sigma(psi_0) + IntervalMatrix Jg = Jf * (sigma.permutation_matrix().template cast()) * psi_0.diff(X); + + Vector z = Y.mid(); + // A is an approximation of the Jacobian of g at the center of X + Matrix A = (Jf_tild * sigma.permutation_matrix() * (psi_0.diff(X.mid()).mid())); + + // Maximum error computation + double rho = error_peibos(Y, z, Jg, A, X); + + // Inflation of the parallelepiped + Matrix A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); + + return Parallelepiped(z, A_inf); +} +// [peibos-capd-2-end] + +// [peibos-capd-1-beg] +std::map> reach_set_(const std::map>>& peibos_output) +{ + std::map> output; + + for (const auto& [time,vec] : peibos_output) + { + for (const auto& [key, z, Jf] : vec) + { + auto p = parallelepiped_inclusion_(z, Jf, Jf.mid(), key.psi_0, key.sigma, key.box); + + output[time].push_back(p); + } + } + + return output; +} +// [peibos-capd-1-end] int main() { + set_threads_used(max_threads()); + // Equation of the pendulum with friction // [codac-capd-2-beg] capd::IMap vectorField("par:l,g;var:t,w;fun:w,-sin(t)*g/l - 0.5*w;"); @@ -59,7 +101,7 @@ int main() // [codac-capd-8-beg] - auto tdomain = create_tdomain(codac2::Interval(0,20),0.05, true); // true to have gates + auto tdomain = create_tdomain(Interval(0,20),0.05, true); // true to have gates auto codac_tube = to_codac(solution, tdomain); // [codac-capd-8-end] @@ -73,17 +115,51 @@ int main() std::cout << "\n\nafter time=" << T << " the image is: " << result; std::cout << "\ndiam(image): " << diam(result) << std::endl << std::endl; - codac2::DefaultFigure::set_axes(codac2::axis(0,{-2,1.5}),codac2::axis(1,{-2,3})); + DefaultFigure::set_axes(axis(0,{-2,1.5}),axis(1,{-2,3})); - codac2::DefaultFigure::draw_tube(codac_tube, codac2::ColorMap::blue_tube()); - codac2::DefaultFigure::draw_tube(codac_tube, codac2::Color::black()); + DefaultFigure::draw_tube(codac_tube, ColorMap::blue_tube()); + DefaultFigure::draw_tube(codac_tube, Color::black()); for (float t=0.;t<20.;t+=0.05) - codac2::DefaultFigure::draw_box(codac2::to_codac(solution(t)), {codac2::Color::none(), codac2::Color::orange(0.5)}); + DefaultFigure::draw_box(to_codac(solution(t)), {Color::none(), Color::orange(0.5)}); - codac2::DefaultFigure::draw_box(codac2::to_codac(c),codac2::Color::green()); - codac2::DefaultFigure::draw_box(codac2::to_codac(result),codac2::Color::red()); + DefaultFigure::draw_box(to_codac(c),Color::green()); + DefaultFigure::draw_box(to_codac(result),Color::red()); // [codac-capd-9-end] } + + { + // [peibos-capd-3-beg] + capd::IMap vectorField_pend("par:l,g;var:t,w;fun:w,-sin(t)*5 - 0.5*w;"); + + VectorVar X_2d(1); + AnalyticFunction psi0_pend ({X_2d},{0.1*X_2d[0],0.1}); + OctaSym id_2d ({1,2}); + OctaSym s ({-2,1}); + + auto peibos_output_pend = PEIBOS(vectorField_pend, 20., 0.2, psi0_pend, {id_2d,s,s*s,s.invert()}, 0.02, {-M_PI/2.,0.}); + + auto m_v_par_2d_pend = reach_set(peibos_output_pend); + // [peibos-capd-3-end] + } + + { + // [peibos-capd-4-beg] + capd::IMap vectorField_lorenz("par:sigma,rho,beta;var:x1,x2,x3;fun:10*(x2-x1),28*x1-x2-x1*x3,-2.6*x3+x1*x2;"); + vectorField_lorenz.setParameter("sigma", 10.); + vectorField_lorenz.setParameter("rho", 28.); + vectorField_lorenz.setParameter("beta", 8/3); + + VectorVar X_3d(2); + AnalyticFunction psi0_lorenz ({X_3d},{1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))}); + OctaSym id_3d ({1,2,3}); + OctaSym s1 ({-2,1,3}); + OctaSym s2 ({3,2,-1}); + + auto peibos_output_lorenz = PEIBOS(vectorField_lorenz, 0.1, 0.05, psi0_lorenz, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 0.1); + + auto m_v_par_lorenz = reach_set(peibos_output_lorenz); + // [peibos-capd-4-end] + } } \ No newline at end of file diff --git a/doc/manual/manual/functions/peibos/peibos.rst b/doc/manual/manual/functions/peibos/peibos.rst index 6afde2fb2..6a36de1ec 100644 --- a/doc/manual/manual/functions/peibos/peibos.rst +++ b/doc/manual/manual/functions/peibos/peibos.rst @@ -8,7 +8,7 @@ The PEIBOS tool provides a way to compute the Parallelepipedic Enclosure of the Let us consider an initial set :math:`\mathbb{X}_0 \subset \mathbb{R}^n` with its boundary :math:`\partial \mathbb{X}_0`. Considering a function :math:`\mathbf{f}:\mathbb{R}^n \to \mathbb{R}^p`, :math:`n \leq p`, the PEIBOS tool allows to compute the set :math:`\mathbf{Y}=\left\{ \mathbf{f}(\mathbf{x}) \mid \mathbf{x} \in \partial \mathbb{X}_0 \right\}`. - +.. _subsec-functions-peibos-gnomonic-atals: Gnomonic atlas -------------- diff --git a/doc/manual/manual/visualization/colors.rst b/doc/manual/manual/visualization/colors.rst index 9cb20760f..fced45d0c 100644 --- a/doc/manual/manual/visualization/colors.rst +++ b/doc/manual/manual/visualization/colors.rst @@ -308,7 +308,7 @@ A string starting with ``"w:"`` followed by a float can be passed to define the Z-value ~~~~~~~~ -**Warning: the Z-value has been added to VIBes since PR #150, a release is being prepared** +**Warning: the Z-value has been added to VIBes since version v0.3.4, please update to this version (or higher) for the Z-value to work** A string starting with ``"z:"`` followed by a float can be passed to define the Z-value (default is 0) diff --git a/examples/11_peibos/main.cpp b/examples/11_peibos/main.cpp index d87f75b11..03059c81e 100644 --- a/examples/11_peibos/main.cpp +++ b/examples/11_peibos/main.cpp @@ -5,6 +5,8 @@ using namespace codac2; int main() { + set_threads_used(max_threads()); + // 2D example of the PEIBOS algorithm VectorVar y_2d(2); double a = 1.4; double b = 0.3; @@ -49,7 +51,7 @@ int main() figure_3d_proj.set_window_properties({25,600},{500,500}); figure_3d_proj.set_axes({0,{-1.5,2.5}}, {1,{-2,2}}); - auto v_par_3d = PEIBOS(f_3d, psi0_3d, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 0.2, true); + auto v_par_3d = PEIBOS(f_3d, psi0_3d, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 0.05, true); for (const auto& p : v_par_3d) { diff --git a/examples/11_peibos/main.m b/examples/11_peibos/main.m index 6db6169f8..349361fbb 100644 --- a/examples/11_peibos/main.m +++ b/examples/11_peibos/main.m @@ -49,7 +49,7 @@ s1 = OctaSym(int64([-2, 1, 3])); s2 = OctaSym(int64([3, 2, -1])); -v_par_3d = PEIBOS(f_3d,psi0_3d,{id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()},0.2,true); +v_par_3d = PEIBOS(f_3d,psi0_3d,{id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()},0.05,true); figure_3d = Figure3D("Conform"); figure_3d.draw_axes(); diff --git a/examples/11_peibos/main.py b/examples/11_peibos/main.py index 95927228a..e6e3d29a0 100644 --- a/examples/11_peibos/main.py +++ b/examples/11_peibos/main.py @@ -2,6 +2,7 @@ import numpy as np if __name__=="__main__": + set_threads_used(max_threads()) # 2D example of the PEIBOS algorithm @@ -46,7 +47,7 @@ figure_3d_proj.set_window_properties([25,600],[500,500]) figure_3d_proj.set_axes(axis(0,[-1.5,2.5]), axis(1,[-2,2])) - v_par_3d = PEIBOS(f_3d,psi0_3d,[id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()],0.2,True) + v_par_3d = PEIBOS(f_3d,psi0_3d,[id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()],0.05,True) for p in v_par_3d: figure_3d.draw_parallelepiped(p,Color.green(0.5)) diff --git a/examples/12_peibos_capd/main.cpp b/examples/12_peibos_capd/main.cpp index 6b0b1e9a1..6ebd82118 100644 --- a/examples/12_peibos_capd/main.cpp +++ b/examples/12_peibos_capd/main.cpp @@ -6,34 +6,36 @@ using namespace codac2; int main() { + set_threads_used(max_threads()); + + ColorMap cmap = ColorMap::rainbow(); + VectorVar X_2d(1); + VectorVar X_3d(2); + // Julien's example capd::IMap vectorField_2d("var:x1,x2; fun:1, sin(x1);"); double tf_2d = 8.0; + double dt_2d = 2.0; - VectorVar X_2d(1); AnalyticFunction psi0_2d ({X_2d},{0.5*X_2d[0],0.5}); OctaSym id_2d ({1,2}); OctaSym s ({-2,1}); - auto peibos_output = PEIBOS(vectorField_2d, tf_2d, psi0_2d, {id_2d,s,s*s,s.invert()}, 0.1, true); + auto peibos_output = PEIBOS(vectorField_2d, tf_2d, dt_2d, psi0_2d, {id_2d,s,s*s,s.invert()}, 0.1, true); Figure2D output ("julien_3_4",GraphicOutput::VIBES); output.set_axes(axis(0,{-1,10}),axis(1,{-1,3})); - output.set_window_properties({50,100},{800,800}); + output.set_window_properties({25,100},{800,800}); - ColorMap cmap = ColorMap::rainbow(); - vector time_samples_2d {1e-20,2.0,4.0,6.0,8.0}; + auto m_v_par_2d = reach_set(peibos_output); - for (double t : time_samples_2d) - { - auto v_par_2d = reach_set(peibos_output,t); + for (auto& [t, v_par_2d] : m_v_par_2d) for (const auto& par: v_par_2d) output.draw_parallelepiped(par, cmap.color(t/tf_2d)); - } // Pendulum @@ -43,19 +45,47 @@ int main() vectorField_pend.setParameter("g",capd::interval(10.)); double tf_pend = 20.0; + double dt_pend = 0.2; - AnalyticFunction psi0_pend ({X_2d},{1e-2*X_2d[0],1e-2}); + AnalyticFunction psi0_pend ({X_2d},{0.1*X_2d[0],0.1}); - auto peibos_output_pend = PEIBOS(vectorField_pend, tf_pend, psi0_pend, {id_2d,s,s*s,s.invert()}, 0.1, {-M_PI/2.,0.}, true); + auto peibos_output_pend = PEIBOS(vectorField_pend, tf_pend, dt_pend, psi0_pend, {id_2d,s,s*s,s.invert()}, 0.02, {-M_PI/2.,0.}, true); Figure2D output_pend ("Pendulum",GraphicOutput::VIBES | GraphicOutput::IPE); - output_pend.set_axes(axis(0,{-2,2}),axis(1,{-4,4})); - output_pend.set_window_properties({50,100},{800,800}); + output_pend.set_axes(axis(0,{-2,2}),axis(1,{-3,3})); + output_pend.set_window_properties({800,100},{800,800}); + + auto m_v_par_2d_pend = reach_set(peibos_output_pend); - for (double t = 1e-20; t <= tf_pend; t += 0.05) - { - auto v_par_2d = reach_set(peibos_output_pend,t); + for (auto& [t, v_par_2d] : m_v_par_2d_pend) for (const auto& par: v_par_2d) output_pend.draw_parallelepiped(par, cmap.color(t/tf_pend)); - } + + // Lorenz + + Figure3D figure_lorenz ("Lorenz"); + ColorMap cmap_lorenz = ColorMap::rainbow(0.5); + + capd::IMap vectorField_lorenz("par:sigma,rho,beta;var:x1,x2,x3;fun:10*(x2-x1),28*x1-x2-x1*x3,-2.6*x3+x1*x2;"); + vectorField_lorenz.setParameter("sigma", 10.); + vectorField_lorenz.setParameter("rho", 28.); + vectorField_lorenz.setParameter("beta", 8/3); + + double tf_lorenz = 0.1; + double dt_lorenz = 0.05; + + AnalyticFunction psi0_lorenz ({X_3d},{1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))}); + + OctaSym id_3d ({1,2,3}); + OctaSym s1 ({-2,1,3}); + OctaSym s2 ({3,2,-1}); + + auto peibos_output_lorenz = PEIBOS(vectorField_lorenz, tf_lorenz, dt_lorenz, psi0_lorenz, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 0.1, true); + + auto m_v_par_lorenz = reach_set(peibos_output_lorenz); + + for (auto& [t, v_par_lorenz] : m_v_par_lorenz) + for (const auto& par: v_par_lorenz) + figure_lorenz.draw_parallelepiped(par, {cmap_lorenz.color(t/tf_lorenz),"t:"+to_string(t)}); + } \ No newline at end of file diff --git a/python/src/core/CMakeLists.txt b/python/src/core/CMakeLists.txt index f29275756..1aaef6974 100644 --- a/python/src/core/CMakeLists.txt +++ b/python/src/core/CMakeLists.txt @@ -116,6 +116,7 @@ tools/codac2_py_Approx.cpp tools/codac2_py_RobotSimulator.cpp tools/codac2_py_serialization.cpp + tools/codac2_py_threading.cpp tools/codac2_py_transformations.cpp tools/codac2_py_trunc.cpp diff --git a/python/src/core/codac2_py_core.cpp b/python/src/core/codac2_py_core.cpp index 35884a859..b11422f4c 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -152,6 +152,7 @@ void export_SepWrapper(py::module& m, py::class_& sep); void export_Approx(py::module& m); void export_RobotSimulator(py::module& m); void export_serialization(py::module& m); +void export_threading(py::module& m); void export_transformations(py::module& m); void export_trunc(py::module& m); @@ -317,6 +318,7 @@ PYBIND11_MODULE(_core, m) // tools export_Approx(m); export_serialization(m); + export_threading(m); export_transformations(m); export_trunc(m); export_RobotSimulator(m); diff --git a/python/src/core/tools/codac2_py_threading.cpp b/python/src/core/tools/codac2_py_threading.cpp new file mode 100644 index 000000000..bfd2f29f1 --- /dev/null +++ b/python/src/core/tools/codac2_py_threading.cpp @@ -0,0 +1,32 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include "codac2_py_threading_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_threading(py::module& m) +{ + m.def("max_threads", &codac2::max_threads, + INT_MAX_THREADS); + + m.def("set_threads_used", &codac2::set_threads_used, + VOID_SET_THREADS_USED_INT, + "n"_a); + + m.def("get_threads_used", &codac2::get_threads_used, + INT_GET_THREADS_USED); +} \ No newline at end of file diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index ffec45e23..39c9c5c18 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -255,6 +255,8 @@ ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_RobotSimulator.h ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_serialization.h ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_template_tools.h + ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_threading.h + ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_threading.cpp ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_transformations.cpp ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_transformations.h ${CMAKE_CURRENT_SOURCE_DIR}/tools/codac2_trunc.cpp @@ -309,7 +311,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}/tools ${CMAKE_CURRENT_SOURCE_DIR}/trajectory ) - target_link_libraries(${PROJECT_NAME}-core PUBLIC Ibex::ibex Eigen3::Eigen) + target_link_libraries(${PROJECT_NAME}-core PUBLIC Ibex::ibex Eigen3::Eigen Threads::Threads) ################################################################################ diff --git a/src/core/functions/analytic/codac2_AnalyticExpr.h b/src/core/functions/analytic/codac2_AnalyticExpr.h index cc27110b2..3adad48de 100644 --- a/src/core/functions/analytic/codac2_AnalyticExpr.h +++ b/src/core/functions/analytic/codac2_AnalyticExpr.h @@ -31,7 +31,7 @@ namespace codac2 virtual std::pair output_shape() const = 0; const T& init_value(ValuesMap& v, const T& x) const - { + { auto& p = v[unique_id()]; if(!p) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 8cbfc70d8..f0bc32323 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -12,6 +12,10 @@ #include "codac2_peibos.h" #include "codac2_peibos_tools.h" #include "codac2_OctaSym_operator.h" +#include "codac2_threading.h" + +#include +#include using namespace codac2; @@ -50,28 +54,66 @@ namespace codac2 assert_release (m < psi_0.output_size()); assert_release (Sigma.size() > 0 && (int) Sigma[0].size() == psi_0.output_size() && "no generator given or wrong dimension of generator (must match output size of psi_0)"); - clock_t t_start = clock(); - - std::vector output; + auto start_time = std::chrono::high_resolution_clock::now(); std::vector boxes; double true_eps = split(IntervalVector::constant(m,{-1,1}), epsilon, boxes); + int nthreads = get_threads_used(); + std::vector> thread_outputs(nthreads); + + struct WorkItem { const OctaSym* sigma; const IntervalVector* box; }; + std::vector work; + work.reserve(Sigma.size() * boxes.size()); for (const auto& sigma : Sigma) + for (const auto& box : boxes) + work.push_back({&sigma, &box}); + + std::vector local_output; + VectorVar x(m); + + auto worker = [&](int start, int end, int tid) { - VectorVar x(m); - AnalyticFunction g_i ({x}, f(sigma(psi_0(x))+offset)); + auto& local_output = thread_outputs[tid]; + for (int i = start; i < end; ++i) + { + const auto& sigma = *work[i].sigma; + const auto& X = *work[i].box; + + AnalyticFunction g_i ({x}, f(sigma(psi_0(x))+offset)); + local_output.push_back(g_i.parallelepiped_eval(X)); + } + }; + + std::vector threads; + int chunk_size = (int(work.size()) + nthreads - 1) / nthreads; - for (const auto& X : boxes) - output.push_back(g_i.parallelepiped_eval(X)); + for (int t = 0; t < nthreads; ++t) + { + int start = t * chunk_size; + int end = std::min(start + chunk_size, (int)work.size()); + if (start >= end) break; + threads.emplace_back(worker, start, end, t); } + for (auto& th : threads) th.join(); + + std::vector output; + output.reserve(Sigma.size() * boxes.size()); + + for (auto& vec : thread_outputs) + for (auto& el : vec) + output.emplace_back(std::move(el)); + + if (verbose) { printf("\nPEIBOS statistics:\n"); printf("------------------\n"); printf("Real epsilon: %.4f\n", true_eps); - printf("Computation time: %.4fs\n\n", (double)(clock()-t_start)/CLOCKS_PER_SEC); + printf("Number of thread used: %d\n", nthreads); + std::chrono::duration elapsed = std::chrono::high_resolution_clock::now() - start_time; + printf("Computation time: %.4fs\n\n", elapsed.count()); } return output; diff --git a/src/core/tools/codac2_threading.cpp b/src/core/tools/codac2_threading.cpp new file mode 100644 index 000000000..a95f438dd --- /dev/null +++ b/src/core/tools/codac2_threading.cpp @@ -0,0 +1,40 @@ +/** + * codac2_threading.cpp + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_threading.h" + +using namespace std; +using namespace codac2; + +namespace +{ + int threads_used = 1; +} + +namespace codac2 +{ + // FOR SIMON : when is this resolved ? + // i.e. if I compile the library on a machine (say a github one), and used it on my personnal machine + // will I get the hardware_concurrency of my machine or the github one ? + // And what if I just write const int max_threads = std::thread::hardware_concurrency(); ? + int max_threads() + { + return std::thread::hardware_concurrency(); + } + + void set_threads_used(int n) + { + threads_used = (n > 1 ? std::min(n,max_threads()) : 1); + } + + int get_threads_used() + { + return threads_used; + } +} \ No newline at end of file diff --git a/src/core/tools/codac2_threading.h b/src/core/tools/codac2_threading.h new file mode 100644 index 000000000..dfbfbad7e --- /dev/null +++ b/src/core/tools/codac2_threading.h @@ -0,0 +1,25 @@ +/** + * \file codac2_threading.h + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include + +namespace codac2_threading +{ + inline int threads_used = 1; +} + +namespace codac2 +{ + int max_threads(); + void set_threads_used(int n); + int get_threads_used(); + +} \ No newline at end of file diff --git a/src/extensions/capd/codac2_peibos_capd.cpp b/src/extensions/capd/codac2_peibos_capd.cpp index 8ed99bf03..0105f66ad 100644 --- a/src/extensions/capd/codac2_peibos_capd.cpp +++ b/src/extensions/capd/codac2_peibos_capd.cpp @@ -10,103 +10,140 @@ #include #include "codac2_AnalyticFunction.h" #include "codac2_peibos_capd.h" +#include "codac2_threading.h" + +#include +#include using namespace std; namespace codac2 { - vector>> PEIBOS(const capd::IMap& i_map, double tf, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, bool verbose) + std::map> PEIBOS(const capd::IMap& i_map, double tf, double dt, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, bool verbose) { - return PEIBOS(i_map, tf, psi_0, Sigma, epsilon, Vector::zero(psi_0.output_size()), verbose); + return PEIBOS(i_map, tf, dt, psi_0, Sigma, epsilon, Vector::zero(psi_0.output_size()), verbose); } - vector>> PEIBOS(const capd::IMap& i_map, double tf, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, const Vector& offset, bool verbose) + std::map> PEIBOS(const capd::IMap& i_map, double tf, double dt, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, const Vector& offset, bool verbose) { - int m = psi_0.input_size(); - int n = psi_0.output_size(); + std::vector time_points; + for (double t = 0.; t <= tf; t += dt) + time_points.push_back(t); + + if (time_points.back() < tf) + time_points.push_back(tf); - assert_release(offset.size() == n); - assert_release(m < n); - assert_release(Sigma.size() > 0 && (int) Sigma[0].size() == n); + int m = psi_0.input_size(); - clock_t t_start = clock(); + assert_release(offset.size() == psi_0.output_size()); + assert_release(m < psi_0.output_size()); + assert_release(Sigma.size() > 0 && (int) Sigma[0].size() == psi_0.output_size()); - vector>> output; - - // CAPD solver setup - capd::IMap g (i_map); - capd::IOdeSolver solver(g, 30); - - capd::ITimeMap timeMap(solver); - capd::ITimeMap timeMap_punct(solver); + auto start_time = std::chrono::high_resolution_clock::now(); capd::interval initialTime(0.); capd::interval finalTime(tf); vector boxes; double true_eps = split(Interval(-1.,1.)*IntervalVector::Ones(m), epsilon, boxes); - + + int nthreads = get_threads_used(); + std::vector>> thread_outputs(nthreads); + + struct WorkItem { const OctaSym* sigma; const IntervalVector* box; }; + std::vector work; + work.reserve(Sigma.size() * boxes.size()); for (const auto& sigma : Sigma) + for (const auto& box : boxes) + work.push_back({&sigma, &box}); + + auto worker = [&](int start, int end, int tid) { - for (const auto& X : boxes) + auto& local_output = thread_outputs[tid]; + + capd::IMap g(i_map); + capd::IOdeSolver solver(g, 30); + + capd::ITimeMap timeMap(solver); + capd::ITimeMap timeMap_punct(solver); + + for (int i = start; i < end; ++i) { + const auto& sigma = *work[i].sigma; + const auto& X = *work[i].box; - PEIBOS_CAPD_Key key {X, psi_0, sigma, offset}; + PEIBOS_CAPD_Key key{X, psi_0, sigma, offset}; - // To get the flow function and its Jacobian (monodromy matrix) for [x] IntervalVector Y = sigma(psi_0.eval(X)) + offset; - - capd::IMatrix monodromyMatrix(n,n); - capd::ITimeMap::SolutionCurve solution(initialTime); + capd::ITimeMap::SolutionCurve solution(initialTime); capd::IVector c = to_capd(Y); - capd::C1Rect2Set s(c); timeMap(finalTime, s, solution); - // To get the flow function and its Jacobian (monodromy matrix) for x_hat auto xc = X.mid(); auto yc = (sigma(psi_0.eval(xc)) + offset).mid(); - - capd::IMatrix monodromyMatrix_punct(n,n); capd::ITimeMap::SolutionCurve solution_punct(initialTime); capd::IVector c_punct = to_capd(IntervalVector(yc)); - capd::C1Rect2Set s_punct(c_punct); - timeMap_punct(finalTime, s_punct, solution_punct); + timeMap_punct(finalTime, s_punct, solution_punct); + + for (auto t : time_points) + { + IntervalVector z = to_codac(solution_punct(t)); + IntervalMatrix Jf = to_codac(solution.derivative(t)); + local_output[t].emplace_back(key, z, Jf); + } + } + }; - output.push_back(make_pair(key, make_pair(solution, solution_punct))); + std::vector threads; + int chunk_size = (int(work.size()) + nthreads - 1) / nthreads; - } + for (int t = 0; t < nthreads; ++t) + { + int start = t * chunk_size; + int end = std::min(start + chunk_size, (int)work.size()); + if (start >= end) break; + threads.emplace_back(worker, start, end, t); } + + for (auto& th : threads) th.join(); + + std::map> output; + + for (auto& vec : thread_outputs) + for (auto t : time_points) + for (auto& item : vec[t]) + output[t].push_back(item); + if (verbose) { printf("\nPEIBOS statistics:\n"); printf("------------------\n"); printf("Real epsilon: %.4f\n", true_eps); - printf("Computation time: %.4fs\n\n", (double)(clock()-t_start)/CLOCKS_PER_SEC); + printf("Number of thread used: %d\n", nthreads); + std::chrono::duration elapsed = std::chrono::high_resolution_clock::now() - start_time; + printf("Computation time: %.4fs\n\n", elapsed.count()); } return output; } - vector reach_set(const vector>>& peibos_output, double t) + std::map> reach_set(const std::map>& peibos_output) { - vector output; + std::map> output; - for (const auto& [key,flow_pair] : peibos_output) + for (const auto& [time,vec] : peibos_output) { - const auto& [flow, flow_punct] = flow_pair; - - IntervalVector z = to_codac(flow_punct(t)); - auto Jf_tild = (to_codac(flow_punct.derivative(t))).mid(); - auto Jf = to_codac(flow.derivative(t)); - - auto p = parallelepiped_inclusion(z, Jf, Jf_tild, key.psi_0, key.sigma, key.box); + for (const auto& [key, z, Jf] : vec) + { + auto p = parallelepiped_inclusion(z, Jf, Jf.mid(), key.psi_0, key.sigma, key.box); - output.push_back(p); + output[time].push_back(p); + } } - + return output; } diff --git a/src/extensions/capd/codac2_peibos_capd.h b/src/extensions/capd/codac2_peibos_capd.h index 955c7b3cd..9aab7cfd0 100644 --- a/src/extensions/capd/codac2_peibos_capd.h +++ b/src/extensions/capd/codac2_peibos_capd.h @@ -10,6 +10,7 @@ #pragma once #include +#include #include "codac2_peibos.h" #include "codac2_capd.h" #include "codac2_OctaSym.h" @@ -17,16 +18,71 @@ namespace codac2 { + /** + * \struct PEIBOS_CAPD_Key + * \brief Key structure for PEIBOS CAPD results. + * + * A PEIBOS_CAPD_Key represents \f$\mathbf{y}= \sigma(\psi_0(\mathbf{\text{box}})) + \text{offset}\f$. + */ struct PEIBOS_CAPD_Key { + /// @brief The box \f$[\mathbf{x}]\f$ IntervalVector box; + /// @brief The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas AnalyticFunction psi_0; + /// @brief The symmetry operator \f$\sigma\f$ to construct the atlas OctaSym sigma; + /// @brief The offset to add to \f$\sigma(\psi_0([x]))\f$ Vector offset; }; - std::vector>> PEIBOS(const capd::IMap& i_map, double tf, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, bool verbose = false); - std::vector>> PEIBOS(const capd::IMap& i_map, double tf, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, const Vector& offset, bool verbose = false); + using T = std::tuple; - std::vector reach_set(const std::vector>>& peibos_output, double t); + /** + * \brief PEIBOS algorithm using CAPD for guaranteed ODE propagation. + * + * \param i_map The CAPD interval map representing the ODE. + * \param tf Final time for the propagation. + * \param dt Time step for the output map. + * \param psi_0 The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param Sigma The set of symmetry operators \f$\sigma\f$ to construct the atlas + * \param epsilon The maximum diameter of the boxes to split \f$[-1,1]^m\f$ before computing the parallelepiped inclusions (each box is called "box" below) + * \param verbose If true, print the time taken to compute the parallelepiped inclusions with other statistics + * + * \return A timed map of PEIBOS CAPD results. At each time \f$t\f$, the value is a vector of tuples. Each tuple contains: + * \li A PEIBOS_CAPD_Key representing \f$\mathbf{x}(0)= \sigma(\psi_0(\mathbf{\text{box}}))\f$ + * \li The interval vector \f$\mathbf{z}\f$ containing the image \f$\bar{\mathbf{x}}(t))\f$ + * \li The interval Jacobian matrix \f$\mathbf{J_f}\f$ containing \f$D\mathbf{\left[x\right]}(t)\f$ + */ + std::map> PEIBOS(const capd::IMap& i_map, double tf, double dt, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, bool verbose = false); + + /** + * \brief PEIBOS algorithm using CAPD for guaranteed ODE propagation. + * + * \param i_map The CAPD interval map representing the ODE. + * \param tf Final time for the propagation. + * \param dt Time step for the output map. + * \param psi_0 The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param Sigma The set of symmetry operators \f$\sigma\f$ to construct the atlas + * \param epsilon The maximum diameter of the boxes to split \f$[-1,1]^m\f$ before computing the parallelepiped inclusions (each box is called "box" below) + * \param offset The offset to add to \f$\sigma(\psi_0([-1,1]^m))\f$ (used to translate the initial manifold) + * \param verbose If true, print the time taken to compute the parallelepiped inclusions with other statistics + * + * \return A timed map of PEIBOS CAPD results. At each time \f$t\f$, the value is a vector of tuples. Each tuple contains: + * \li A PEIBOS_CAPD_Key representing \f$\mathbf{x}(0)= \sigma(\psi_0(\mathbf{\text{box}})) + \text{offset}\f$ + * \li The interval vector \f$\mathbf{z}\f$ containing the image \f$\bar{\mathbf{x}}(t))\f$ + * \li The interval Jacobian matrix \f$\mathbf{J_f}\f$ containing \f$D\mathbf{\left[x\right]}(t)\f$ + */ + std::map> PEIBOS(const capd::IMap& i_map, double tf, double dt, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, const Vector& offset, bool verbose = false); + + + /** + * \brief Compute the reach set parallelepipeds from the PEIBOS CAPD output. + * + * \param peibos_output The output of the PEIBOS CAPD algorithm. + * + * \return A timed map of reach set parallelepipeds. At each time \f$t\f$, the value is a vector of Parallelepipeds enclosing the reach set at time \f$t\f$. + * The function \ref parallelepiped_inclusion is used to compute each Parallelepiped from the PEIBOS CAPD output. + */ + std::map> reach_set(const std::map>& peibos_output); } \ No newline at end of file diff --git a/src/graphics/figures/codac2_Figure2D.cpp b/src/graphics/figures/codac2_Figure2D.cpp index 794df3b56..c476ab806 100644 --- a/src/graphics/figures/codac2_Figure2D.cpp +++ b/src/graphics/figures/codac2_Figure2D.cpp @@ -278,10 +278,13 @@ void Figure2D::draw_parallelepiped(const Parallelepiped& p, const StylePropertie auto a1 = p.A.col(0), a2 = p.A.col(1); - draw_polygon(vector({ - Vector(p.z+a1+a2), Vector(p.z-a1+a2), - Vector(p.z-a1-a2), Vector(p.z+a1-a2) - }), style); + if (a1.isZero() || a2.isZero()) + draw_polyline(vector({p.z-a1-a2,p.z+a1+a2}), style); + else + draw_polygon(vector({ + Vector(p.z+a1+a2), Vector(p.z-a1+a2), + Vector(p.z-a1-a2), Vector(p.z+a1-a2) + }), style); } void Figure2D::draw_pie(const Vector& c, const Interval& r, const Interval& theta, const StyleProperties& style) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5185c4be8..a9165f1cd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -119,6 +119,7 @@ set(CODAC_LIBRARIES ${PROJECT_NAME}-core ${PROJECT_NAME}-graphics) if (WITH_CAPD) list(APPEND SRC_TESTS extensions/capd/codac2_tests_capd + extensions/capd/codac2_tests_peibos_capd ../doc/manual/manual/extensions/capd/src ) set(CODAC_LIBRARIES ${CODAC_LIBRARIES} ${PROJECT_NAME}-capd capd::capd) diff --git a/tests/extensions/capd/codac2_tests_peibos_capd.cpp b/tests/extensions/capd/codac2_tests_peibos_capd.cpp new file mode 100644 index 000000000..ffa1eaaa8 --- /dev/null +++ b/tests/extensions/capd/codac2_tests_peibos_capd.cpp @@ -0,0 +1,26 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("PEIBOS_CAPD") +{ + { + capd::IMap vectorField_pend("var:t,w;fun:w,-sin(t) - 0.5*w;"); + VectorVar X(1); + AnalyticFunction psi0_pend ({X},{0.1*X[0],0.1}); + auto peibos_output_pend = PEIBOS(vectorField_pend, 1.0, 0.5, psi0_pend, {OctaSym({1,2})}, 0.2, {-M_PI/2.,0.}); + auto m_v_par_2d_pend = reach_set(peibos_output_pend); + } +} \ No newline at end of file diff --git a/tests/extensions/capd/codac2_tests_peibos_capd.py b/tests/extensions/capd/codac2_tests_peibos_capd.py new file mode 100644 index 000000000..63a8b0cdf --- /dev/null +++ b/tests/extensions/capd/codac2_tests_peibos_capd.py @@ -0,0 +1,19 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2026 +# \author Maël Godard +# \copyright Copyright 2024 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +from codac import * + +class TestPEIBOSCAPD(unittest.TestCase): + + def test_peibos_capd(self): + self.assertTrue(True) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file