diff --git a/.github/workflows/dockercentos.yml b/.github/workflows/dockercentos.yml index e81b929ac..71b2da456 100644 --- a/.github/workflows/dockercentos.yml +++ b/.github/workflows/dockercentos.yml @@ -22,7 +22,7 @@ jobs: - { img: 'lebarsfa/manylinux_2_28_aarch64-for-codac', os: ubuntu-24.04-arm, shell: bash, arch: aarch64, bitness: 64, runtime: manylinux_2_28, docker_flags: '--platform linux/arm64', desc: 'CentOS manylinux_2_28 aarch64' } name: ${{ matrix.cfg.desc }} steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 with: submodules: true fetch-depth: 0 diff --git a/.github/workflows/dockermatrix.yml b/.github/workflows/dockermatrix.yml index bf9ccccad..7660553b3 100644 --- a/.github/workflows/dockermatrix.yml +++ b/.github/workflows/dockermatrix.yml @@ -29,7 +29,7 @@ jobs: #- { img: 'lebarsfa/pi-64:focal-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: arm64 , bitness: 64, runtime: focal , cmake_flags: '-fPIC', deb: true, desc: 'Ubuntu 20.04 arm64' } - { img: 'lebarsfa/amd64:trixie-for-codac' , os: ubuntu-latest , shell: bash, arch: amd64 , bitness: 64, runtime: trixie , cmake_flags: '-fPIC', deb: true, desc: 'Debian Trixie amd64' } - { img: 'lebarsfa/pi-64:trixie-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: arm64 , bitness: 64, runtime: trixie , cmake_flags: '-fPIC', deb: true, desc: 'Debian Trixie arm64' } - #- { img: 'lebarsfa/pi:trixie-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: armhf , bitness: 32, runtime: trixie , cmake_flags: '-fPIC', deb: true, desc: 'Raspbian Trixie armv6hf' } # No base image available... + - { img: 'lebarsfa/pi:trixie-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: armhf , bitness: 32, runtime: trixie , cmake_flags: '-fPIC', deb: true, desc: 'Raspbian Trixie armv6hf' } - { img: 'lebarsfa/amd64:bookworm-for-codac' , os: ubuntu-latest , shell: bash, arch: amd64 , bitness: 64, runtime: bookworm , cmake_flags: '-fPIC', deb: true, desc: 'Debian Bookworm amd64' } - { img: 'lebarsfa/pi-64:bookworm-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: arm64 , bitness: 64, runtime: bookworm , cmake_flags: '-fPIC', deb: true, desc: 'Debian Bookworm arm64' } - { img: 'lebarsfa/pi:bookworm-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: armhf , bitness: 32, runtime: bookworm , cmake_flags: '-fPIC', deb: true, desc: 'Raspbian Bookworm armv6hf' } @@ -42,7 +42,7 @@ jobs: #- { img: 'lebarsfa/pi:buster-for-codac' , os: ubuntu-24.04-arm, shell: bash, arch: armhf , bitness: 32, runtime: buster , cmake_flags: '-fPIC', deb: true, desc: 'Raspbian Buster armv6hf' } name: ${{ matrix.cfg.desc }} steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 with: submodules: true fetch-depth: 0 diff --git a/.github/workflows/macosmatrix.yml b/.github/workflows/macosmatrix.yml index da66789bb..683f9c234 100644 --- a/.github/workflows/macosmatrix.yml +++ b/.github/workflows/macosmatrix.yml @@ -34,7 +34,7 @@ jobs: - { os: macos-15-intel , shell: bash, arch: x86_64 , runtime: sequoia , cmake_flags: '-fPIC', trgt: '10.16', cpcfg: '-macosx_10_16_x86_64', py_v_maj: 3, py_v_min: 8 , desc: 'macOS Sequoia Python 3.8 x86_64' } name: ${{ matrix.cfg.desc }} steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 with: submodules: true fetch-depth: 0 @@ -54,7 +54,7 @@ jobs: # if: runner.os=='macOS' - run: brew install catch2 # Issues with binary packages when cross-compiling... if: (runner.os=='macOS')&&(matrix.cfg.cross!=true) - - run: brew install graphviz ; brew install --formula doxygen ; python -m pip install --upgrade pip ; pip install --upgrade wheel setuptools sphinx breathe sphinx_rtd_theme sphinx-tabs sphinx-issues sphinx-reredirects furo sphinx-math-dollar sympy + - run: brew install graphviz ; brew install --formula doxygen ; python -m pip install --upgrade pip ; pip install --upgrade wheel setuptools sphinx breathe sphinx_rtd_theme sphinx-tabs sphinx-issues sphinx-reredirects furo sphinx-math-dollar sphinx_togglebutton sympy if: runner.os=='macOS' - run: | wget https://github.com/lebarsfa/ibex-lib/releases/download/ibex-2.8.9.20250626/ibex_${{ matrix.cfg.arch }}_${{ matrix.cfg.runtime }}.zip --no-check-certificate -nv diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d85e07323..7d5cfe493 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -20,7 +20,7 @@ jobs: - { os: ubuntu-22.04 , gcc_v: 11, py_v_maj: 3, py_v_min: 10, with_capd: ON , desc: 'Ubuntu 22.04 GCC 11 Python 3.10 CAPD tests' } name: ${{ matrix.cfg.desc }} steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 with: submodules: true fetch-depth: 0 diff --git a/.github/workflows/unixmatrix.yml b/.github/workflows/unixmatrix.yml index e6436dfdf..ddd637283 100644 --- a/.github/workflows/unixmatrix.yml +++ b/.github/workflows/unixmatrix.yml @@ -10,7 +10,7 @@ jobs: runs-on: ubuntu-latest if: (github.event_name!='pull_request')&&((github.ref_name=='codac1')||(github.ref_name=='codac2')||(github.ref_name=='codac2_codac4matlab')) steps: - - uses: softprops/action-gh-release@v2 + - uses: softprops/action-gh-release@v3 id: create_release env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} @@ -68,7 +68,7 @@ jobs: - { os: macos-14 , shell: bash, arch: x86_64 , bitness: 64, runtime: sonoma , cmake_params: '-D CMAKE_SYSTEM_NAME=Darwin -D CMAKE_OSX_ARCHITECTURES=x86_64 -D CMAKE_OSX_DEPLOYMENT_TARGET=14.0', cmake_flags: '-fPIC', cross: true, desc: 'macOS Sonoma x86_64 (cross)' } name: ${{ matrix.cfg.desc }} steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 with: submodules: true fetch-depth: 0 diff --git a/.github/workflows/vcmatrix.yml b/.github/workflows/vcmatrix.yml index 62b58767f..3039da799 100644 --- a/.github/workflows/vcmatrix.yml +++ b/.github/workflows/vcmatrix.yml @@ -35,7 +35,7 @@ jobs: - { os: windows-2022 , shell: cmd , arch: x64 , runtime: vc17 , cmake_params: '-G "Visual Studio 17" -T v143 -A x64' , cpcfg: '-win_amd64', py_v_maj: 3, py_v_min: 8 , desc: 'Windows Visual Studio 2022 x64 Python 3.8' } name: ${{ matrix.cfg.desc }} steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 with: submodules: true fetch-depth: 0 @@ -56,7 +56,7 @@ jobs: if: runner.os=='Windows' #- run: choco install -y -r --no-progress eigen --version=3.4.0.20240224 ${{ matrix.cfg.choco_flags }} # if: runner.os=='Windows' - - run: choco install -y -r --no-progress graphviz doxygen.install & python -m pip install --upgrade pip & pip install --upgrade wheel setuptools sphinx breathe sphinx-issues sphinx-tabs sphinx_rtd_theme sphinx-reredirects furo sphinx-math-dollar sympy + - run: choco install -y -r --no-progress graphviz doxygen.install & python -m pip install --upgrade pip & pip install --upgrade wheel setuptools sphinx breathe sphinx-issues sphinx-tabs sphinx_rtd_theme sphinx-reredirects furo sphinx-math-dollar sphinx_togglebutton sympy if: runner.os=='Windows' - run: | wget https://github.com/lebarsfa/ibex-lib/releases/download/ibex-2.8.9.20250626/ibex.2.8.9.20250626.nupkg --no-check-certificate -nv diff --git a/README.md b/README.md index f6e21a60f..b6c20bc7a 100644 --- a/README.md +++ b/README.md @@ -81,8 +81,8 @@ This list is in alphabetical order by surname. - [Peter Franek](http://www.cs.cas.cz/~franek) - [Maël Godard](https://godardma.github.io) - [Nuwan Herath M.](https://members.loria.fr/NHerathMudiyanselage) -- [Luc Jaulin](https://www.ensta-bretagne.fr/jaulin) -- [Fabrice Le Bars](https://www.ensta-bretagne.fr/lebars) +- [Luc Jaulin](https://webperso.ensta.fr/jaulin) +- [Fabrice Le Bars](https://webperso.ensta.fr/lebars) - [Morgan Louédec](https://morgan-louedec.fr) - [Damien Massé](https://lab-sticc.univ-brest.fr/~dmasse) - [Bertrand Neveu](http://imagine.enpc.fr/~neveub) @@ -91,7 +91,7 @@ This list is in alphabetical order by surname. - [Simon Rohou](http://simon-rohou.fr/research) - [Joris Tillet](https://perso.ensta-paris.fr/~tillet/) - [Gilles Trombettoni](https://www.lirmm.fr/~trombetton) -- [Christophe Viel](https://www.ensta-bretagne.fr/viel/) +- [Christophe Viel](https://webperso.ensta.fr/viel/) - [Raphael Voges](https://raphael-voges.de) We appreciate all contributions, whether they are code, documentation, bug reports, or suggestions. If you believe you should be listed here and are not, please contact us to update the list. diff --git a/doc/api/Doxyfile.in b/doc/api/Doxyfile.in index 26b09d2c1..c4c4ed42f 100644 --- a/doc/api/Doxyfile.in +++ b/doc/api/Doxyfile.in @@ -2502,7 +2502,7 @@ HIDE_UNDOC_RELATIONS = YES # set to NO # The default value is: NO. -HAVE_DOT = YES +HAVE_DOT = NO # The DOT_NUM_THREADS specifies the number of dot invocations doxygen is allowed # to run in parallel. When set to 0 doxygen will base this on the number of diff --git a/doc/manual/development/changelog.rst b/doc/manual/development/changelog.rst index 4539c6e7d..1002ecf90 100644 --- a/doc/manual/development/changelog.rst +++ b/doc/manual/development/changelog.rst @@ -3,12 +3,19 @@ Changelog ========= -Pull Request from godardma (17/04/26) -------------------------------------- +Upcoming version +**************** +Pull Request #379 from godardma (27/04) +--------------------------------------- -Commit 52b81c8 ([cmake] warning for Doxygen version) ----------------------------------------------------- +``Zonotope`` and ``Parallelepiped`` center is now ``c`` instead of ``z`` to avoid ambiguity. + +Version 2.0.2 +************* + +Commit ca1f6f4 +-------------- set_axes change for Figure2D ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -16,6 +23,12 @@ set_axes change for Figure2D Since this PR, a bounding box can be passed to the set_axes method to specify the ranges of the x and y axes. The old method with ``axis(id,Interval)`` still works +Version 2.0.0 +************* + +Commit 52b81c8 ([cmake] warning for Doxygen version) +---------------------------------------------------- + Python binding build requirement ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -24,6 +37,8 @@ If this version is not available, the API documentation will not be generated, which may cause issues when building the Python binding from source. This has no impact on the C++ build nor on downloading/installing the Python packages. +Pre-Release +*********** Pull Request from godardma (15/10) ---------------------------------- diff --git a/doc/manual/development/info_dev.rst b/doc/manual/development/info_dev.rst index 0b55724d2..e99bf5a24 100644 --- a/doc/manual/development/info_dev.rst +++ b/doc/manual/development/info_dev.rst @@ -8,10 +8,11 @@ Building the Manual with Sphinx To build this manual using Sphinx, follow these steps: -1. **Requirements**: Install `Sphinx `_ together with the extensions used by Codac, if they are not already installed. +1. **Requirements**: Install `Sphinx `_ together with the extensions used by Codac, if they are not already installed. You will need also the ``dot`` command from `Graphviz `_. .. code-block:: bash + sudo apt-get install sphinx-common graphviz pip install sphinx sphinx_rtd_theme furo sphinx-math-dollar sphinx_tabs breathe sphinx_togglebutton 2. **Building the manual**: After configuring the CMake project, the manual can be build with: @@ -41,8 +42,8 @@ If you simply want to use the latest Codac release in Python, you can download t sudo apt-get install -y g++ gcc cmake git flex bison - - a supported version of Python (>=3.6). - - a recent `Doxygen `_ version (for instance, release 1.13.0 or newest). On Linux systems, latest releases are not available as Debian packages, so we advice to install Doxygen from the sources: + - a supported version of Python (>=3.8). + - a recent `Doxygen `_ version (for instance, release 1.16.1 or newest). On Linux systems, latest releases are not available as Debian packages, so we advise to install Doxygen from the sources: .. code-block:: bash @@ -95,12 +96,12 @@ If you simply want to use the latest Codac release in Python, you can download t Note that you will then have to ``import codac2`` instead of ``import codac`` in your Python scripts. - In addition to the ``-fPIC`` options, you will have to configure ``WITH_PYTHON=ON`` and ``PYBIND11_FINDPYTHON=OFF``. Note that CMake will automatically get the `pybind11 `_ files required for the binding. + In addition to the ``-fPIC`` options, you will have to configure ``WITH_PYTHON=ON`` and ``PYBIND11_FINDPYTHON=OFF``. Note that CMake will automatically get the `pybind11 `_ files required for the binding. Also, you will have to configure ``BUILD_TESTS=ON`` if you want to run the unit tests. .. code-block:: bash mkdir build ; cd build - cmake -DCMAKE_CXX_FLAGS="-fPIC" -DCMAKE_C_FLAGS="-fPIC" -DWITH_PYTHON=ON -DPYBIND11_FINDPYTHON=OFF -DCMAKE_INSTALL_PREFIX=$HOME/codac/build_install -DCMAKE_PREFIX_PATH="$HOME/ibex-lib/build_install;$HOME/doxygen/build_install" -DCMAKE_BUILD_TYPE=Release .. + cmake -DCMAKE_CXX_FLAGS="-fPIC" -DCMAKE_C_FLAGS="-fPIC" -DWITH_PYTHON=ON -DPYBIND11_FINDPYTHON=OFF -DBUILD_TESTS=ON -DCMAKE_INSTALL_PREFIX=$HOME/codac/build_install -DCMAKE_PREFIX_PATH="$HOME/ibex-lib/build_install;$HOME/doxygen/build_install" -DCMAKE_BUILD_TYPE=Release .. make ; make install 4. **Configure your Python environment**: @@ -135,7 +136,7 @@ If you simply want to use the latest Codac release in Python, you can download t cd $HOME/codac/examples/03_sivia python main.py - Note that before executing the example, you will have to launch the VIBes viewer. + Note that before executing the example, you will have to launch the `VIBes viewer <../visualization/vibes.html>`_. You should obtain a graphical output corresponding to a set inversion. diff --git a/doc/manual/index.rst b/doc/manual/index.rst index 8c38197fd..f217d3970 100644 --- a/doc/manual/index.rst +++ b/doc/manual/index.rst @@ -255,7 +255,7 @@ User manual * CtcCross / CtcNoCross * Shape contractors * CtcCtcBoundary - * CtcWrapper + * :ref:`sec-ctc-shape-ctcwrapper` * CtcImage * CtcDiscreteSet * Temporal contractors @@ -317,8 +317,8 @@ User manual * :ref:`sec-zonotope` * Polyhedron -* :ref:`sec-actions` - * :ref:`sec-actions-octasym` +* :ref:`sec-tools` + * :ref:`sec-tools-octasym` * :ref:`sec-ellipsoids` * :ref:`sec-ellipsoids-intro` diff --git a/doc/manual/manual/contractors/index.rst b/doc/manual/manual/contractors/index.rst index 0267838da..9b2537821 100644 --- a/doc/manual/manual/contractors/index.rst +++ b/doc/manual/manual/contractors/index.rst @@ -10,6 +10,7 @@ Contractors, separators CtcDist CtcPolar CtcVisible + CtcWrapper .. What are contractors? .. The Ctc class @@ -109,7 +110,7 @@ Overview of contractors and separators * - ``CtcCtcBoundary`` - ``SepCtcBoundary`` - * - ``CtcWrapper`` + * - :ref:`CtcWrapper ` - ``SepWrapper`` * - ``CtcImage`` diff --git a/doc/manual/manual/contractors/shape/ctcwrapper.rst b/doc/manual/manual/contractors/shape/ctcwrapper.rst new file mode 100644 index 000000000..17c7a625a --- /dev/null +++ b/doc/manual/manual/contractors/shape/ctcwrapper.rst @@ -0,0 +1,40 @@ +.. _sec-ctc-shape-ctcwrapper: + +CtcWrapper +=============================== + + Main author: `Maël Godard `_ + +The CtcWrapper is a contractor to contract a box with respect to a set represented by a Codac object. +Currently supported objects are: + +- :ref:`IntervalVector ` +- :ref:`Parallelepiped ` +- PavingOut + +Methods +------- + +Below is an example of the use of the CtcWrapper with a Parallelepiped. + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [ctcparallelepiped-1-beg] + :end-before: [ctcparallelepiped-1-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [ctcparallelepiped-1-beg] + :end-before: [ctcparallelepiped-1-end] + :dedent: 4 + +.. admonition:: Technical documentation + + See the `C++ API documentation of this class <../../api/html/classcodac2_1_1_ctc_wrapper.html>`_. \ No newline at end of file diff --git a/doc/manual/manual/contractors/shape/index.rst b/doc/manual/manual/contractors/shape/index.rst new file mode 100644 index 000000000..793066cca --- /dev/null +++ b/doc/manual/manual/contractors/shape/index.rst @@ -0,0 +1,9 @@ +Geometric contractors +===================== + +.. toctree:: + + CtcCtcBoundary + ctcwrapper.rst + CtcImage + CtcDiscreteSet \ No newline at end of file diff --git a/doc/manual/manual/contractors/shape/src.cpp b/doc/manual/manual/contractors/shape/src.cpp new file mode 100644 index 000000000..ab9c1eee3 --- /dev/null +++ b/doc/manual/manual/contractors/shape/src.cpp @@ -0,0 +1,32 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include + +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("CtcWrapper - Parallelepiped - manual") +{ + { + // [ctcparallelepiped-1-beg] + IntervalVector x ({{0,5}, {0,5}}); + Parallelepiped p (Vector({1.5,2.8}), Matrix({{0.5,0.4},{0,0.2}})); + CtcWrapper c(p); + c.contract(x); + // x = [ [0.599999, 2.40001] ; [2.59999, 3] ] + // [ctcparallelepiped-1-end] + + CHECK(Approx(x) == IntervalVector({{0.6, 2.4}, {2.6, 3}})); + } +} diff --git a/doc/manual/manual/contractors/shape/src.py b/doc/manual/manual/contractors/shape/src.py new file mode 100644 index 000000000..08344e907 --- /dev/null +++ b/doc/manual/manual/contractors/shape/src.py @@ -0,0 +1,30 @@ +#!/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 sys, os +import unittest +import math +from codac import * + +class TestCtcShapeManual(unittest.TestCase): + + def tests_CtcWrapper_Parallelepiped_manual(test): + + # [ctcparallelepiped-1-beg] + x = IntervalVector([[0,5], [0,5]]) + p = Parallelepiped(Vector([1.5,2.8]), Matrix([[0.5,0.4],[0,0.2]])) + c = CtcWrapper_Parallelepiped(p) + x = c.contract(x) + # x = [ [0.599999, 2.40001] ; [2.59999, 3] ] + # [ctcparallelepiped-1-end] + + test.assertTrue(Approx(x) == IntervalVector([[0.6, 2.4], [2.6, 3]])) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/doc/manual/manual/extensions/capd/capd.rst b/doc/manual/manual/extensions/capd/capd.rst index 70db69213..40c23731e 100644 --- a/doc/manual/manual/extensions/capd/capd.rst +++ b/doc/manual/manual/extensions/capd/capd.rst @@ -12,7 +12,6 @@ To use CAPD with Codac, you first need to install the CAPD library. You can find 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 --------------------------------------- @@ -52,6 +51,18 @@ The header of the ``codac-capd`` extension is not included by default. You need #include #include +Furthermore, you need to link the extension to your project, for instance by updating your ``CMakeLists.txt`` with:: + + add_executable(${PROJECT_NAME} main.cpp) + target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS}) + target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS} ${CODAC_CAPD_INCLUDE_DIR}) + target_link_libraries(${PROJECT_NAME} PRIVATE + ${CODAC_LIBRARIES} + ${CODAC_CAPD_LIBRARY} # linking to the codac-capd extension + capd::capd # linking to CAPD + Ibex::ibex + ) + You can use the functions ``to_capd`` and ``to_codac`` to convert between CAPD and Codac objects as follows: .. tabs:: @@ -200,4 +211,4 @@ The result is the following figure, with in green the initial set (:math:`t=0s`) .. figure:: img/pendulum_result.png :width: 500px - Result of the CAPD integration of the pendulum, enclosed in a Codac tube. + Result of the CAPD integration of the pendulum, enclosed in a Codac tube. \ No newline at end of file diff --git a/doc/manual/manual/extensions/capd/img/pendulum_peibos.png b/doc/manual/manual/extensions/capd/img/pendulum_peibos.png index 927185ce1..0841cc200 100644 Binary files a/doc/manual/manual/extensions/capd/img/pendulum_peibos.png 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 17ccff750..c030baa19 100644 --- a/doc/manual/manual/extensions/capd/index.rst +++ b/doc/manual/manual/extensions/capd/index.rst @@ -3,212 +3,7 @@ CAPD (rigorous numerics in dynamical systems) ============================================= - 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 - -Furthermore, you need to link the extension to your project, for instance by updating your ``CMakeLists.txt`` with:: - - add_executable(${PROJECT_NAME} main.cpp) - target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS}) - target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS} ${CODAC_CAPD_INCLUDE_DIR}) - target_link_libraries(${PROJECT_NAME} PRIVATE - ${CODAC_LIBRARIES} - ${CODAC_CAPD_LIBRARY} # linking to the codac-capd extension - capd::capd # linking to CAPD - Ibex::ibex - ) - -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 index 4d5c92216..fb80ca017 100644 --- a/doc/manual/manual/extensions/capd/peibos_capd.rst +++ b/doc/manual/manual/extensions/capd/peibos_capd.rst @@ -3,6 +3,10 @@ PEIBOS-CAPD =========== + Main author: `Maël Godard `_ + +The content from this page supports multi-threading, see :ref:`here ` for details. + 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`. diff --git a/doc/manual/manual/functions/peibos/peibos.rst b/doc/manual/manual/functions/peibos/peibos.rst index 6a36de1ec..f2afcf83f 100644 --- a/doc/manual/manual/functions/peibos/peibos.rst +++ b/doc/manual/manual/functions/peibos/peibos.rst @@ -3,6 +3,10 @@ PEIBOS ====== + Main author: `Maël Godard `_ + +The content from this page supports multi-threading, see :ref:`here ` for details. + The PEIBOS tool provides a way to compute the Parallelepipedic Enclosure of the Image of the BOundary of a Set. Let us consider an initial set :math:`\mathbb{X}_0 \subset \mathbb{R}^n` with its boundary :math:`\partial \mathbb{X}_0`. diff --git a/doc/manual/manual/geometry/zonotope.png b/doc/manual/manual/geometry/zonotope.png index fb74b6718..cc0584583 100644 Binary files a/doc/manual/manual/geometry/zonotope.png 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 2fe1f1354..a0471be64 100644 --- a/doc/manual/manual/geometry/zonotope.rst +++ b/doc/manual/manual/geometry/zonotope.rst @@ -16,13 +16,13 @@ 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`. +The center is a :class:`Vector`, noted :math:`\mathbf{c}`, and the generators are stored in a :class:`Matrix`, noted :math:`\mathbf{A}`. Each column of the matrix corresponds to a generator. The resulting zonotope is: .. math:: - Z = z + A \cdot [-1,1]^m + Z = \mathbf{c} + \mathbf{A} \cdot [-1,1]^m Where :math:`m` is the number of generators (i.e., the number of columns of the matrix :math:`A`). @@ -37,24 +37,24 @@ It can be constructed in Codac as follows: .. code-tab:: py - z = Vector([1,2]) + c = Vector([1,2]) A = Matrix([[1,0,2],[0,1,1]]) - Z = Zonotope(z, A) + Z = Zonotope(c, A) .. code-tab:: c++ - Vector z({1,2}); + Vector c({1,2}); Matrix A({{1,0,2},{0,1,1}}); - Zonotope Z(z, A); + Zonotope Z(c, A); .. code-tab:: matlab - z = Vector({1,2}); + c = Vector({1,2}); A = Matrix({{1,0,2},{0,1,1}}); - Z = Zonotope(z,A); + Z = Zonotope(c,A); The resulting zonotope is represented in the figure below: diff --git a/doc/manual/manual/installation/cpp.rst b/doc/manual/manual/installation/cpp.rst index 9f0399474..4bf0ae41d 100644 --- a/doc/manual/manual/installation/cpp.rst +++ b/doc/manual/manual/installation/cpp.rst @@ -14,6 +14,10 @@ This page describes how to install the C++ library on Linux, Windows, and macOS. Linux Installation ------------------ +.. admonition:: Linux distributions + + Ubuntu or Debian-like systems are implied in most of the instructions (especially the ``sudo apt`` commands to install packages). + .. Install from packages (latest release) .. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -68,7 +72,7 @@ Steps These prerequisites can be installed directly with the command:: - sudo apt-get install -y g++ gcc cmake git + sudo apt-get install -y build-essential cmake git 2. **Install the IBEX dependency**: @@ -145,7 +149,7 @@ Steps cmake .. make - Before executing the example, you will have to launch the VIBes viewer. Then, + Before executing the example, you will have to launch the `VIBes viewer <../visualization/vibes.html>`_. Then, .. code-block:: bash @@ -163,7 +167,11 @@ Install from packages (latest release) Using MinGW ~~~~~~~~~~~ -Check https://community.chocolatey.org/packages/codac. +.. Check https://community.chocolatey.org/packages/codac. + +Install `Chocolatey package manager `_, run `choco install -y ibex cmake make qtcreator` in PowerShell and then download and extract *e.g.* ``codac_standalone_x64_mingw13.zip`` (for MinGW 13) from https://github.com/codac-team/codac/releases/latest, launch Qt Creator and choose Open Project, open ``example\CMakelists.txt``, ensure Desktop is selected and click Configure Project (might be hidden behind notifications at the bottom-right), wait 10 s then click on the big bottom-left green Run button, and finally check that the graphical output appears. + +Note that in order to obtain graphical outputs, you will have to download and run https://github.com/ENSTABretagneRobotics/VIBES/releases/latest/download/VIBes-viewer_x86.exe before running the project. Using Visual Studio @@ -171,7 +179,7 @@ Using Visual Studio Download and extract *e.g.* ``codac_standalone_x64_vc17.zip`` (for Visual Studio 2022) from https://github.com/codac-team/codac/releases/latest, open ``example\CMakelists.txt``, choose ``x64-Release`` configuration in Visual Studio (instead of ``x64-Debug``), double-click on ``main.cpp`` in the Solution Explorer and then click on the green Start button, and finally check that the graphical output appears. -Note that in order to obtain graphical outputs, you will have to download and run https://github.com/ENSTABretagneRobotics/VIBES/releases/download/0.2.3/VIBes-0.2.3-win32.exe before running the project. +Note that in order to obtain graphical outputs, you will have to download and run https://github.com/ENSTABretagneRobotics/VIBES/releases/latest/download/VIBes-viewer_x86.exe before running the project. Install from the sources (latest developments) @@ -190,15 +198,16 @@ Optionally, for Python binding (*e.g.* ``choco install python --version=3.10.4`` .. code-block:: bash - choco install doxygen.install --version=1.16 + choco install doxygen.install --version=1.16.1 choco install graphviz python -m pip install --upgrade pip pip install --upgrade wheel setuptools sphinx sphinx_rtd_theme furo sphinx-math-dollar sphinx_tabs breathe sphinx_togglebutton -The logic to follow will then be similar to Linux (note that for Visual Studio, commands such as ``make install`` need to be replaced with something similar to: +The logic to follow will then be similar to Linux. You might want to right-click in your desired folder and choose ``Git Bash Here`` to run the commands related to Git and compilation. Note that for Visual Studio, commands such as ``make`` and ``make install`` need to be replaced with something similar to: .. code-block:: bash + cmake --build . --config Release cmake --build . --config Release --target install .. warning:: @@ -221,7 +230,7 @@ Install `Homebrew package manager `_ and then build tools: brew install wget autoconf automake libtool brew install --cask cmake -Download and extract *e.g.* ``codac_standalone_arm64_monterey.zip`` from ``_, then in ``example`` folder run: +Download and extract *e.g.* ``codac_standalone_arm64_sonoma.zip`` from https://github.com/codac-team/codac/releases/latest, then in ``example`` folder run: .. code-block:: bash @@ -229,12 +238,12 @@ Download and extract *e.g.* ``codac_standalone_arm64_monterey.zip`` from ``_. +One of the graphical outputs supported in Codac is the `VIBes viewer <../visualization/vibes.html>`_. VIBes uses a client/server architecture and reads its configuration from the ``.vibes.json`` file located in your home directory (typically ``$HOME``). This file is automatically generated. If you launch MATLAB with ``sudo``, the ``HOME`` variable points to ``/root``, and VIBes will not find the ``.vibes.json`` file, causing the display to fail. diff --git a/doc/manual/manual/installation/python.rst b/doc/manual/manual/installation/python.rst index 0cdcfeb2c..f9df54ed4 100644 --- a/doc/manual/manual/installation/python.rst +++ b/doc/manual/manual/installation/python.rst @@ -8,7 +8,19 @@ Python installation In case you want to use Codac only with Python, then the installation procedure is simply: 1. **Dependencies**: Ensure you have a supported version of Python (>=3.8). -2. **Install via pip**: Use the following command to install Codac: +2. **Install pip**: For Ubuntu or Debian-like systems: + + .. code-block:: bash + + sudo apt install python3-pip + + For other systems, you can usually install pip with something similar to: + + .. code-block:: bash + + python -m pip install pip + +3. **Install via pip**: Use the following command to install Codac: .. code-block:: bash @@ -26,6 +38,7 @@ You can also update and test your Codac Python package: * **Testing the installation**: In order to verify that your installation is working properly in Python, you can run: .. code-block:: bash + pip install numpy sympy python -m unittest discover codac.tests @@ -34,7 +47,11 @@ Troubleshooting Depending on your configuration, you may encounter difficulties when installing the Codac Python package. Here are a few tips to help you. -.. admonition:: Debian Bookworm and possibly other configurations +.. admonition:: ``python`` vs ``python3``, ``pip`` vs ``pip3`` commands + + | You may need to use ``python3`` instead of ``python`` command depending on your configuration. Same applies to the ``pip`` command, which may be ``pip3`` in some cases. + +.. admonition:: Debian Bookworm, Ubuntu 24.04 and possibly other configurations | ``pip install ...`` or ``pip3 install ...`` commands may only work inside `virtual environments `_ or with ``--break-system-packages`` parameter, *e.g.*: @@ -42,6 +59,10 @@ Depending on your configuration, you may encounter difficulties when installing pip3 install --break-system-packages codac +.. admonition:: Python package dependencies for Ubuntu or Debian-like systems + + | If available, prefer ``sudo apt install python3-mypackage`` instead of ``pip install mypackage`` commands to install Python package dependencies for Ubuntu or Debian-like systems, to avoid breaking the system packages. However, not all the packages are available with ``apt``. + .. admonition:: macOS Big Sur and later (x86_64) | Python wheels provided for Catalina or earlier `may need to be installed with something similar to `_: diff --git a/doc/manual/manual/intervals/Interval_class.rst b/doc/manual/manual/intervals/Interval_class.rst index 48b24989a..39f49a7d7 100644 --- a/doc/manual/manual/intervals/Interval_class.rst +++ b/doc/manual/manual/intervals/Interval_class.rst @@ -243,6 +243,20 @@ All standard arithmetic operations are supported in Codac, both element-wise and z = x / y % [1, 3] +Interval comparisons +-------------------- + +The comparison operators ``<`` and ``>`` are also overloaded for intervals. +They return a :ref:`BoolInterval value ` rather than a classical Boolean, in order to encode the set of feasible truth values when the exact operands are only known within interval bounds. + +For instance, ``x < y`` returns: + +- ``BoolInterval.EMPTY`` if one of the two intervals is empty; +- ``BoolInterval.TRUE`` if :math:`\forall a\in[x],\forall b\in[y],\ a` + + - :ref:`PEIBOS for ODE integration ` + +How it works +------------ + +By default, Codac uses only one thread for its computations. +A getter ``nb_threads`` and a setter ``set_nb_threads`` are available to get and change the number of threads. + +In addition, the function ``max_threads`` provides the maximum number of threads usable on the machine. + +To use a given number of threads in a function, the setter just needs to be called before running the desired function. An example is provided below + +.. tabs:: + + .. code-tab:: py + + set_nb_threads(max_threads()) # using as many threads as possible + PEIBOS(...) + + + .. code-tab:: c++ + + set_nb_threads(max_threads()); // using as many threads as possible + PEIBOS(...); + + .. code-tab:: matlab + + set_nb_threads(max_threads()); % using as many threads as possible + PEIBOS(...); + \ No newline at end of file diff --git a/doc/manual/manual/visualization/functions.rst b/doc/manual/manual/visualization/functions.rst index afa66ccc2..b3efcc917 100644 --- a/doc/manual/manual/visualization/functions.rst +++ b/doc/manual/manual/visualization/functions.rst @@ -56,9 +56,9 @@ for more information. .. code-tab:: matlab - DefaultFigure().draw_box(IntervalVector({{2.2,2.5},{2.2,2.5}}),Color().red()); % Red edge - DefaultFigure().draw_box(IntervalVector({{2.2,2.5},{2.2,2.5}}),StyleProperties({Color().red()})); % Red edge - DefaultFigure().draw_box(IntervalVector({{2.2,2.5},{2.2,2.5}}),StyleProperties({Color().red(),Color().blue()})); % Red edge and blue fill + fig.draw_box(IntervalVector({{2.2,2.5},{2.2,2.5}}),Color().red()); % Red edge + fig.draw_box(IntervalVector({{2.2,2.5},{2.2,2.5}}),StyleProperties({Color().red()})); % Red edge + fig.draw_box(IntervalVector({{2.2,2.5},{2.2,2.5}}),StyleProperties({Color().red(),Color().blue()})); % Red edge and blue fill A line style and a layer can be added to the StyleProperties object, see :ref:`the dedicated page ` for more information. diff --git a/examples/06_graphics_3D/main.cpp b/examples/06_graphics_3D/main.cpp index 3736d7950..a77f85495 100644 --- a/examples/06_graphics_3D/main.cpp +++ b/examples/06_graphics_3D/main.cpp @@ -38,7 +38,9 @@ int main() Figure3D fig_examples("3D examples"); - fig_examples.draw_axes(1.0); + fig_examples.draw_axes(); + fig_examples.draw_axes(0.5); + fig_examples.draw_axes(2.0,{0.5,0.5,0.5}); fig_examples.draw_triangle({1,0,0},{0,1,0},{0,0,1},{ Color::dark_green(0.5), "triangle1" }); fig_examples.draw_triangle({2,0,0},{{-1,0,0},{0,1,1},{0,0,-1}}, {1,0,0},{0,1,0},{0,0,1},Color::purple(0.5)); diff --git a/examples/06_graphics_3D/main.py b/examples/06_graphics_3D/main.py index af87fcd5d..4ebf27354 100644 --- a/examples/06_graphics_3D/main.py +++ b/examples/06_graphics_3D/main.py @@ -22,7 +22,9 @@ fig_sep.draw_paving(p_sep) fig_examples = Figure3D("3D examples") -fig_examples.draw_axes(1.0) +fig_examples.draw_axes() +fig_examples.draw_axes(0.5) +fig_examples.draw_axes(2.0,[0.5,0.5,0.5]) fig_examples.draw_triangle([1,0,0],[0,1,0],[0,0,1],StyleProperties(Color.dark_green(0.5),"triangle")) fig_examples.draw_triangle([2,0,0], Matrix([[2,0,0],[-1,0,0],[0,1,1]]), diff --git a/examples/11_peibos/main.m b/examples/11_peibos/main.m index be52dbcd8..c5fe7acfc 100644 --- a/examples/11_peibos/main.m +++ b/examples/11_peibos/main.m @@ -4,6 +4,9 @@ [currentDir, ~, ~] = fileparts(currentFilePath); cd(currentDir); +%% +set_nb_threads(max_threads()); + %% % 2D example of the PEIBOS algorithm diff --git a/packages/choco/codac/codac.nuspec b/packages/choco/codac/codac.nuspec index 4dcf8ecec..b4d71e6bd 100644 --- a/packages/choco/codac/codac.nuspec +++ b/packages/choco/codac/codac.nuspec @@ -18,8 +18,8 @@ Codac is a library providing tools for constraint programming over reals, trajectories and sets. ## Quick start -- After [installing Chocolatey](https://chocolatey.org/install), run `choco install -y codac cmake make qtcreator` in PowerShell and then launch Qt Creator and choose Open Project, open `test_codac\CMakelists.txt` (get the folder from https://www.ensta-bretagne.fr/lebars/Share/test_codac.zip), ensure Desktop is selected and click Configure Project, wait 10 s then click on the big bottom green Run button, finally check that "My first tube:Tube [0, 10]" appears. -- Optionally, download and run https://github.com/ENSTABretagneRobotics/VIBES/releases/download/0.2.3/VIBes-0.2.3-win32.exe before running the project, and check that a tube appears in VIBes window. +- After [installing Chocolatey](https://chocolatey.org/install), run `choco install -y codac cmake make qtcreator` in PowerShell and then launch Qt Creator and choose Open Project, open `test_codac\CMakelists.txt` (get the folder from https://www.ensta-bretagne.fr/lebars/Share/test_codac.zip), ensure Desktop is selected and click Configure Project (might be hidden behind notifications at the bottom-right), wait 10 s then click on the big bottom-left green Run button, finally check that "Computation time: XXXs, XXX boxes" appears. +- Download and run https://github.com/ENSTABretagneRobotics/VIBES/releases/latest/download/VIBes-viewer_x86.exe before running the project, and check that the graphical output appears. ## Troubleshooting - Check that all the packages and their dependencies were installed, if one failed (e.g. due to network-related errors) try to reinstall it using `--force` or try a previous version... diff --git a/python/src/core/CMakeLists.txt b/python/src/core/CMakeLists.txt index d569f74bd..9d83a210d 100644 --- a/python/src/core/CMakeLists.txt +++ b/python/src/core/CMakeLists.txt @@ -110,6 +110,7 @@ separators/codac2_py_SepInter.cpp separators/codac2_py_SepInverse.cpp separators/codac2_py_SepNot.cpp + separators/codac2_py_SepPolarCart.cpp separators/codac2_py_SepPolygon.cpp separators/codac2_py_SepProj.cpp separators/codac2_py_SepQInter.cpp diff --git a/python/src/core/codac2_py_core.cpp b/python/src/core/codac2_py_core.cpp index 370630623..afe9b27f3 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -148,6 +148,7 @@ void export_SepCtcPair(py::module& m, py::class_& pysep); void export_SepInter(py::module& m, py::class_& sep); void export_SepInverse(py::module& m, py::class_& sep); void export_SepNot(py::module& m, py::class_& sep); +void export_SepPolarCart(py::module& m, py::class_& pysep); void export_SepPolygon(py::module& m, py::class_& sep); void export_SepProj(py::module& m, py::class_& sep); void export_SepQInter(py::module& m, py::class_& sep); @@ -331,6 +332,7 @@ PYBIND11_MODULE(_core, m) export_SepInter(m,py_sep); export_SepInverse(m,py_sep); export_SepNot(m,py_sep); + export_SepPolarCart(m,py_sep); export_SepPolygon(m,py_sep); export_SepProj(m,py_sep); export_SepQInter(m,py_sep); @@ -368,4 +370,4 @@ PYBIND11_MODULE(_core, m) // Extension > sympy auto ms = m.def_submodule("_sympy"); // to keep a dedicated namespace export_sympy(ms); -} \ No newline at end of file +} diff --git a/python/src/core/contractors/codac2_py_CtcWrapper.cpp b/python/src/core/contractors/codac2_py_CtcWrapper.cpp index d07306bad..4d3b6155e 100644 --- a/python/src/core/contractors/codac2_py_CtcWrapper.cpp +++ b/python/src/core/contractors/codac2_py_CtcWrapper.cpp @@ -34,6 +34,18 @@ void export_CtcWrapper(py::module& m, py::class_,pyCtcIn CONTRACT_METHODS(exported_ctcwrapper_intervalvector, CtcWrapper, VOID_CTCWRAPPER_YX_CONTRACT_X_REF_CONST) + py::class_> exported_ctcwrapper_parallelepiped(m, "CtcWrapper_Parallelepiped", pyctc, CTCWRAPPER_MAIN); + exported_ctcwrapper_parallelepiped + + .def(py::init(), + CTCWRAPPER_YX_CTCWRAPPER_CONST_Y_REF, + "y"_a) + + ; + + CONTRACT_METHODS(exported_ctcwrapper_parallelepiped, CtcWrapper, + VOID_CTCWRAPPER_YX_CONTRACT_X_REF_CONST) + py::class_> exported_ctcwrapper_pavingout(m, "CtcWrapper_PavingOut", pyctc, CTCWRAPPER_MAIN); exported_ctcwrapper_pavingout diff --git a/python/src/core/domains/interval/codac2_py_Interval.cpp b/python/src/core/domains/interval/codac2_py_Interval.cpp index 145e9ad4c..2e92ca484 100644 --- a/python/src/core/domains/interval/codac2_py_Interval.cpp +++ b/python/src/core/domains/interval/codac2_py_Interval.cpp @@ -76,6 +76,14 @@ py::class_ export_Interval(py::module& m) BOOL_INTERVAL_OPERATORNEQ_CONST_INTERVAL_REF_CONST, "x"_a) + .def(py::self < py::self, + BOOLINTERVAL_INTERVAL_OPERATOR__CONST_INTERVAL_REF_CONST, + "x"_a) + + .def(py::self > py::self, + BOOLINTERVAL_INTERVAL_OPERATOR_CONST_INTERVAL_REF_CONST, + "x"_a) + .def("lb", &Interval::lb, DOUBLE_INTERVAL_LB_CONST) diff --git a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp index 0cfd384bd..d3181936b 100644 --- a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp +++ b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp @@ -31,6 +31,10 @@ void export_Parallelepiped(py::module& m) PARALLELEPIPED_PARALLELEPIPED_CONST_VECTOR_REF_CONST_MATRIX_REF, "z"_a, "A"_a) + .def_static("empty", &Parallelepiped::empty, + STATIC_PARALLELEPIPED_PARALLELEPIPED_EMPTY_INDEX, + "n"_a) + .def("vertices", &Parallelepiped::vertices, VECTOR_VECTOR_PARALLELEPIPED_VERTICES_CONST) @@ -42,5 +46,15 @@ void export_Parallelepiped(py::module& m) BOOLINTERVAL_PARALLELEPIPED_IS_SUPERSET_CONST_INTERVALVECTOR_REF_CONST, "x"_a) + .def("__and__", &Parallelepiped::operator&, + PARALLELEPIPED_PARALLELEPIPED_OPERATORINTER_CONST_INTERVALVECTOR_REF_CONST, + "x"_a) + + .def("__rand__", [](const Parallelepiped& p, const IntervalVector& x) + { + return x & p; + }, + INTERVALVECTOR_OPERATORINTER_CONST_INTERVALVECTOR_REF_CONST_PARALLELEPIPED_REF, + "x"_a) ; } diff --git a/python/src/core/domains/zonotope/codac2_py_Zonotope.cpp b/python/src/core/domains/zonotope/codac2_py_Zonotope.cpp index 7b0977919..89f25b699 100644 --- a/python/src/core/domains/zonotope/codac2_py_Zonotope.cpp +++ b/python/src/core/domains/zonotope/codac2_py_Zonotope.cpp @@ -29,6 +29,16 @@ void export_Zonotope(py::module& m) ZONOTOPE_ZONOTOPE_CONST_VECTOR_REF_CONST_MATRIX_REF, "z"_a, "A"_a) + .def_static("empty", &Zonotope::empty, + STATIC_ZONOTOPE_ZONOTOPE_EMPTY_INDEX, + "n"_a) + + .def("is_empty", &Zonotope::is_empty, + BOOL_ZONOTOPE_IS_EMPTY_CONST) + + .def("set_empty", &Zonotope::set_empty, + VOID_ZONOTOPE_SET_EMPTY) + .def("box", &Zonotope::box, INTERVALVECTOR_ZONOTOPE_BOX_CONST) @@ -39,8 +49,20 @@ void export_Zonotope(py::module& m) ZONOTOPE_ZONOTOPE_PROJ_CONST_VECTOR_INDEX_REF_CONST, "indices"_a) - .def_readwrite("z", &Zonotope::z, - VECTOR_ZONOTOPE_Z) + .def("__add__", &Zonotope::operator+, + ZONOTOPE_ZONOTOPE_OPERATORPLUS_CONST_ZONOTOPE_REF, + py::is_operator()) + + .def("__repr__", [](const Zonotope& z) { + std::ostringstream stream; + stream << z; + return string(stream.str()); + }, + OSTREAM_REF_OPERATOROUT_OSTREAM_REF_CONST_ZONOTOPE_REF) + + + .def_readwrite("c", &Zonotope::c, + VECTOR_ZONOTOPE_C) .def_readwrite("A", &Zonotope::A, MATRIX_ZONOTOPE_A) diff --git a/python/src/core/separators/codac2_py_SepPolarCart.cpp b/python/src/core/separators/codac2_py_SepPolarCart.cpp new file mode 100644 index 000000000..3ad0fd89d --- /dev/null +++ b/python/src/core/separators/codac2_py_SepPolarCart.cpp @@ -0,0 +1,39 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include +#include "codac2_py_Sep.h" +#include "codac2_py_SepPolarCart_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_SepPolarCart(py::module& m, py::class_& pysep) +{ + py::class_ exported(m, "SepPolarCart", pysep, SEPPOLARCART_MAIN); + exported + + .def(py::init([](const SepBase& s1) { + return std::make_unique(s1.copy()); + }), + SEPPOLARCART_SEPPOLARCART_CONST_S_REF, + "s1"_a) + + .def("separate", &SepPolarCart::separate, + BOXPAIR_SEPPOLARCART_SEPARATE_CONST_INTERVALVECTOR_REF_CONST, + "x"_a) + + ; +} \ No newline at end of file diff --git a/python/src/graphics/figures/codac2_py_Figure3D.cpp b/python/src/graphics/figures/codac2_py_Figure3D.cpp index 24dbf165f..620e9ecdc 100644 --- a/python/src/graphics/figures/codac2_py_Figure3D.cpp +++ b/python/src/graphics/figures/codac2_py_Figure3D.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "codac2_py_Figure3D_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): #include "codac2_py_matlab.h" @@ -35,8 +36,8 @@ void export_Figure3D(py::module& m) CONST_STRING_REF_FIGURE3D_NAME_CONST) .def("draw_axes", &Figure3D::draw_axes, - VOID_FIGURE3D_DRAW_AXES_DOUBLE, - "size"_a=1.0) + VOID_FIGURE3D_DRAW_AXES_DOUBLE_CONST_VECTOR_REF, + "size"_a=1.0, "origin"_a=Vector::zero(3)) // Geometric shapes .def("draw_triangle", (void(Figure3D::*)(const Vector &c, const Matrix &A, const Vector &p1, const Vector &p2, const Vector &p3, const StyleProperties &s))&Figure3D::draw_triangle, diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index d06db4fd8..db02e9526 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -241,6 +241,8 @@ ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepInter.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepInverse.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepNot.h + ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolarCart.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolarCart.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolygon.cpp ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolygon.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepProj.cpp diff --git a/src/core/contractors/codac2_CtcWrapper.h b/src/core/contractors/codac2_CtcWrapper.h index 83b7745fd..278e6ef93 100644 --- a/src/core/contractors/codac2_CtcWrapper.h +++ b/src/core/contractors/codac2_CtcWrapper.h @@ -29,7 +29,7 @@ namespace codac2 void contract(X& x) const { assert_release(x.size() == this->size()); - x = _y & x; + x = x & _y; } protected: diff --git a/src/core/domains/interval/codac2_Interval.h b/src/core/domains/interval/codac2_Interval.h index 34e082421..09dca1ac5 100644 --- a/src/core/domains/interval/codac2_Interval.h +++ b/src/core/domains/interval/codac2_Interval.h @@ -21,6 +21,7 @@ #include "codac2_Domain.h" #include "codac2_assert.h" #include "codac2_TypeInfo.h" +#include "codac2_BoolInterval.h" namespace codac2 { @@ -160,6 +161,40 @@ namespace codac2 */ bool operator!=(const Interval& x) const; + /** + * \brief Comparison (strict less-than) between this and x + * + * The returned ``BoolInterval`` encloses the truth value of + * \f$ts\f$ for \f$t\in[\mathrm{this}]\f$ and \f$s\in[x]\f$. + * + * \note Returns: + * - ``BoolInterval::EMPTY`` if this or x is empty + * - ``BoolInterval::TRUE`` iff \f$\mathrm{lb}([\mathrm{this}])>\mathrm{ub}([x])\f$ + * - ``BoolInterval::FALSE`` iff \f$\mathrm{lb}([x])\geq\mathrm{ub}([\mathrm{this}])\f$ + * - ``BoolInterval::UNKNOWN`` otherwise + * + * \param x interval to be compared with + * \return interval Boolean result + */ + BoolInterval operator>(const Interval& x) const; + /** * \brief Returns the lower bound of this * diff --git a/src/core/domains/interval/codac2_Interval_impl.h b/src/core/domains/interval/codac2_Interval_impl.h index bcdcef9b8..8e2bd0cfc 100644 --- a/src/core/domains/interval/codac2_Interval_impl.h +++ b/src/core/domains/interval/codac2_Interval_impl.h @@ -107,6 +107,22 @@ namespace codac2 return !(*this == x); } + inline BoolInterval Interval::operator<(const Interval& x) const + { + if(is_empty() || x.is_empty()) + return BoolInterval::EMPTY; + if(this->ub() < x.lb()) + return BoolInterval::TRUE; + if(x.ub() <= this->lb()) + return BoolInterval::FALSE; + return BoolInterval::UNKNOWN; + } + + inline BoolInterval Interval::operator>(const Interval& x) const + { + return x < *this; + } + inline double Interval::lb() const { return gaol::interval::left(); diff --git a/src/core/domains/paving/codac2_Paving.cpp b/src/core/domains/paving/codac2_Paving.cpp index f7eb538d8..c14c2142f 100644 --- a/src/core/domains/paving/codac2_Paving.cpp +++ b/src/core/domains/paving/codac2_Paving.cpp @@ -74,6 +74,11 @@ namespace codac2 return list(); }; + IntervalVector operator&(const IntervalVector& x, const PavingOut& p) + { + return p & x; + }; + // PavingInOut class diff --git a/src/core/domains/paving/codac2_Paving.h b/src/core/domains/paving/codac2_Paving.h index 3bd70a2a7..25f5c8f39 100644 --- a/src/core/domains/paving/codac2_Paving.h +++ b/src/core/domains/paving/codac2_Paving.h @@ -173,6 +173,8 @@ namespace codac2 static const NodeValue_ outer, outer_complem; }; + IntervalVector operator&(const IntervalVector& x, const PavingOut& p); + class PavingInOut; using PavingInOut_Node = PavingNode; diff --git a/src/core/domains/zonotope/codac2_Parallelepiped.cpp b/src/core/domains/zonotope/codac2_Parallelepiped.cpp index ef75b58fa..79594db6f 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped.cpp +++ b/src/core/domains/zonotope/codac2_Parallelepiped.cpp @@ -10,55 +10,106 @@ #include "codac2_Parallelepiped.h" #include "codac2_inversion.h" -using namespace codac2; +using namespace std; -Parallelepiped::Parallelepiped(const Vector& z_, const Matrix& A_) - : Zonotope(z_, A_) +namespace codac2 { - assert_release(A.cols() <= z.size() && "too many vectors, you are describing a zonotope"); -} + Parallelepiped::Parallelepiped(const Vector& c_, const Matrix& A_) + : Zonotope(c_, A_) + { + assert_release(A.cols() <= c.size() && "too many vectors, you are describing a zonotope"); + } -void generate_vertices(Index i, Index n, const Vector& z, const Matrix& A, std::vector& L_v) -{ - if (i == n) + Parallelepiped Parallelepiped::empty(Index n) { - L_v.push_back(z); + Parallelepiped p(Vector::Constant(n, std::numeric_limits::quiet_NaN()), Matrix::Identity(n,n)); + p.empty_flag = true; + return p; } - else if (i& L_v) { - generate_vertices(i+1, n, z + A.col(i), A, L_v); - generate_vertices(i+1, n, z - A.col(i), A, L_v); + if (i == n) + { + L_v.push_back(c); + } + else if (i Parallelepiped::vertices() const -{ - std::vector L_v; - generate_vertices(0, z.size(),z,A,L_v); - return L_v; -} + std::vector Parallelepiped::vertices() const + { + std::vector L_v; + generate_vertices(0, size(),c,A,L_v); + return L_v; + } -BoolInterval Parallelepiped::contains(const Vector& v) const -{ - return is_superset(v.template cast()); -} + BoolInterval Parallelepiped::contains(const Vector& v) const + { + return is_superset(v.template cast()); + } -BoolInterval Parallelepiped::is_superset(const IntervalVector& x) const -{ - assert_release(A.rows() == A.cols() && "Matrix A must be square to check containment."); - assert_release(x.size() == z.size() && "Point dimension must match parallelepiped dimension."); + BoolInterval Parallelepiped::is_superset(const IntervalVector& x) const + { + assert_release(A.rows() == A.cols() && "Matrix A must be square to check containment."); + assert_release(x.size() == size() && "Point dimension must match parallelepiped dimension."); - IntervalVector B = inverse_enclosure(A)*(x - z); - IntervalVector IV = IntervalVector::constant(A.cols(),{-1,1}); + IntervalVector B = inverse_enclosure(A)*(x - c); + IntervalVector IV = IntervalVector::constant(A.cols(),{-1,1}); - if (!(B.intersects(IV))) - return BoolInterval::FALSE; + if (!(B.intersects(IV))) + return BoolInterval::FALSE; - else - { - if (B.is_subset(IV)) - return BoolInterval::TRUE; else - return BoolInterval::UNKNOWN; + { + if (B.is_subset(IV)) + return BoolInterval::TRUE; + else + return BoolInterval::UNKNOWN; + } } -} \ No newline at end of file + + Parallelepiped Parallelepiped::operator&(const IntervalVector& x) const + { + assert_release(A.rows() == A.cols() && "Matrix A must be square to compute intersection."); + assert_release(x.size() == size() && "Box dimension must match parallelepiped dimension."); + + if (is_empty() || x.is_empty()) + return Parallelepiped::empty(size()); + + IntervalVector x_p = inverse_enclosure(A)*(x-c); + x_p &= IntervalVector::constant(c.size(),Interval(-1,1)); + + if (x_p.is_empty()) + return Parallelepiped::empty(size()); + + Vector c_i = c + A*x_p.mid(); + + Matrix A_i (A); + for (int i = 0; i < c.size(); ++i) + A_i.col(i) *= x_p[i].rad(); + + return Parallelepiped(c_i,A_i); + } + + IntervalVector operator&(const IntervalVector& x, const Parallelepiped& p) + { + assert_release(p.A.rows() == p.A.cols() && "Matrix A must be square to compute intersection."); + assert_release(x.size() == p.size() && "Box dimension must match parallelepiped dimension."); + + if (p.is_empty() || x.is_empty()) + return IntervalVector::empty(x.size()); + + IntervalVector x_p = inverse_enclosure(p.A)*(x-p.c); + x_p &= IntervalVector::constant(p.c.size(),Interval(-1,1)); + + if (x_p.is_empty()) + return IntervalVector::empty(x.size()); + + IntervalVector x_i = p.c + p.A*x_p; + return x_i & x; + } +} diff --git a/src/core/domains/zonotope/codac2_Parallelepiped.h b/src/core/domains/zonotope/codac2_Parallelepiped.h index 19b6d9338..9dfbfd81a 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped.h @@ -21,9 +21,9 @@ namespace codac2 /** * \class Parallelepiped * - * \brief Class representing a parallelepiped \f$\mathbf{z} + \mathbf{A}\cdot[-1,1]^m\f$ + * \brief Class representing a parallelepiped \f$\mathbf{c} + \mathbf{A}\cdot[-1,1]^m\f$ * - * This class represents a parallelepiped in n-dimensional space, defined by a center point \f$\mathbf{z}\f$ and a shape matrix \f$\mathbf{A}\f$. + * This class represents a parallelepiped in n-dimensional space, defined by a center point \f$\mathbf{c}\f$ and a shape matrix \f$\mathbf{A}\f$. * * A parallelepiped is a special case of a zonotope where the shape matrix \f$\mathbf{A}\f$ has \f$m\f$ columns with \f$m \leqslant n\f$. */ @@ -34,10 +34,18 @@ namespace codac2 /** * \brief Constructs a n-parallelepiped object with a given center and shape matrix * - * \param z Center of the parallelepiped (n-dimensional vector) + * \param c Center of the parallelepiped (n-dimensional vector) * \param A Shape matrix of the parallelepiped (\f$n\times m\f$ matrix with \f$m \leqslant n\f$) */ - Parallelepiped(const Vector& z, const Matrix& A); + Parallelepiped(const Vector& c, const Matrix& A); + + /** + * \brief Constructs an empty n-parallelepiped + * + * \param n Dimension of the parallelepiped + * \return A new Parallelepiped object representing an empty parallelepiped in n-dimensional space + */ + static Parallelepiped empty(Index n); /** * \brief Computes the vertices of the parallelepiped @@ -64,5 +72,23 @@ namespace codac2 */ BoolInterval is_superset(const IntervalVector& x) const; + /** + * \brief Computes the intersection of the parallelepiped with a given box. The matrix A has to be square and invertible. + * + * \param x The box to intersect with + * + * \return A new Parallelepiped representing the intersection of the original parallelepiped with the box. + */ + Parallelepiped operator&(const IntervalVector& x) const; }; + + /** + * \brief Computes the intersection of a box with a parallelepiped. The matrix A has to be square and invertible. + * + * \param x The box to intersect with + * \param p The parallelepiped to intersect with + * + * \return A new IntervalVector representing the intersection of the box with the parallelepiped. + */ + IntervalVector operator&(const IntervalVector& x, const Parallelepiped& p); } diff --git a/src/core/domains/zonotope/codac2_Zonotope.cpp b/src/core/domains/zonotope/codac2_Zonotope.cpp index 9091e7f4d..0e1844d43 100644 --- a/src/core/domains/zonotope/codac2_Zonotope.cpp +++ b/src/core/domains/zonotope/codac2_Zonotope.cpp @@ -9,33 +9,76 @@ #include "codac2_Zonotope.h" -using namespace codac2; - -Zonotope::Zonotope(const Vector& z_, const Matrix& A_) - : z(z_), A(A_) +namespace codac2 { - assert_release(z.size() == A.rows()); -} + Zonotope::Zonotope(const Vector& c_, const Matrix& A_) + : c(c_), A(A_) + { + assert_release(c.size() == A.rows()); + } -IntervalVector Zonotope::box() const -{ - return z + (A.template cast())*IntervalVector::constant(A.cols(),{-1,1}); -} + Zonotope Zonotope::empty(Index n) + { + Zonotope z(Vector::Constant(n, std::numeric_limits::quiet_NaN()), Matrix::Identity(n,n)); + z.empty_flag = true; + return z; + } -Zonotope Zonotope::proj(const std::vector& indices) const -{ - assert_release(*std::min_element(indices.begin(), indices.end()) >= 0 && "indices out of range"); - assert_release(*std::max_element(indices.begin(), indices.end()) <= z.size() && "indices out of range"); + bool Zonotope::is_empty() const + { + return empty_flag; + } - Matrix A_cropped (indices.size(), A.cols()); - Vector z_cropped (indices.size()); + void Zonotope::set_empty() + { + empty_flag = true; + c = Vector::Constant(c.size(), std::numeric_limits::quiet_NaN()); + } - for (size_t i = 0; i < indices.size(); ++i) + IntervalVector Zonotope::box() const { - A_cropped.row(i) = A.row(indices[i]); - z_cropped[i] = z[indices[i]]; + return c + (A.template cast())*IntervalVector::constant(A.cols(),{-1,1}); } - return Zonotope(z_cropped, A_cropped); -} + Zonotope Zonotope::proj(const std::vector& indices) const + { + assert_release(*std::min_element(indices.begin(), indices.end()) >= 0 && "indices out of range"); + assert_release(*std::max_element(indices.begin(), indices.end()) <= c.size() && "indices out of range"); + Matrix A_cropped (indices.size(), A.cols()); + Vector c_cropped (indices.size()); + + for (size_t i = 0; i < indices.size(); ++i) + { + A_cropped.row(i) = A.row(indices[i]); + c_cropped[i] = c[indices[i]]; + } + + return Zonotope(c_cropped, A_cropped); + } + + Zonotope Zonotope::operator+(const Zonotope& zonotope) + { + assert_release(c.size() == zonotope.c.size() && "Zonotopes must have the same dimension"); + + Vector c_sum = c + zonotope.c; + Matrix A_sum (A.rows(), A.cols() + zonotope.A.cols()); + + A_sum << A, zonotope.A; + + return Zonotope(c_sum, A_sum); + } + + std::ostream& operator<<(std::ostream& str, const Zonotope& z) + { + str << "{ "; + + if (z.is_empty()) + str << "empty " << (z.size()) << "-d zonotope"; + else + str << "Zonotope center: \n" << z.c.transpose() << "\nZonotope shape matrix: \n" << z.A; + + str << " }"; + return str; + } +} diff --git a/src/core/domains/zonotope/codac2_Zonotope.h b/src/core/domains/zonotope/codac2_Zonotope.h index 863a2681b..bd9e6e043 100644 --- a/src/core/domains/zonotope/codac2_Zonotope.h +++ b/src/core/domains/zonotope/codac2_Zonotope.h @@ -18,11 +18,11 @@ namespace codac2 { /** * \class Zonotope - * \brief Class representing a zonotope \f$\mathbf{z} + \mathbf{A}\cdot[-1,1]^m\f$ + * \brief Class representing a zonotope \f$\mathbf{c} + \mathbf{A}\cdot[-1,1]^m\f$ * - * This class represents a zonotope in n-dimensional space, defined by a center point \f$\mathbf{z}\f$ and a shape matrix \f$\mathbf{A}\f$. + * This class represents a zonotope in n-dimensional space, defined by a center point \f$\mathbf{c}\f$ and a shape matrix \f$\mathbf{A}\f$. * - * The vector \f$\mathbf{z}\f$ and each column of the matrix \f$\mathbf{A}\f$ must have the same dimension \f$n\f$, but the matrix \f$\mathbf{A}\f$ can have any number of columns \f$m\f$. + * The vector \f$\mathbf{c}\f$ and each column of the matrix \f$\mathbf{A}\f$ must have the same dimension \f$n\f$, but the matrix \f$\mathbf{A}\f$ can have any number of columns \f$m\f$. */ class Zonotope { @@ -31,10 +31,39 @@ namespace codac2 /** * \brief Constructs a n-zonotope object with a given center and shape matrix * - * \param z Center of the zonotope (n-dimensional vector) + * \param c Center of the zonotope (n-dimensional vector) * \param A Shape matrix of the zonotope (\f$n\times m\f$ matrix) */ - Zonotope(const Vector& z, const Matrix& A); + Zonotope(const Vector& c, const Matrix& A); + + /** + * \brief Constructs an empty n-zonotope + * + * \param n Dimension of the zonotope + * + * \return A new Zonotope object representing an empty zonotope in n-dimensional space + */ + static Zonotope empty(Index n); + + /** + * \brief Checks if the zonotope is empty + * + * \return True if the zonotope is empty, false otherwise + */ + bool is_empty() const; + + /** + * \brief Sets the zonotope to be empty + */ + void set_empty(); + + /** + * \brief Outputs the size of the zonotope + * + * \return The size of the zonotope, defined as the size of its center vector + */ + + Index size() const { return c.size(); } /** * \brief Computes the axis-aligned bounding box of the zonotope @@ -52,15 +81,37 @@ namespace codac2 */ Zonotope proj(const std::vector& indices) const; + /** + * \brief Computes the Minkowski sum of this Zonotope with another Zonotope + * + * \param zonotope The other Zonotope to add to this one + * + * \return A new Zonotope object representing the Minkowski sum of the two Zonotopes + */ + Zonotope operator+(const Zonotope& zonotope); + /** * \brief Center of the zonotope */ - Vector z; + Vector c; /** * \brief Shape matrix of the zonotope */ Matrix A; - + + /** + * \brief Flag indicating whether the zonotope is empty + */ + bool empty_flag = false; }; + + /** + * \brief Stream output operator for ``Zonotope``. + * + * \param str Output stream. + * \param z The zonotope to print. + * \return The output stream with the zonotope information. + */ + std::ostream& operator<<(std::ostream& str, const Zonotope& z); } diff --git a/src/core/functions/analytic/codac2_AnalyticFunction.h b/src/core/functions/analytic/codac2_AnalyticFunction.h index 337a64d03..0fe6fd645 100644 --- a/src/core/functions/analytic/codac2_AnalyticFunction.h +++ b/src/core/functions/analytic/codac2_AnalyticFunction.h @@ -220,17 +220,17 @@ namespace codac2 "Parallelepiped evaluation requires at least one input."); IntervalVector Y = this->eval(((typename Wrapper::Domain)(x)).mid()...); - Vector z = Y.mid(); + Vector c = Y.mid(); Matrix A = this->diff(((typename Wrapper::Domain)(x)).mid()...).mid(); // Maximum error computation - double rho = error_peibos(Y, z, this->diff(x...), A, cart_prod(x...)); + double rho = error_peibos(Y, c, this->diff(x...), A, cart_prod(x...)); // Inflation of the parallelepiped Matrix A_inf = inflate_flat_parallelepiped(A, (cart_prod(x...).template cast()).rad(), rho); - return Parallelepiped(z, A_inf); + return Parallelepiped(c, A_inf); } Index output_size() const diff --git a/src/core/matrices/codac2_IntvFullPivLU.cpp b/src/core/matrices/codac2_IntvFullPivLU.cpp index a8e41431a..4820cf9a8 100644 --- a/src/core/matrices/codac2_IntvFullPivLU.cpp +++ b/src/core/matrices/codac2_IntvFullPivLU.cpp @@ -410,8 +410,9 @@ IntervalMatrix IntvFullPivLU::cokernel() const { IntervalMatrix IntvFullPivLU::solve(const IntervalMatrix &rhs) const { assert_release(rhs.rows()==matrixLU_.rows()); IntervalMatrix prhs = _LU.permutationP()*rhs; + Index dim = std::min(matrixLU_.rows(),matrixLU_.cols()); /* inverse L */ - for (Index r = 0; r=0;r--) { for (Index r1=r+1;r1(); for (Index r = dim-1;r>=0;r--) { IntervalRow rw=U.row(r); diff --git a/src/core/separators/codac2_SepPolarCart.cpp b/src/core/separators/codac2_SepPolarCart.cpp index dbb726d12..ac011ccc3 100644 --- a/src/core/separators/codac2_SepPolarCart.cpp +++ b/src/core/separators/codac2_SepPolarCart.cpp @@ -1,74 +1,115 @@ /** - * codac2_SepPolarCart.cpp + * \file codac2_SepPolarCart.cpp * ---------------------------------------------------------------------------- * \date 2026 - * \author Benoît Desrochers, (Simon Rohou) + * \author Simon Rohou, from a former implementation of Benoît Desrochers * \copyright Copyright 2026 Codac Team * \license GNU Lesser General Public License (LGPL) */ #include "codac2_SepPolarCart.h" #include "codac2_CtcPolar.h" +#include "codac2_cart_prod.h" +#include "codac2_hull.h" using namespace codac2; -BoxPair SepPolarCart::separate(const IntervalVector& x) const +/** + * Contracts x with the Cartesian projection of the complement of a polar box. + * + * For a polar box [rho]x[theta], the complement is represented as the union + * of four polar boxes: + * - rho <= rho.lb(), with theta in [theta] + * - rho >= rho.ub(), with theta in [theta] + * - theta before theta.lb(), with rho in R+ + * - theta after theta.ub(), with rho in R+ + * + * The angular complement is computed in the 2*pi-wide window centered on [theta]. + */ +IntervalVector contract_polar_box_complement( + const CtcPolar& ctc_polar, + const IntervalVector& x, + const IntervalVector& x_pol) { assert_release(x.size() == 2); + assert_release(x_pol.size() == 2); - CtcPolar ctc_polar; - IntervalVector x_cart(2), x_pol(x); + if(x.is_empty()) + return x; - ctc_polar.contract(x_cart[0], x_cart[1], x_pol[0], x_pol[1]); - auto x_cart_sep = _sep.front()->separate(x_cart); - BoxPair x_polar_sep { x, x }; - ctc_polar.contract(x_cart_sep.inner[0], x_cart_sep.inner[1], x_polar_sep.inner[0], x_polar_sep.inner[1]); - ctc_polar.contract(x_cart_sep.outer[0], x_cart_sep.outer[1], x_polar_sep.outer[0], x_polar_sep.outer[1]); + // Complement of the empty polar box: no point is removed by the inner + // contractor, so the Cartesian box is left unchanged. + if(x_pol.is_empty()) + return x; - return x_polar_sep; -} + const Interval& rho = x_pol[0]; + const Interval& theta = x_pol[1]; + std::list parts; -SepPolarXY::SepPolarXY(Interval rho, Interval theta) : rho(rho), theta(theta), Sep(2) { - rho_m = Interval(0, rho.lb()); - rho_p = Interval(rho.ub(), POS_INFINITY); - double limit = theta.mid() - M_PI; - theta_m = Interval(limit, theta.lb()); - theta_p = Interval(theta.ub(), limit + 2*M_PI); - cmpl = Interval(0, 2*M_PI); -} + auto add_polar_part = [&ctc_polar,&parts](const IntervalVector& x, + const Interval& rho, + const Interval& theta) + { + if(rho.is_empty() || theta.is_empty()) + return; + Interval rho_(rho), theta_(theta); + IntervalVector x_(x); + ctc_polar.contract(x_[0],x_[1],rho_,theta_); + parts.push_back(x_); + }; -void SepPolarXY::contractOut(IntervalVector &x_out){ - Interval th = this->theta; - Interval r = this->rho; - this->ctc.contract(x_out[0], x_out[1], r, th); - if(x_out[0].is_empty() || x_out[1].is_empty()) - x_out.set_empty(); -} + // Radial complement, restricted to the current angular sector. The angular + // outside parts below cover the remaining angles. + if(std::isfinite(rho.lb()) && rho.lb() >= 0) + add_polar_part(x, {0,rho.lb()}, theta); + + if(std::isfinite(rho.ub())) + add_polar_part(x, {rho.ub(),oo}, theta); + // Angular complement. If [theta] already covers at least one full turn, + // there is no angular complement in the 2*pi-periodic sense. + if(std::isfinite(theta.lb()) && std::isfinite(theta.ub()) + && theta.ub() - theta.lb() < 2*PI) + { + double limit = theta.mid()-PI; + add_polar_part(x, {0,oo}, {limit,theta.lb()}); + add_polar_part(x, {0,oo}, {theta.ub(),limit+2*PI}); + } -void SepPolarXY::contractIn(IntervalVector &x_in){ - Interval x1(x_in[0]); Interval y1(x_in[1]); - Interval x2(x_in[0]); Interval y2(x_in[1]); - Interval x3(x_in[0]); Interval y3(x_in[1]); - Interval x4(x_in[0]); Interval y4(x_in[1]); - - Interval ALLREALS1 = Interval::POS_REALS; - Interval ALLREALS2 = Interval::POS_REALS; - Interval cmpl1(cmpl); - Interval cmpl2(cmpl); - Interval theta_m_tmp(theta_m); - Interval theta_p_tmp(theta_p); - Interval rho_m_tmp(rho_m); - Interval rho_p_tmp(rho_p); - - this->ctc.contract(x1, y1, ALLREALS1, theta_m_tmp); - this->ctc.contract(x2, y2, ALLREALS2, theta_p_tmp); - this->ctc.contract(x3, y3, rho_m_tmp, cmpl1); - this->ctc.contract(x4, y4, rho_p_tmp, cmpl2); - x_in[0] &= (x1 | x2 | x3 | x4); - x_in[1] &= (y1 | y2 | y3 | y4); - if(x_in[0].is_empty() || x_in[1].is_empty()) - x_in.set_empty(); + return x & hull(parts); } + +BoxPair SepPolarCart::separate(const IntervalVector& x) const +{ + assert_release(x.size() == 2); + + BoxPair x_cart_sep { x, x }; + + if(x.is_empty()) + return x_cart_sep; + + CtcPolar ctc_polar; + IntervalVector x_cart(x), x_pol(2); + + // Cartesian input -> polar enclosure + ctc_polar.contract(x_cart[0], x_cart[1], x_pol[0], x_pol[1]); + + if(x_cart.is_empty() || x_pol.is_empty()) + return x_cart_sep; + + // Separation in the polar space + const BoxPair x_pol_sep = _sep.front()->separate(x_pol); + + // Outer contraction: keep the Cartesian points whose polar coordinates may + // belong to the polar separator output + IntervalVector copy_x_pol_sep = x_pol_sep.outer.subvector(0,1); + ctc_polar.contract(x_cart_sep.outer[0], x_cart_sep.outer[1], copy_x_pol_sep[0], copy_x_pol_sep[1]); + + // Inner contraction: remove points belonging to the polar set by contracting + // with the complement of the polar outer box + x_cart_sep.inner = contract_polar_box_complement(ctc_polar, x_cart_sep.inner, x_pol_sep.outer); + + return x_cart_sep; +} \ No newline at end of file diff --git a/src/core/separators/codac2_SepPolarCart.h b/src/core/separators/codac2_SepPolarCart.h index 27651702f..e4f6b6e9f 100644 --- a/src/core/separators/codac2_SepPolarCart.h +++ b/src/core/separators/codac2_SepPolarCart.h @@ -2,7 +2,7 @@ * \file codac2_SepPolarCart.h * ---------------------------------------------------------------------------- * \date 2026 - * \author Benoît Desrochers, (Simon Rohou) + * \author Benoît Desrochers, Simon Rohou * \copyright Copyright 2026 Codac Team * \license GNU Lesser General Public License (LGPL) */ @@ -14,7 +14,15 @@ namespace codac2 { - // For separating a Cartesian box from a Sep expressed in polar coordinates + /** + * \class SepPolarCart + * \brief For separating a Cartesian box from a separator expressed in polar + * coordinates \f$(\rho,\theta)\f$. + * + * This separator is the inverse counterpart of SepCartPolar: the wrapped + * separator acts on polar boxes, while this class exposes a separator on + * Cartesian boxes \f$(x,y)\f$. + */ class SepPolarCart : public Sep { public: @@ -33,4 +41,4 @@ namespace codac2 const Collection _sep; }; -} \ No newline at end of file +} diff --git a/src/graphics/figures/codac2_Figure2D.cpp b/src/graphics/figures/codac2_Figure2D.cpp index 53045ebc1..15c9503c8 100644 --- a/src/graphics/figures/codac2_Figure2D.cpp +++ b/src/graphics/figures/codac2_Figure2D.cpp @@ -240,55 +240,61 @@ void Figure2D::draw_polygon(const Polygon& x, const StyleProperties& style) void Figure2D::draw_zonotope(const Zonotope& z, const StyleProperties& style) { - map sides; - for (int i=0; i < z.A.cols(); i++) { - auto u = z.A.col(i); - assert_release(u.size()==2); - if (u==Vector::zero(2)) continue; - double theta = std::atan2(u[1],u[0]); - Vector v(u); - if (theta<=0.0) { theta=theta+PI; v=-v; } + if (!z.is_empty()) + { + map sides; + for (int i=0; i < z.A.cols(); i++) { + auto u = z.A.col(i); + assert_release(u.size()==2); + if (u==Vector::zero(2)) continue; + double theta = std::atan2(u[1],u[0]); + Vector v(u); + if (theta<=0.0) { theta=theta+PI; v=-v; } // Theta in ]0,PI] , v[1]>=0 and if v[1]=0, v[0]<0 - auto try_insert=sides.insert({theta,v}); - if (try_insert.second==false) { - (try_insert.first)->second += v; - } - } - vector vertices; - Vector point=z.z; - // Start from v[1] maximum (and v[0] min for horizontal side) - for (const auto& a : sides) { - point+=a.second; - } - // Turn anticlockwise : first half - for (const auto& a : sides) { - vertices.push_back(point); - point-=2*a.second; - } - // Turn anticlockwise : second half - for (const auto& a : sides) { - vertices.push_back(point); - point+=2*a.second; - } - for(const auto& output_fig : _output_figures) + auto try_insert=sides.insert({theta,v}); + if (try_insert.second==false) { + (try_insert.first)->second += v; + } + } + vector vertices; + Vector point=z.c; + // Start from v[1] maximum (and v[0] min for horizontal side) + for (const auto& a : sides) { + point+=a.second; + } + // Turn anticlockwise : first half + for (const auto& a : sides) { + vertices.push_back(point); + point-=2*a.second; + } + // Turn anticlockwise : second half + for (const auto& a : sides) { + vertices.push_back(point); + point+=2*a.second; + } + for(const auto& output_fig : _output_figures) output_fig->draw_polygon(vertices,style); + } } void Figure2D::draw_parallelepiped(const Parallelepiped& p, const StyleProperties& style) { - assert_release(p.A.is_squared() && p.A.rows() == p.z.size()); - assert_release(p.z.size() == 2); - - auto a1 = p.A.col(0), a2 = p.A.col(1); + assert_release(p.A.is_squared() && p.A.rows() == p.c.size()); + assert_release(p.c.size() == 2); - if (a1.isZero() || a2.isZero()) - draw_polyline({p.z-a1-a2,p.z+a1+a2}, style); - else - draw_polygon({ - p.z+a1+a2, p.z-a1+a2, - p.z-a1-a2, p.z+a1-a2 - }, style); + if (!p.is_empty()) + { + auto a1 = p.A.col(0), a2 = p.A.col(1); + + if (a1.isZero() || a2.isZero()) + draw_polyline({p.c-a1-a2,p.c+a1+a2}, style); + else + draw_polygon({ + p.c+a1+a2, p.c-a1+a2, + p.c-a1-a2, p.c+a1-a2 + }, style); + } } void Figure2D::draw_pie(const Vector& c, const Interval& r, const Interval& theta, const StyleProperties& style) diff --git a/src/graphics/figures/codac2_Figure2D.h b/src/graphics/figures/codac2_Figure2D.h index 9a950ea18..3b920c734 100644 --- a/src/graphics/figures/codac2_Figure2D.h +++ b/src/graphics/figures/codac2_Figure2D.h @@ -310,7 +310,7 @@ namespace codac2 void draw_parallelepiped(const Parallelepiped& p, const StyleProperties& style = StyleProperties()); /** - * \brief Draws a zonotope z+sum_i [-1,1] A_i on the figure + * \brief Draws a zonotope c+sum_i [-1,1] A_i on the figure * * \param z Zonotope to draw (center and shape matrix) * \param style Style of the zonotope (edge color and fill color) @@ -856,15 +856,15 @@ namespace codac2 } /** - * \brief Draws a zonotope z+sum_i [-1,1] A_i on the figure + * \brief Draws a zonotope c+sum_i [-1,1] A_i on the figure * - * \param z Zonotope to draw (center and shape matrix) + * \param c Zonotope to draw (center and shape matrix) * \param style Style of the zonotope (edge color and fill color) */ - static void draw_zonotope(const Zonotope& z, const StyleProperties& style = StyleProperties()) + static void draw_zonotope(const Zonotope& c, const StyleProperties& style = StyleProperties()) { auto_init(); - selected_fig()->draw_zonotope(z,style); + selected_fig()->draw_zonotope(c,style); } /** diff --git a/src/graphics/figures/codac2_Figure3D.cpp b/src/graphics/figures/codac2_Figure3D.cpp index 86d7fcdc0..ca0aee3ef 100644 --- a/src/graphics/figures/codac2_Figure3D.cpp +++ b/src/graphics/figures/codac2_Figure3D.cpp @@ -105,19 +105,19 @@ void Figure3D::draw_parallelogram(const Vector &c, const Matrix &A, void Figure3D::draw_parallelepiped(const Parallelepiped& p, const StyleProperties& style) { - assert_release(p.z.size() == 3); + assert_release(p.c.size() == 3); assert_release(p.A.rows() == 3 && p.A.cols() == 3); this->set_style_internal(style); - size_t ip0 = this->move_write_v(p.z,p.A,Vector({-1,-1,-1})); - size_t ip1 = this->move_write_v(p.z,p.A,Vector({-1,-1,1})); - size_t ip2 = this->move_write_v(p.z,p.A,Vector({-1,1,-1})); - size_t ip3 = this->move_write_v(p.z,p.A,Vector({-1,1,1})); - size_t ip4 = this->move_write_v(p.z,p.A,Vector({1,-1,-1})); - size_t ip5 = this->move_write_v(p.z,p.A,Vector({1,-1,1})); - size_t ip6 = this->move_write_v(p.z,p.A,Vector({1,1,-1})); - size_t ip7 = this->move_write_v(p.z,p.A,Vector({1,1,1})); + size_t ip0 = this->move_write_v(p.c,p.A,Vector({-1,-1,-1})); + size_t ip1 = this->move_write_v(p.c,p.A,Vector({-1,-1,1})); + size_t ip2 = this->move_write_v(p.c,p.A,Vector({-1,1,-1})); + size_t ip3 = this->move_write_v(p.c,p.A,Vector({-1,1,1})); + size_t ip4 = this->move_write_v(p.c,p.A,Vector({1,-1,-1})); + size_t ip5 = this->move_write_v(p.c,p.A,Vector({1,-1,1})); + size_t ip6 = this->move_write_v(p.c,p.A,Vector({1,1,-1})); + size_t ip7 = this->move_write_v(p.c,p.A,Vector({1,1,1})); _file << "f " << ip0 << " " << ip1 << " " << ip3 << " " << ip2 << "\n"; _file << "f " << ip4 << " " << ip5 << " " << ip7 << " " << ip6 << "\n"; @@ -137,7 +137,7 @@ void Figure3D::draw_box(const IntervalVector& x, const StyleProperties& style) } void Figure3D::draw_zonotope(const Zonotope& z, const StyleProperties& style) { - assert_release(z.z.size() == 3); + assert_release(z.c.size() == 3); Matrix id = Matrix::Identity(3,3); this->set_style_internal(style); lock_style=true; @@ -186,8 +186,8 @@ void Figure3D::draw_zonotope(const Zonotope& z, const StyleProperties& style) { R1 -= Ak; } } - this->draw_parallelogram(z.z,id,R1+R2,Ai,Aj,style); - this->draw_parallelogram(z.z,id,-R1+R2,Ai,Aj,style); + this->draw_parallelogram(z.c,id,R1+R2,Ai,Aj,style); + this->draw_parallelogram(z.c,id,-R1+R2,Ai,Aj,style); } } lock_style=false; @@ -214,10 +214,11 @@ void Figure3D::draw_arrow(const Vector& c, const Matrix &A, } -void Figure3D::draw_axes(double size) +void Figure3D::draw_axes(double size, const Vector& origin) { + assert_release(origin.size()==3); const std::string name="axes"; - Vector z = Vector::zero(3); + Vector z = origin; // X axis Matrix AX = Matrix::Identity(3,3); draw_arrow(z,size*AX,StyleProperties(Color::red(),name)); diff --git a/src/graphics/figures/codac2_Figure3D.h b/src/graphics/figures/codac2_Figure3D.h index 3887bdee1..442293832 100644 --- a/src/graphics/figures/codac2_Figure3D.h +++ b/src/graphics/figures/codac2_Figure3D.h @@ -116,7 +116,7 @@ namespace codac2 void draw_parallelepiped(const Parallelepiped& p, const StyleProperties& style = { Color::dark_gray(0.5) }); /** - * \brief Draws a zonotope z+sum_i [-1,1] A_i on the figure + * \brief Draws a zonotope c+sum_i [-1,1] A_i on the figure * * \param z Zonotope to draw (center and shape matrix) * \param style Style of the zonotope (edge color) @@ -146,8 +146,9 @@ namespace codac2 * \brief Draws the (x,y,z) axes on the figure in red, green and blue * * \param size Size of the axes + * \param origin Origin of the axes */ - void draw_axes(double size = 1.0); + void draw_axes(double size = 1.0, const Vector& origin = Vector::Zero(3)); /** * \brief Draws a parametric surface diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c38934659..4ea53368d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -54,8 +54,9 @@ list(APPEND SRC_TESTS # listing files without extension core/domains/interval/codac2_tests_IntervalMatrix core/domains/interval/codac2_tests_IntervalVector ../doc/manual/manual/intervals/src - core/domains/zonotope/codac2_tests_Parallelepiped core/domains/zonotope/codac2_tests_Parallelepiped_eval + core/domains/zonotope/codac2_tests_Parallelepiped + core/domains/zonotope/codac2_tests_Zonotope core/domains/tube/codac2_tests_TDomain core/domains/tube/codac2_tests_Slice core/domains/tube/codac2_tests_Slice_polygon diff --git a/tests/core/domains/interval/codac2_tests_Interval.cpp b/tests/core/domains/interval/codac2_tests_Interval.cpp index 7dfe8ba2c..e98d6ca55 100644 --- a/tests/core/domains/interval/codac2_tests_Interval.cpp +++ b/tests/core/domains/interval/codac2_tests_Interval.cpp @@ -281,6 +281,27 @@ TEST_CASE("Interval - other tests") CHECK(Approx(a.second,1e-7) == Interval(2,10)); } +TEST_CASE("Interval - comparison") +{ + CHECK((Interval(0,1) < Interval(2,3)) == BoolInterval::TRUE); + CHECK((Interval(0,1) < Interval(1,2)) == BoolInterval::UNKNOWN); + CHECK((Interval(0,1) < Interval(0,1)) == BoolInterval::UNKNOWN); + CHECK((Interval(2,3) < Interval(0,1)) == BoolInterval::FALSE); + CHECK((Interval(2,3) < Interval(-oo,2)) == BoolInterval::FALSE); + CHECK((Interval(-oo,2) < Interval(3,oo)) == BoolInterval::TRUE); + CHECK((Interval::empty() < Interval(0,1)) == BoolInterval::EMPTY); + CHECK((Interval(0,1) < Interval::empty()) == BoolInterval::EMPTY); + + CHECK((Interval(2,3) > Interval(0,1)) == BoolInterval::TRUE); + CHECK((Interval(1,2) > Interval(0,1)) == BoolInterval::UNKNOWN); + CHECK((Interval(0,1) > Interval(0,1)) == BoolInterval::UNKNOWN); + CHECK((Interval(0,1) > Interval(2,3)) == BoolInterval::FALSE); + CHECK((Interval(-oo,2) > Interval(2,3)) == BoolInterval::FALSE); + CHECK((Interval(3,oo) > Interval(-oo,2)) == BoolInterval::TRUE); + CHECK((Interval::empty() > Interval(0,1)) == BoolInterval::EMPTY); + CHECK((Interval(0,1) > Interval::empty()) == BoolInterval::EMPTY); +} + #if 0 // Tests from the IBEX lib that are not considered in this file: diff --git a/tests/core/domains/interval/codac2_tests_Interval.py b/tests/core/domains/interval/codac2_tests_Interval.py index 858a88932..9a505e469 100644 --- a/tests/core/domains/interval/codac2_tests_Interval.py +++ b/tests/core/domains/interval/codac2_tests_Interval.py @@ -249,6 +249,26 @@ def test_interval(self): b = Interval(1,3) & Interval(6,7) # [b] is empty c = a+b self.assertTrue(c.is_empty()) + + # Comparisons + + self.assertTrue((Interval(0,1) < Interval(2,3)) == BoolInterval.TRUE) + self.assertTrue((Interval(0,1) < Interval(1,2)) == BoolInterval.UNKNOWN) + self.assertTrue((Interval(0,1) < Interval(0,1)) == BoolInterval.UNKNOWN) + self.assertTrue((Interval(2,3) < Interval(0,1)) == BoolInterval.FALSE) + self.assertTrue((Interval(2,3) < Interval(-oo,2)) == BoolInterval.FALSE) + self.assertTrue((Interval(-oo,2) < Interval(3,oo)) == BoolInterval.TRUE) + self.assertTrue((Interval.empty() < Interval(0,1)) == BoolInterval.EMPTY) + self.assertTrue((Interval(0,1) < Interval.empty()) == BoolInterval.EMPTY) + + self.assertTrue((Interval(2,3) > Interval(0,1)) == BoolInterval.TRUE) + self.assertTrue((Interval(1,2) > Interval(0,1)) == BoolInterval.UNKNOWN) + self.assertTrue((Interval(0,1) > Interval(0,1)) == BoolInterval.UNKNOWN) + self.assertTrue((Interval(0,1) > Interval(2,3)) == BoolInterval.FALSE) + self.assertTrue((Interval(-oo,2) > Interval(2,3)) == BoolInterval.FALSE) + self.assertTrue((Interval(3,oo) > Interval(-oo,2)) == BoolInterval.TRUE) + self.assertTrue((Interval.empty() > Interval(0,1)) == BoolInterval.EMPTY) + self.assertTrue((Interval(0,1) > Interval.empty()) == BoolInterval.EMPTY) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp index 2c5766d4b..68f9b5950 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp @@ -8,6 +8,7 @@ */ #include +#include #include using namespace std; @@ -26,8 +27,42 @@ TEST_CASE("Parallelepiped") CHECK((p.is_superset(IntervalVector({{0.,5.},{2.,7.},{4.,9.}})))==BoolInterval::UNKNOWN); Zonotope z = p.proj({2,1,0}); - CHECK(z.z == Vector({4,2,0})); + CHECK(z.c == Vector({4,2,0})); CHECK(z.A == Matrix({{0.,1.,1.},{0.,1.,0.},{0.5,0.,0.}})); CHECK(z.box() == IntervalVector({{2.,6.},{1.,3.},{-0.5,0.5}})); + + Parallelepiped p1 ({1.5,2.8},Matrix({{0.5,0.4},{0,0.2}})); + + IntervalVector x0 ({{0,5},{0,5}}); + IntervalVector x1 ({{0,0.5},{0,0.5}}); + IntervalVector x2 ({{1.7,5},{0,3}}); + + CHECK(Approx(x0 & p1) == IntervalVector({{0.6,2.4},{2.6,3}})); + CHECK(Approx((p1 & x0).c) == p1.c); + CHECK(Approx((p1 & x0).A) == p1.A); + + CHECK((x1 & p1).is_empty()); + CHECK((p1 & x1).is_empty()); + + CHECK(Approx(x2 & p1) == IntervalVector({{1.7,2.4},{2.6,3}})); + CHECK(Approx((p1 & x2).c) == Vector({1.65, 2.8})); + CHECK(Approx((p1 & x2).A) == Matrix({{0.35,0.4},{0,0.2}})); + + Parallelepiped p2 ({-5,-4,3,8},Matrix({{0.5,0.4,0.1,0},{0,0.2,0.3,0},{0,0,0.1,0.4},{0,0,0,0.2}})); + Parallelepiped p_empty = Parallelepiped::empty(4); + + IntervalVector x3 ({{-10,10},{-10,10},{-10,10},{-10,10}}); + IntervalVector x4 ({{-0.5,0.5},{-0.5,0.5},{-0.5,0.5},{-0.5,0.5}}); + IntervalVector x_empty = IntervalVector::empty(4); + + CHECK(Approx(x3 & p2) == IntervalVector({{-6,-4},{-4.5,-3.5},{2.5,3.5},{7.8,8.2}})); + CHECK(Approx((p2 & x3).c) == Vector({-5,-4,3,8})); + CHECK(Approx((p2 & x3).A) == Matrix({{0.5,0.4,0.1,0},{0,0.2,0.3,0},{0,0,0.1,0.4},{0,0,0,0.2}})); + + CHECK((x4 & p2).is_empty()); + CHECK((p2 & x4).is_empty()); + + CHECK((x_empty & p2).is_empty()); + CHECK((p_empty & x3).is_empty()); } diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py index 75b4cd987..bc201cc0d 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py @@ -27,9 +27,44 @@ def test_parallelepiped(self): self.assertTrue((p.is_superset(IntervalVector([[0.,5.],[2.,7.],[4.,9.]])))==BoolInterval.UNKNOWN) z = p.proj([2,1,0]) - self.assertTrue(z.z == Vector([4,2,0])) + self.assertTrue(z.c == Vector([4,2,0])) self.assertTrue(z.A == Matrix([[0,1,1],[0,1,0],[0.5,0,0]])) self.assertTrue(z.box() == IntervalVector([[2,6],[1,3],[-0.5,0.5]])) + + + p1 = Parallelepiped([1.5,2.8],Matrix([[0.5,0.4],[0,0.2]])) + + x0 = IntervalVector([[0,5],[0,5]]) + x1 = IntervalVector([[0,0.5],[0,0.5]]) + x2 = IntervalVector([[1.7,5],[0,3]]) + + self.assertTrue(Approx(x0 & p1) == IntervalVector([[0.6,2.4],[2.6,3]])) + self.assertTrue(Approx((p1 & x0).c) == p1.c) + self.assertTrue(Approx((p1 & x0).A) == p1.A) + + self.assertTrue((x1 & p1).is_empty()) + self.assertTrue((p1 & x1).is_empty()) + + self.assertTrue(Approx(x2 & p1) == IntervalVector([[1.7,2.4],[2.6,3]])) + self.assertTrue(Approx((p1 & x2).c) == Vector([1.65, 2.8])) + self.assertTrue(Approx((p1 & x2).A) == Matrix([[0.35,0.4],[0,0.2]])) + + p2 = Parallelepiped([-5,-4,3,8],Matrix([[0.5,0.4,0.1,0],[0,0.2,0.3,0],[0,0,0.1,0.4],[0,0,0,0.2]])) + p_empty = Parallelepiped.empty(4) + + x3 = IntervalVector([[-10,10],[-10,10],[-10,10],[-10,10]]) + x4 = IntervalVector([[-0.5,0.5],[-0.5,0.5],[-0.5,0.5],[-0.5,0.5]]) + x_empty = IntervalVector.empty(4) + + self.assertTrue(Approx(x3 & p2) == IntervalVector([[-6,-4],[-4.5,-3.5],[2.5,3.5],[7.8,8.2]])) + self.assertTrue(Approx((p2 & x3).c) == Vector([-5,-4,3,8])) + self.assertTrue(Approx((p2 & x3).A) == Matrix([[0.5,0.4,0.1,0],[0,0.2,0.3,0],[0,0,0.1,0.4],[0,0,0,0.2]])) + + self.assertTrue((x4 & p2).is_empty()) + self.assertTrue((p2 & x4).is_empty()) + + self.assertTrue((x_empty & p2).is_empty()) + self.assertTrue((p_empty & x3).is_empty()) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp index 0ad038b38..9dd1705b2 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp @@ -28,17 +28,17 @@ TEST_CASE("Parallelepiped_eval") auto p1a = f1.parallelepiped_eval(Interval(-0.1,0.1)); auto p1b = f1.parallelepiped_eval(1.0); - CHECK(Approx(p1a.z,1e-6) == Vector({0,0})); + CHECK(Approx(p1a.c,1e-6) == Vector({0,0})); CHECK(Approx(p1a.A,1e-6) == Matrix({{0.12,0},{0,0.02}})); - CHECK(Approx(p1b.z,1e-6) == Vector({1,1})); + CHECK(Approx(p1b.c,1e-6) == Vector({1,1})); CHECK(Approx(p1b.A,1e-6) == Matrix({{0,0},{0,0}})); auto pa = f2.parallelepiped_eval(Interval(-0.1,0.1), Interval(-0.1,0.1)); auto pb = f2.parallelepiped_eval(1.0, Interval(-1,1)); - CHECK(Approx(pa.z,1e-6) == Vector({0,0,0})); + CHECK(Approx(pa.c,1e-6) == Vector({0,0,0})); CHECK(Approx(pa.A,1e-6) == Matrix({{0.14,0,0},{0,0.14,0},{0,0,0.04}})); - CHECK(Approx(pb.z,1e-6) == Vector({1,0,1})); + CHECK(Approx(pb.c,1e-6) == Vector({1,0,1})); CHECK(Approx(pb.A,1e-5) == Matrix({{0.894428,0,1.78886},{0,3,0},{1.78886,0,-0.894427}})); @@ -55,7 +55,7 @@ TEST_CASE("Parallelepiped_eval") auto p2 = f2.parallelepiped_eval(X0,Y0); auto p3 = f3.parallelepiped_eval(IntervalVector({X0,Y0})); - CHECK(Approx(p2.z,1e-6) == p3.z); + CHECK(Approx(p2.c,1e-6) == p3.c); CHECK(Approx(p2.A,1e-6) == p3.A); y0+=dx; } diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py index c7dcd18ed..5cac5ffe7 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py @@ -27,17 +27,17 @@ def test_parallelepiped_eval(self): p1a = f1.parallelepiped_eval(Interval(-0.1,0.1)) p1b = f1.parallelepiped_eval(1.0) - self.assertTrue(Approx(p1a.z,1e-6)==Vector([0.0,0.0])) + self.assertTrue(Approx(p1a.c,1e-6)==Vector([0.0,0.0])) self.assertTrue(Approx(p1a.A,1e-6)==Matrix([[0.12,0.0],[0.0,0.02]])) - self.assertTrue(Approx(p1b.z,1e-6)==Vector([1.0,1.0])) + self.assertTrue(Approx(p1b.c,1e-6)==Vector([1.0,1.0])) self.assertTrue(Approx(p1b.A,1e-6)==Matrix([[0.0,0.0],[0.0,0.0]])) pa = f2.parallelepiped_eval(Interval(-0.1,0.1), Interval(-0.1,0.1)) pb = f2.parallelepiped_eval(1.0,Interval(-1,1)) - self.assertTrue(Approx(pa.z,1e-6)==Vector([0,0,0])) + self.assertTrue(Approx(pa.c,1e-6)==Vector([0,0,0])) self.assertTrue(Approx(pa.A,1e-6)==Matrix([[0.14,0,0],[0,0.14,0],[0,0,0.04]])) - self.assertTrue(Approx(pb.z,1e-6)==Vector([1,0,1])) + self.assertTrue(Approx(pb.c,1e-6)==Vector([1,0,1])) self.assertTrue(Approx(pb.A,1e-5)==Matrix([[0.894428,0,1.78886],[0,3,0],[1.78886,0,-0.894427]])) @@ -52,7 +52,7 @@ def test_parallelepiped_eval(self): p2 = f2.parallelepiped_eval(X0,Y0) p3 = f3.parallelepiped_eval(IntervalVector([X0,Y0])) - self.assertTrue(Approx(p2.z,1e-6)==p3.z) + self.assertTrue(Approx(p2.c,1e-6)==p3.c) self.assertTrue(Approx(p2.A,1e-6)==p3.A) y0 += dx x0 += dx diff --git a/tests/core/domains/zonotope/codac2_tests_Zonotope.cpp b/tests/core/domains/zonotope/codac2_tests_Zonotope.cpp new file mode 100644 index 000000000..60da5d7ea --- /dev/null +++ b/tests/core/domains/zonotope/codac2_tests_Zonotope.cpp @@ -0,0 +1,28 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("Zonotope") +{ + Zonotope Z1 (Vector({2,1}), Matrix({{0.2 ,0.08}, {0.04,0.18}})); + Zonotope Z2 (Vector({2,0.5}), Matrix({{-0.2}, {0.1}})); + auto Zs = Z1+Z2; + + Vector c ({4,1.5}); + Matrix A ({{0.2, 0.08, -0.2}, + {0.04, 0.18, 0.1}}); + CHECK(Zs.c == c); + CHECK(Zs.A == A); +} + diff --git a/tests/core/domains/zonotope/codac2_tests_Zonotope.py b/tests/core/domains/zonotope/codac2_tests_Zonotope.py new file mode 100644 index 000000000..46ddd969a --- /dev/null +++ b/tests/core/domains/zonotope/codac2_tests_Zonotope.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2025 +# \author Maël Godard +# \copyright Copyright 2024 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +from codac import * +import sys +import math + +class TestZonotope(unittest.TestCase): + + def test_zonotope(self): + + Z1 = Zonotope(Vector([2,1]), Matrix([[0.2 ,0.08], [0.04,0.18]])) + Z2 = Zonotope(Vector([2,0.5]), Matrix([[-0.2], [0.1]])) + Zs = Z1+Z2 + + c = Vector([4,1.5]) + A = Matrix([[0.2 ,0.08, -0.2], + [0.04, 0.18, 0.1]]) + self.assertTrue(Zs.c == c) + self.assertTrue(Zs.A == A) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/core/matrices/codac2_tests_IntFullPivLU.cpp b/tests/core/matrices/codac2_tests_IntFullPivLU.cpp index 2e05e6fce..4fe00f00b 100644 --- a/tests/core/matrices/codac2_tests_IntFullPivLU.cpp +++ b/tests/core/matrices/codac2_tests_IntFullPivLU.cpp @@ -146,3 +146,68 @@ TEST_CASE("IntvFullPivLU") } } + +TEST_CASE("IntvFullPivLU solve tall matrix") +{ + Matrix A({ + { 1, 0 }, + { 0, 1 }, + { 1, 1 }, + { 2,-1 }, + {-1, 2 }, + }); + + Matrix X({ + { 2,-1 }, + {-1, 3 }, + }); + + Matrix rhs_mid = A*X; + IntervalMatrix rhs(rhs_mid); + + IntvFullPivLU lu(A); + + SECTION("solve(rhs) with rows >= cols + 2") + { + IntervalMatrix sol = lu.solve(rhs); + + CHECK(!sol.is_empty()); + CHECK(sol.contains(X)); + CHECK((A.template cast() * sol).contains(rhs_mid)); + } + + SECTION("solve(rhs, B) with rows >= cols + 2") + { + IntervalMatrix B = IntervalMatrix::Constant(2,2,{-10.0,10.0}); + + lu.solve(rhs, B); + + CHECK(!B.is_empty()); + CHECK(B.contains(X)); + CHECK((A.template cast() * B).contains(rhs_mid)); + } + + SECTION("solve(rhs) detects inconsistent rhs") + { + Matrix rhs_bad_mid = rhs_mid; + rhs_bad_mid(4,0) += 1.0; // breaks consistency on the extra row + + IntervalMatrix rhs_bad(rhs_bad_mid); + IntervalMatrix sol = lu.solve(rhs_bad); + + CHECK(sol.is_empty()); + } + + SECTION("solve(rhs, B) detects inconsistent rhs") + { + Matrix rhs_bad_mid = rhs_mid; + rhs_bad_mid(4,0) += 1.0; // breaks consistency on the extra row + + IntervalMatrix rhs_bad(rhs_bad_mid); + IntervalMatrix B = IntervalMatrix::Constant(2,2,{-10.0,10.0}); + + lu.solve(rhs_bad, B); + + CHECK(B.is_empty()); + } +} \ No newline at end of file diff --git a/tests/core/matrices/codac2_tests_IntFullPivLU.py b/tests/core/matrices/codac2_tests_IntFullPivLU.py index e658f3967..5d8209279 100644 --- a/tests/core/matrices/codac2_tests_IntFullPivLU.py +++ b/tests/core/matrices/codac2_tests_IntFullPivLU.py @@ -13,72 +13,150 @@ class TestIntvFullPivLU(unittest.TestCase): + def test_IntvFullPivLU_1(self): + M = Matrix([ + [ 1, -4, 3, 7 ], + [ 2, 1, -4, 6 ], + [ 5, 2, 1 , 9 ], + [ -1, 0, 3, 2 ] + ]) -# bindings not implemented - def test_IntvFullPivLU_1(self): - M = Matrix([ - [ 1, -4, 3, 7 ], - [ 2, 1, -4, 6 ], - [ 5, 2, 1 , 9 ], - [ -1, 0, 3, 2 ] - ]) + LUdec = IntvFullPivLU(M) + self.assertTrue(LUdec.is_injective()==BoolInterval.TRUE) + self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE) + self.assertTrue(LUdec.rank()==Interval(4)) + self.assertTrue((LUdec.determinant()+602).mag()<=1e-10) + self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10) + I1 = LUdec.solve(IntervalMatrix.eye(4,4)); + self.assertTrue((I1*M-Matrix.eye(4,4)).norm().ub()<1e-10) - LUdec = IntvFullPivLU(M) - self.assertTrue(LUdec.is_injective()==BoolInterval.TRUE) - self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE) - self.assertTrue(LUdec.rank()==Interval(4)) - self.assertTrue((LUdec.determinant()+602).mag()<=1e-10) - self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10) - I1 = LUdec.solve(IntervalMatrix.eye(4,4)); - self.assertTrue((I1*M-Matrix.eye(4,4)).norm().ub()<1e-10) + def test_IntvFullPivLU_2(self): + M = Matrix([ + [ 1, -4, 6, 7 ], + [ 2, 1, 3, 6 ], + [ 5, 2, 8 , 9 ], + [ -1, 0, -2, 2 ] + ]) - def test_IntvFullPivLU_2(self): - M = Matrix([ - [ 1, -4, 6, 7 ], - [ 2, 1, 3, 6 ], - [ 5, 2, 8 , 9 ], - [ -1, 0, -2, 2 ] - ]) + LUdec = IntvFullPivLU(M) + self.assertTrue(LUdec.is_injective()==BoolInterval.UNKNOWN) + self.assertTrue(LUdec.is_surjective()==BoolInterval.UNKNOWN) + self.assertTrue(LUdec.rank()==Interval([3,4])) + self.assertTrue((LUdec.determinant()).mag()<=1e-10) + self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10) + K = LUdec.kernel() + self.assertTrue(K.cols()==1) + self.assertTrue((M*K).norm().ub()<1e-10) + coK = LUdec.cokernel() + self.assertTrue(coK.rows()==1) + self.assertTrue((coK*M).norm().ub()<1e-10) + Im = LUdec.image(M) + self.assertTrue(Im.cols()==3) + coIm = LUdec.coimage(M) + self.assertTrue(coIm.rows()==3) - LUdec = IntvFullPivLU(M) - self.assertTrue(LUdec.is_injective()==BoolInterval.UNKNOWN) - self.assertTrue(LUdec.is_surjective()==BoolInterval.UNKNOWN) - self.assertTrue(LUdec.rank()==Interval([3,4])) - self.assertTrue((LUdec.determinant()).mag()<=1e-10) - self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10) - K = LUdec.kernel() - self.assertTrue(K.cols()==1) - self.assertTrue((M*K).norm().ub()<1e-10) - coK = LUdec.cokernel() - self.assertTrue(coK.rows()==1) - self.assertTrue((coK*M).norm().ub()<1e-10) - Im = LUdec.image(M) - self.assertTrue(Im.cols()==3) - coIm = LUdec.coimage(M) - self.assertTrue(coIm.rows()==3) + def test_IntvFullPivLU_3(self): + M = Matrix([ + [ 1, -4, 6, 7, 6 ], + [ 2, 1, 3, 6, -2 ], + [ 5, 2, 8 , 9, -1 ] + ]) + LUdec = IntvFullPivLU(M) + self.assertTrue(LUdec.is_injective()==BoolInterval.FALSE) + self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE) + self.assertTrue(LUdec.rank()==Interval(3)) + self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10) + K = LUdec.kernel() + self.assertTrue(K.cols()==2) + self.assertTrue((M*K).norm().ub()<1e-10) + Im = LUdec.image(M) + self.assertTrue(Im.cols()==3) + I1 = LUdec.solve(IntervalMatrix.eye(3,3)); + self.assertTrue((M*I1-Matrix.eye(3,3)).norm().ub()<1e-10) + A = IntervalMatrix([ [1], [[-20,20]], [2], [[-20,20]], [[-20,20]] ]) + B = IntervalMatrix([ [2.0], [1.0], [4.0] ]) + LUdec.solve(B,A) + self.assertTrue((M*A-B).norm().ub()<=1e-10) - def test_IntvFullPivLU_3(self): - M = Matrix([ - [ 1, -4, 6, 7, 6 ], - [ 2, 1, 3, 6, -2 ], - [ 5, 2, 8 , 9, -1 ] - ]) - LUdec = IntvFullPivLU(M) - self.assertTrue(LUdec.is_injective()==BoolInterval.FALSE) - self.assertTrue(LUdec.is_surjective()==BoolInterval.TRUE) - self.assertTrue(LUdec.rank()==Interval(3)) - self.assertTrue((LUdec.reconstructed_matrix()-M).norm().ub()<=1e-10) - K = LUdec.kernel() - self.assertTrue(K.cols()==2) - self.assertTrue((M*K).norm().ub()<1e-10) - Im = LUdec.image(M) - self.assertTrue(Im.cols()==3) - I1 = LUdec.solve(IntervalMatrix.eye(3,3)); - self.assertTrue((M*I1-Matrix.eye(3,3)).norm().ub()<1e-10) - A = IntervalMatrix([ [1], [[-20,20]], [2], [[-20,20]], [[-20,20]] ]) - B = IntervalMatrix([ [2.0], [1.0], [4.0] ]) - LUdec.solve(B,A) - self.assertTrue((M*A-B).norm().ub()<=1e-10) + def test_IntvFullPivLU_solve_tall_matrix(self): + A = Matrix([ + [ 1, 0 ], + [ 0, 1 ], + [ 1, 1 ], + [ 2,-1 ], + [-1, 2 ], + ]) + + X = Matrix([ + [ 2,-1 ], + [-1, 3 ], + ]) + + rhs_mid = A*X + rhs = IntervalMatrix(rhs_mid) + + lu = IntvFullPivLU(A) + + sol = lu.solve(rhs) + + self.assertTrue(not sol.is_empty()) + self.assertTrue(sol.contains(X)) + self.assertTrue((A*sol).contains(rhs_mid)) + + def test_IntvFullPivLU_solve_tall_matrix_with_box(self): + + A = Matrix([ + [ 1, 0 ], + [ 0, 1 ], + [ 1, 1 ], + [ 2,-1 ], + [-1, 2 ], + ]) + + X = Matrix([ + [ 2,-1 ], + [-1, 3 ], + ]) + + rhs_mid = A*X + rhs = IntervalMatrix(rhs_mid) + + lu = IntvFullPivLU(A) + B = IntervalMatrix.constant(2,2,[-10,10]) + + lu.solve(rhs, B) + + self.assertTrue(not B.is_empty()) + self.assertTrue(B.contains(X)) + self.assertTrue((A*B).contains(rhs_mid)) + + def test_IntvFullPivLU_solve_tall_matrix_inconsistent_rhs(self): + + A = Matrix([ + [ 1, 0 ], + [ 0, 1 ], + [ 1, 1 ], + [ 2,-1 ], + [-1, 2 ], + ]) + + rhs_bad = IntervalMatrix([ + [[ 2, 2],[-1,-1]], + [[-1,-1],[ 3, 3]], + [[ 1, 1],[ 2, 2]], + [[ 5, 5],[-5,-5]], + [[-3,-3],[ 7, 7]], # should be [-4, 7] for consistency + ]) + + lu = IntvFullPivLU(A) + + sol = lu.solve(rhs_bad) + self.assertTrue(sol.is_empty()) + + B = IntervalMatrix.constant(2,2,[-10,10]) + lu.solve(rhs_bad, B) + self.assertTrue(B.is_empty()) + if __name__ == '__main__': - unittest.main() + unittest.main() \ No newline at end of file diff --git a/tests/core/peibos/codac2_tests_peibos.cpp b/tests/core/peibos/codac2_tests_peibos.cpp index bda717048..755b741e7 100644 --- a/tests/core/peibos/codac2_tests_peibos.cpp +++ b/tests/core/peibos/codac2_tests_peibos.cpp @@ -69,12 +69,12 @@ TEST_CASE("Peibos") CHECK(v_par_3d.size() == 6); - CHECK(Approx(v_par_3d[0].z,1e-6) == Vector({1.,0.,0.})); - CHECK(Approx(v_par_3d[1].z,1e-6) == Vector({0.,1.,0.})); - CHECK(Approx(v_par_3d[2].z,1e-6) == Vector({-1.,0.,0.})); - CHECK(Approx(v_par_3d[3].z,1e-6) == Vector({0.,-1.,0.})); - CHECK(Approx(v_par_3d[4].z,1e-6) == Vector({0.,0.,-1.})); - CHECK(Approx(v_par_3d[5].z,1e-6) == Vector({0.,0.,1.})); + CHECK(Approx(v_par_3d[0].c,1e-6) == Vector({1.,0.,0.})); + CHECK(Approx(v_par_3d[1].c,1e-6) == Vector({0.,1.,0.})); + CHECK(Approx(v_par_3d[2].c,1e-6) == Vector({-1.,0.,0.})); + CHECK(Approx(v_par_3d[3].c,1e-6) == Vector({0.,-1.,0.})); + CHECK(Approx(v_par_3d[4].c,1e-6) == Vector({0.,0.,-1.})); + CHECK(Approx(v_par_3d[5].c,1e-6) == Vector({0.,0.,1.})); double a = 4.35066; diff --git a/tests/core/peibos/codac2_tests_peibos.py b/tests/core/peibos/codac2_tests_peibos.py index 66042d8e0..f4f12191e 100644 --- a/tests/core/peibos/codac2_tests_peibos.py +++ b/tests/core/peibos/codac2_tests_peibos.py @@ -63,12 +63,12 @@ def test_peibos(self): self.assertTrue(len(v_par_3d) == 6) - self.assertTrue(Approx(v_par_3d[0].z,1e-6) == Vector([1.,0.,0.])) - self.assertTrue(Approx(v_par_3d[1].z,1e-6) == Vector([0.,1.,0.])) - self.assertTrue(Approx(v_par_3d[2].z,1e-6) == Vector([-1,0.,0.])) - self.assertTrue(Approx(v_par_3d[3].z,1e-6) == Vector([0.,-1.,0.])) - self.assertTrue(Approx(v_par_3d[4].z,1e-6) == Vector([0.,0.,-1.])) - self.assertTrue(Approx(v_par_3d[5].z,1e-6) == Vector([0.,0.,1.])) + self.assertTrue(Approx(v_par_3d[0].c,1e-6) == Vector([1.,0.,0.])) + self.assertTrue(Approx(v_par_3d[1].c,1e-6) == Vector([0.,1.,0.])) + self.assertTrue(Approx(v_par_3d[2].c,1e-6) == Vector([-1,0.,0.])) + self.assertTrue(Approx(v_par_3d[3].c,1e-6) == Vector([0.,-1.,0.])) + self.assertTrue(Approx(v_par_3d[4].c,1e-6) == Vector([0.,0.,-1.])) + self.assertTrue(Approx(v_par_3d[5].c,1e-6) == Vector([0.,0.,1.])) a = 4.35066 diff --git a/tests/test_codac/CMakeLists.txt b/tests/test_codac/CMakeLists.txt index 6303b8f72..0f4479bc9 100644 --- a/tests/test_codac/CMakeLists.txt +++ b/tests/test_codac/CMakeLists.txt @@ -31,7 +31,12 @@ # Compilation + if(FAST_RELEASE) + add_compile_definitions(FAST_RELEASE) + message(STATUS "You are running Codac in fast release mode. (option -DCMAKE_BUILD_TYPE=Release is required)") + endif() + add_executable(${PROJECT_NAME} main.cpp) target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS}) target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS}) - target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES} Ibex::ibex) \ No newline at end of file + target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES} Ibex::ibex) diff --git a/tests/test_codac/main.cpp b/tests/test_codac/main.cpp index de7865754..4004eaf16 100644 --- a/tests/test_codac/main.cpp +++ b/tests/test_codac/main.cpp @@ -1,18 +1,23 @@ + +// Download and run https://github.com/ENSTABretagneRobotics/VIBES/releases/latest +// before to be able to see the figure in VIBes window... + #include +using namespace std; using namespace codac2; -// Download and run https://github.com/ENSTABretagneRobotics/VIBES/releases -// before to be able to see the tube in VIBes window... - int main() { VectorVar x(3); - AnalyticFunction f({x}, vec( - -sqr(x[2])+2*x[2]*sin(x[2]*x[0])+cos(x[2]*x[1]), - 2*x[2]*cos(x[2]*x[0])-sin(x[2]*x[1]) - )); + AnalyticFunction f { {x}, + { + -sqr(x[2])+2*x[2]*sin(x[2]*x[0])+cos(x[2]*x[1]), + 2*x[2]*cos(x[2]*x[0])-sin(x[2]*x[1]) + } + }; - CtcInverse ctc(f, {0.,0.}); - DefaultFigure::pave(IntervalVector({{0,2},{2,4},{0,10}}), ctc, 0.01); -} \ No newline at end of file + CtcInverse ctc(f, {0,0}); + DefaultFigure::pave({{0,2},{2,4},{0,10}}, ctc, 0.004); + return 0; +}