Skip to content

Commit ef09eef

Browse files
SimonRohouSimonRohou
committed
Merge branch 'codac2_dev' of ssh://github.com/codac-team/codac into codac2_dev
2 parents 1ff3b92 + e3a7a51 commit ef09eef

16 files changed

Lines changed: 911 additions & 0 deletions

File tree

‎doc/manual/index.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -243,6 +243,8 @@ User manual
243243
* Analytic contractors
244244
* :ref:`sec-ctc-analytic-ctcinverse`
245245
* CtcInverseNotIn
246+
* Dynamic contractors
247+
* :ref:`sec-ctc-dynamic-ctclohner`
246248
* Geometric contractors
247249
* :ref:`sec-ctc-geom-ctcdist`
248250
* :ref:`sec-ctc-geom-ctcpolar`

‎doc/manual/manual/contractors/dynamic/clohner.png‎

36.7 KB
Loading

‎doc/manual/manual/contractors/dynamic/ctclohner.rst‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
.. _sec-ctc-dynamic-ctclohner:
2+
3+
The CtcLohner contractor
4+
=========================
5+
6+
Main authors: Auguste Bourgois, `Maël Godard <https://godardma.github.io>`_ (porting to Codac2)
7+
8+
Definition
9+
----------
10+
11+
The Lohner contractor uses Lohner's guaranteed integration algorithm to contract a tube according to a differential
12+
equation. This algorithm performs two main steps:
13+
14+
- an estimation of a global enclosure of the system's trajectories over an integration step (*i.e.* the domain of a slice),
15+
which corresponds to the codomain of that slice;
16+
- using this estimation and the input gate (or output gate in the backward mode), the corresponding output gate (or input gate
17+
in the backward mode) is then estimated;
18+
- over a single slice, these steps can be iterated a few times to obtain tighter enclosures for both the gate and the codomain of the slice.
19+
20+
This contractor is supposed to yield better results than the Picard contractor, as long as the tubes are "thin enough".
21+
22+
.. important::
23+
24+
.. math::
25+
26+
\left.\begin{array}{r}\dot{\mathbf{x}}(t)=\mathbf{f}\big(\mathbf{x}(t)\big)\end{array}\right. \longrightarrow \mathcal{C}_{\textrm{Lohner}}\big([\mathbf{x}](\cdot)\big)
27+
28+
.. tabs::
29+
30+
.. code-tab:: py
31+
32+
var = VectorVar(n) # n is the dimension of the system
33+
ctc_lohner = CtcLohner(AnalyticFunction([var], [expr1, expr2, ..., exprn]))
34+
ctc_lohner.contract(x)
35+
36+
.. code-tab:: c++
37+
38+
VectorVar var(n); // n is the dimension of the system
39+
CtcLohner ctc_lohner(AnalyticFunction({x}, {expr1, expr2, ..., exprn}););
40+
ctc_lohner.contract(x);
41+
42+
.. important::
43+
44+
The contractor might throw a runtime error when it cannot find a global enclosure over a specific time step. This usually
45+
happens when the time step is too large, and may therefore be avoided by reducing the discretisation frequency of the tube.
46+
47+
Note that this behaviour may be tackled using automatic step adjustment. However, such a feature is not yet implemented.
48+
49+
50+
Simple example of use
51+
---------------------
52+
53+
In the following example, we consider the two-dimensional system described by:
54+
55+
.. math::
56+
57+
\left(\begin{array}{c}\dot{x}_1\\\dot{x}_2\end{array}\right) = \left(\begin{array}{c}-x_1\\-\sin({x}_2)\end{array}\right).
58+
59+
The Lohner contractor is used for obtaining a two-dimensional tube, considering an initial condition :math:`\mathbf{x}_0\in[0.9,1.1]\times[0.9,1.1]`. The contraction is tested on two tubes :math:`[\mathbf{a}](\cdot)` and :math:`[\mathbf{b}](\cdot)` of different resolutions (width of slices).
60+
61+
.. tabs::
62+
63+
.. code-tab:: py
64+
65+
x0 = IntervalVector([Interval(1,1),Interval(1,1)]) # the box [1,1]×[1,1]..
66+
x0.inflate(0.1) # ..becomes [0.9,1.1]×[0.9,1.1]
67+
68+
# 2d tubes:
69+
ta = create_tdomain([0.,10.0],0.2) # low resolution
70+
a = SlicedTube_IntervalVector(ta, IntervalVector(2))
71+
a.set(x0,0.0) # setting initial value
72+
73+
tb =create_tdomain([0.,10.0],0.01) # high resolution
74+
b = SlicedTube_IntervalVector(tb, IntervalVector(2))
75+
b.set(x0,0.0) # setting initial value
76+
77+
# Defining Lohner contractor from f
78+
x = VectorVar(2)
79+
f = AnalyticFunction([x], [-x[0],-sin(x[1])])
80+
ctc_lohner = CtcLohner(f)
81+
82+
# Contracting the tubes
83+
ctc_lohner.contract(a)
84+
ctc_lohner.contract(b)
85+
86+
# Graphics
87+
fig1 = Figure2D("Lohner_1",GraphicOutput.VIBES|GraphicOutput.IPE)
88+
fig2 = Figure2D("Lohner_2",GraphicOutput.VIBES|GraphicOutput.IPE)
89+
90+
fig1.plot_tube(a[0])
91+
fig1.plot_tube(b[0],StyleProperties([Color.blue(),Color.blue()]))
92+
93+
fig2.plot_tube(a[1])
94+
fig2.plot_tube(b[1],StyleProperties([Color.blue(),Color.blue()]))
95+
96+
.. code-tab:: c++
97+
98+
IntervalVector x0({Interval(1,1),Interval(1,1)}); // the box [1,1]×[1,1]..
99+
x0.inflate(0.1); // ..becomes [0.9,1.1]×[0.9,1.1]
100+
101+
// 2d tubes:
102+
auto ta =create_tdomain({0.,10.0},0.2); // low resolution
103+
SlicedTube a(ta, IntervalVector(2));
104+
a.set(x0,0.0); // setting initial value
105+
106+
auto tb =create_tdomain({0.,10.0},0.01); // high resolution
107+
SlicedTube b(tb, IntervalVector(2));
108+
b.set(x0,0.0); // setting initial value
109+
110+
// Defining Lohner contractor from f
111+
VectorVar x(2);
112+
AnalyticFunction f({x}, {-x[0],-sin(x[1])});
113+
CtcLohner ctc_lohner(f);
114+
115+
// Contracting the tubes
116+
ctc_lohner.contract(a);
117+
ctc_lohner.contract(b);
118+
119+
// Graphics
120+
Figure2D fig1 ("Lohner_1",GraphicOutput::VIBES|GraphicOutput::IPE);
121+
Figure2D fig2 ("Lohner_2",GraphicOutput::VIBES|GraphicOutput::IPE);
122+
123+
fig1.plot_tube(a[0]);
124+
fig1.plot_tube(b[0],{Color::blue(),Color::blue()});
125+
126+
fig2.plot_tube(a[1]);
127+
fig2.plot_tube(b[1],{Color::blue(),Color::blue()});
128+
129+
fig1.set_window_properties({100,50},{1000,400});
130+
fig2.set_window_properties({100,500},{1000,400});
131+
132+
133+
The above code yields the following result for :math:`[a_1](\cdot)` (in gray, large slices) and :math:`[b_1](\cdot)` (in blue, thin slices):
134+
135+
.. figure:: clohner.png
136+
137+
138+
Related content
139+
---------------
140+
141+
.. admonition:: Mathematical documentation
142+
143+
See `Auguste Bourgois' thesis <https://www.ensta-bretagne.fr/jaulin/thesis_auguste.pdf>`_. In particular, check Chapter 4.
144+
145+
.. admonition:: Technical documentation
146+
147+
See the `C++ API documentation of this class <../../../extra_html/api/classcodac2_1_1_ctc_lohner.html>`_.
148+

‎doc/manual/manual/contractors/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ Contractors, separators
66

77
CtcInter <set/ctcinter>
88
CtcInverse <analytic/ctcinverse>
9+
CtcLohner <dynamic/ctclohner>
910
CtcDist <geometric/ctcdist>
1011
CtcPolar <geometric/ctcpolar>
1112

‎examples/14_lohner/CMakeLists.txt‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
# ==================================================================
2+
# codac / basics example - cmake configuration file
3+
# ==================================================================
4+
5+
cmake_minimum_required(VERSION 3.5)
6+
project(codac_example LANGUAGES CXX)
7+
8+
set(CMAKE_CXX_STANDARD 20)
9+
set(CMAKE_CXX_STANDARD_REQUIRED ON)
10+
11+
# Adding Codac
12+
13+
# In case you installed Codac in a local directory, you need
14+
# to specify its path with the CMAKE_PREFIX_PATH option.
15+
# set(CMAKE_PREFIX_PATH "~/codac/build_install")
16+
17+
find_package(CODAC REQUIRED)
18+
message(STATUS "Found Codac version ${CODAC_VERSION}")
19+
20+
# Initializating Ibex
21+
22+
ibex_init_common()
23+
24+
# Compilation
25+
26+
if(FAST_RELEASE)
27+
add_compile_definitions(FAST_RELEASE)
28+
message(STATUS "You are running Codac in fast release mode. (option -DCMAKE_BUILD_TYPE=Release is required)")
29+
endif()
30+
31+
add_executable(${PROJECT_NAME} main.cpp)
32+
target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS})
33+
target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS})
34+
target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES})

‎examples/14_lohner/main.cpp‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
#include <codac>
2+
3+
using namespace std;
4+
using namespace codac2;
5+
6+
int main()
7+
{
8+
IntervalVector x0({Interval(1,1),Interval(1,1)});
9+
x0.inflate(0.1);
10+
11+
auto ta =create_tdomain({0.,10.0},0.2);
12+
SlicedTube a(ta, IntervalVector(2));
13+
a.set(x0,0.0);
14+
15+
auto tb =create_tdomain({0.,10.0},0.01);
16+
SlicedTube b(tb, IntervalVector(2));
17+
b.set(x0,0.0);
18+
19+
VectorVar x(2);
20+
AnalyticFunction f({x}, {-x[0],-sin(x[1])});
21+
CtcLohner ctc_lohner(f);
22+
23+
ctc_lohner.contract(a);
24+
ctc_lohner.contract(b);
25+
26+
Figure2D fig1 ("Lohner_1",GraphicOutput::VIBES|GraphicOutput::IPE);
27+
Figure2D fig2 ("Lohner_2",GraphicOutput::VIBES|GraphicOutput::IPE);
28+
29+
fig1.plot_tube(a[0]);
30+
fig1.plot_tube(b[0],{Color::blue(),Color::blue()});
31+
32+
fig2.plot_tube(a[1]);
33+
fig2.plot_tube(b[1],{Color::blue(),Color::blue()});
34+
35+
fig1.set_window_properties({100,50},{1000,400});
36+
fig2.set_window_properties({100,500},{1000,400});
37+
}

‎examples/14_lohner/main.py‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
from codac import *
2+
3+
x0 = IntervalVector([Interval(1,1),Interval(1,1)])
4+
x0.inflate(0.1)
5+
6+
ta = create_tdomain([0.,10.0],0.2)
7+
a = SlicedTube_IntervalVector(ta, IntervalVector(2))
8+
a.set(x0,0.0)
9+
10+
tb =create_tdomain([0.,10.0],0.01)
11+
b = SlicedTube_IntervalVector(tb, IntervalVector(2))
12+
b.set(x0,0.0)
13+
14+
x = VectorVar(2)
15+
f = AnalyticFunction([x], [-x[0],-sin(x[1])])
16+
ctc_lohner = CtcLohner(f)
17+
18+
ctc_lohner.contract(a)
19+
ctc_lohner.contract(b)
20+
21+
fig1 = Figure2D("Lohner_1",GraphicOutput.VIBES|GraphicOutput.IPE)
22+
fig2 = Figure2D("Lohner_2",GraphicOutput.VIBES|GraphicOutput.IPE)
23+
24+
fig1.plot_tube(a[0])
25+
fig1.plot_tube(b[0],StyleProperties([Color.blue(),Color.blue()]))
26+
27+
fig2.plot_tube(a[1])
28+
fig2.plot_tube(b[1],StyleProperties([Color.blue(),Color.blue()]))
29+
30+
fig1.set_window_properties([100,50],[1000,400])
31+
fig2.set_window_properties([100,500],[1000,400])

‎python/src/core/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
contractors/codac2_py_CtcInverse.h
2929
contractors/codac2_py_CtcInverseNotIn.h
3030
contractors/codac2_py_CtcLazy.cpp
31+
contractors/codac2_py_CtcLohner.cpp
3132
contractors/codac2_py_CtcNot.cpp
3233
contractors/codac2_py_CtcPointCloud.cpp
3334
contractors/codac2_py_CtcPolar.cpp

‎python/src/core/codac2_py_core.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@ void export_CtcIdentity(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcI
4848
void export_CtcInnerOuter(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcIntervalVector>& ctc);
4949
void export_CtcInter(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcIntervalVector>& ctc);
5050
void export_CtcLazy(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcIntervalVector>& ctc);
51+
void export_CtcLohner(py::module& m);
5152
void export_CtcNot(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcIntervalVector>& ctc);
5253
void export_CtcPointCloud(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcIntervalVector>& ctc);
5354
void export_CtcPolar(py::module& m, py::class_<CtcBase<IntervalVector>,pyCtcIntervalVector>& ctc);
@@ -197,6 +198,7 @@ PYBIND11_MODULE(_core, m)
197198
export_CtcInverseNotIn<ScalarType>(m,"CtcInverseNotIn_Interval",py_ctc_iv);
198199
export_CtcInverseNotIn<VectorType>(m,"CtcInverseNotIn_IntervalVector",py_ctc_iv);
199200
export_CtcLazy(m, py_ctc_iv);
201+
export_CtcLohner(m);
200202
export_CtcNot(m, py_ctc_iv);
201203
export_CtcPointCloud(m, py_ctc_iv);
202204
export_CtcPolar(m, py_ctc_iv);

‎python/src/core/contractors/codac2_py_CtcLohner.cpp‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
/**
2+
* Codac binding (core)
3+
* ----------------------------------------------------------------------------
4+
* \date 2026
5+
* \author Auguste Bourgois, Simon Rohou, Maël Godard
6+
* \copyright Copyright 2024 Codac Team
7+
* \license GNU Lesser General Public License (LGPL)
8+
*/
9+
10+
#include <pybind11/pybind11.h>
11+
#include <pybind11/operators.h>
12+
#include <pybind11/stl.h>
13+
#include <codac2_template_tools.h>
14+
#include <codac2_SlicedTube.h>
15+
#include <codac2_CtcLohner.h>
16+
#include "codac2_py_Ctc.h"
17+
#include "codac2_py_CtcLohner_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py):
18+
#include "codac2_py_cast.h"
19+
20+
using namespace std;
21+
using namespace codac2;
22+
namespace py = pybind11;
23+
using namespace pybind11::literals;
24+
25+
void export_CtcLohner(py::module& m)
26+
{
27+
py::class_<CtcLohner> exported(m, "CtcLohner", CTCLOHNER_MAIN);
28+
exported
29+
30+
.def(py::init(
31+
[](const py::object& f, int contractions, double eps)
32+
{
33+
if(!is_instance<AnalyticFunction<VectorType>>(f)) {
34+
assert_release("CtcLohner: invalid function type");
35+
}
36+
return std::make_unique<CtcLohner>(cast<AnalyticFunction<VectorType>>(f), contractions, eps);
37+
}
38+
),
39+
CTCLOHNER_CTCLOHNER_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_INT_DOUBLE,
40+
"f"_a, "contractions"_a = 5, "eps"_a = 0.1)
41+
42+
.def("contract", [](const CtcLohner& c, py::object& tube, TimePropag t_propa) -> py::object&
43+
{
44+
if(!is_instance<SlicedTube<IntervalVector>>(tube)) {
45+
assert_release("contract: invalid tube type");
46+
}
47+
c.contract(cast<SlicedTube<IntervalVector>>(tube), t_propa);
48+
return tube;
49+
},
50+
VOID_CTCLOHNER_CONTRACT_SLICEDTUBE_INTERVALVECTOR_REF_TIMEPROPAG_CONST,
51+
"tube"_a, "t_propa"_a = TimePropag::FWD_BWD)
52+
;
53+
}

0 commit comments

Comments
 (0)