diff --git a/doc/manual/manual/intervals/Interval_class.rst b/doc/manual/manual/intervals/Interval_class.rst index 86d6048c7..b46da7d2c 100644 --- a/doc/manual/manual/intervals/Interval_class.rst +++ b/doc/manual/manual/intervals/Interval_class.rst @@ -143,6 +143,8 @@ You can access key interval properties: x.diam % diameter x.mag % magnitude x.mig % mignitude + x.smag % signed magnitude + x.smig % signed mignitude x.size % dimension (always 1) diff --git a/doc/manual/manual/intervals/src.cpp b/doc/manual/manual/intervals/src.cpp index 87abf39d9..e057d7cf7 100644 --- a/doc/manual/manual/intervals/src.cpp +++ b/doc/manual/manual/intervals/src.cpp @@ -60,6 +60,8 @@ TEST_CASE("Interval class - manual") x.diam(); // diameter x.mag(); // magnitude x.mig(); // mignitude + x.smag(); // signed magnitude + x.smig(); // signed mignitude x.size(); // dimension (always 1) // [interval-class-4-end] } diff --git a/doc/manual/manual/intervals/src.py b/doc/manual/manual/intervals/src.py index b04817bd6..4683a13bd 100644 --- a/doc/manual/manual/intervals/src.py +++ b/doc/manual/manual/intervals/src.py @@ -48,6 +48,8 @@ def tests_Interval_manual(test): x.diam() # diameter x.mag() # magnitude x.mig() # mignitude + x.smag() # signed magnitude + x.smig() # signed mignitude x.size() # dimension (always 1) # [interval-class-4-end] diff --git a/doc/manual/manual/visualization/colors.rst b/doc/manual/manual/visualization/colors.rst index 4206031ff..b6cfc6da6 100644 --- a/doc/manual/manual/visualization/colors.rst +++ b/doc/manual/manual/visualization/colors.rst @@ -73,7 +73,7 @@ Available line styles are: .. code-tab:: c++ fig.draw_box({{2.2,2.5},{2.2,2.5}}, StyleProperties(Color::red(), "..", "layer1")); // Red edge, dotted line, line width of 0.1 and layer1 - // fig.draw_box({{2.2,2.5},{2.2,2.5}}, {Color.red(), "..", "layer1", "0.1"}); //equivalent + // fig.draw_box({{2.2,2.5},{2.2,2.5}}, {Color::red(), "..", "layer1", "0.1"}); //equivalent Colors ------ diff --git a/examples/07_centered_2D/main.cpp b/examples/07_centered_2D/main.cpp index af464ecd1..48ca64520 100644 --- a/examples/07_centered_2D/main.cpp +++ b/examples/07_centered_2D/main.cpp @@ -35,7 +35,10 @@ int main(){ Vector inflationbox = cent.rad() + T.rad()*df.rad(); Matrix v (2,3); v << (T.rad()*df.mid()), Vector({ inflationbox[0], 0.0 }), Vector({ 0.0, inflationbox[1] }); - fig4.draw_zonotope({cent.mid(),v},{Color::red(),Color::yellow(0.1)}); + fig4.draw_zonotope({cent.mid(),v},{{Color::red(),Color::yellow(0.1)},"zonotopes"}); + + Parallelepiped p = f1.parallelepiped_eval(T); + fig4.draw_parallelepiped(p,{{Color::black(),Color::green(0.1)},"parallelepipeds"}); } time = time+dt; } diff --git a/examples/11_peibos/CMakeLists.txt b/examples/11_peibos/CMakeLists.txt new file mode 100644 index 000000000..cf9ef35ce --- /dev/null +++ b/examples/11_peibos/CMakeLists.txt @@ -0,0 +1,34 @@ +# ================================================================== +# codac / basics example - cmake configuration file +# ================================================================== + + cmake_minimum_required(VERSION 3.5) + project(codac_example LANGUAGES CXX) + + set(CMAKE_CXX_STANDARD 20) + set(CMAKE_CXX_STANDARD_REQUIRED ON) + +# Adding Codac + + # In case you installed Codac in a local directory, you need + # to specify its path with the CMAKE_PREFIX_PATH option. + # set(CMAKE_PREFIX_PATH "~/codac/build_install") + + find_package(CODAC REQUIRED) + message(STATUS "Found Codac version ${CODAC_VERSION}") + +# Initializating Ibex + + ibex_init_common() + +# 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}) \ No newline at end of file diff --git a/examples/11_peibos/main.cpp b/examples/11_peibos/main.cpp new file mode 100644 index 000000000..fbbf1ec51 --- /dev/null +++ b/examples/11_peibos/main.cpp @@ -0,0 +1,105 @@ +#include + +using namespace std; +using namespace codac2; + +int main() +{ + // 2D example of the PEIBOS algorithm + VectorVar y_2d(2); + double a = 1.4; double b = 0.3; + AnalyticFunction f_2d({y_2d},{y_2d[1]+1-a*sqr(y_2d[0]),b*y_2d[0]}); + + VectorVar X_2d(1); + AnalyticFunction psi0_2d ({X_2d},{cos(X_2d[0]*PI/4.-PI/2),sin(X_2d[0]*PI/4.-PI/2)}); + + OctaSym id_2d ({1,2}); + OctaSym s ({-2,1}); + + auto v_par_2d = PEIBOS(f_2d, psi0_2d, {id_2d,s,s*s,s.invert()}, 0.2, {-0.2,0.}, true); + + Figure2D figure_2d ("Henon Map", GraphicOutput::VIBES); + figure_2d.set_window_properties({25,50},{500,500}); + figure_2d.set_axes({0,{-1.4,2.2}}, {1,{-0.4,0.3}}); + + for (const auto& p : v_par_2d) + { + + figure_2d.draw_parallelepiped(p, {Color::green(),Color::green(0.5)}); + figure_2d.draw_box(p.box(), {Color::blue()}); + for (const auto& vertice : p.vertices()) + figure_2d.draw_point(vertice, {Color::red(),Color::red(0.5)}); + } + + // 3D example of the PEIBOS algorithm + VectorVar y_3d(3); + AnalyticFunction f_3d({y_3d},{sqr(y_3d[0])-sqr(y_3d[1])+y_3d[0],2*y_3d[0]*y_3d[1]+y_3d[1],y_3d[2]}); + + VectorVar X_3d(2); + AnalyticFunction psi0_3d ({X_3d},{1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))}); + + OctaSym id_3d ({1,2,3}); + OctaSym s1 ({-2,1,3}); + OctaSym s2 ({3,2,-1}); + + Figure3D figure3d ("Conform"); + figure3d.draw_axes(); + + Figure2D figure_3d_proj ("Conform projected", GraphicOutput::VIBES); + figure_3d_proj.set_window_properties({25,600},{500,500}); + figure_3d_proj.set_axes({0,{-1.5,2.5}}, {1,{-2,2}}); + + auto v_par_3d = PEIBOS(f_3d, psi0_3d, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 0.2, true); + + for (const auto& p : v_par_3d) + { + figure3d.draw_parallelepiped(p, Color::green(0.5)); + figure_3d_proj.draw_zonotope(p.proj({0,1}) , {Color::black(),Color::green(0.2)}); + } + + // nD example of the PEIBOS algorithm + + VectorVar y_nd(3); + Matrix rot_matrix_1 ({ {1,0,0}, + {0,1/std::sqrt(2.0),-1/std::sqrt(2.0)}, + {0,1/std::sqrt(2.0),+1/std::sqrt(2.0)} }); + Matrix rot_matrix_2 ({ {1/std::sqrt(2.0),-1/std::sqrt(2.0),0}, + {1/std::sqrt(2.0),+1/std::sqrt(2.0),0}, + {0,0,1} }); + + // Function mapping cube to sphere + AnalyticFunction g_nd ({y_nd}, {y_nd[0]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[1]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[2]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2]))}); + // Apply two rotations + AnalyticFunction f_nd ({y_nd}, rot_matrix_1 * rot_matrix_2 * g_nd(y_nd)); + + VectorVar X_nd(1); + AnalyticFunction psi0_nd ({X_nd},{X_nd[0],1,1}); + + OctaSym id_nd ({1,2,3}); + OctaSym s1_nd ({-2,1,3}); + OctaSym s2_nd ({3,2,-1}); + + Figure3D figure_3d_nd ("Cube on Sphere"); + figure_3d_nd.draw_axes(0.5); + + Figure2D figure_2d_nd_xy ("XY Plane", GraphicOutput::VIBES); + figure_2d_nd_xy.set_window_properties({575,50},{500,500}); + figure_2d_nd_xy.set_axes(axis(0,{-1.,1.}), axis(1,{-1.,1.})); + + Figure2D figure_2d_nd_zy ("ZY Plane", GraphicOutput::VIBES); + figure_2d_nd_zy.set_window_properties({1125,50},{500,500}); + figure_2d_nd_zy.set_axes(axis(0,{-1.,1.}), axis(1,{-1.,1.})); + + // 12 symmetries are needed : id, s1, s1^2, s1^-1, the four same multiplied by s2, and the four same multiplied by s2^2 + auto v_par_nd = PEIBOS(f_nd, psi0_nd, {id_nd,s1_nd,s1_nd*s1_nd,s1_nd.invert(), + s2_nd,s2_nd*s1_nd,s2_nd*s1_nd*s1_nd,s2_nd*s1_nd.invert(), + s2_nd*s2_nd,s2_nd*s2_nd*s1_nd,s2_nd*s2_nd*s1_nd*s1_nd,s2_nd*s2_nd*s1_nd.invert()}, + 0.1, true); + + for (const auto& p : v_par_nd) + { + figure_3d_nd.draw_parallelepiped(p, Color::green(0.5)); + figure_2d_nd_xy.draw_zonotope(p.proj({0,1}), {Color::black(),Color::green(0.2)}); + figure_2d_nd_zy.draw_zonotope(p.proj({2,1}), {Color::black(),Color::green(0.2)}); + } +} \ No newline at end of file diff --git a/examples/11_peibos/main.py b/examples/11_peibos/main.py new file mode 100644 index 000000000..6e4fd52b4 --- /dev/null +++ b/examples/11_peibos/main.py @@ -0,0 +1,88 @@ +from codac import * +import numpy as np + +if __name__=="__main__": + + # 2D example of the PEIBOS algorithm + + y_2d = VectorVar(2) + a,b = 1.4,0.3 + f_2d = AnalyticFunction([y_2d],[y_2d[1]+1-a*sqr(y_2d[0]),b*y_2d[0]]) + + X_2d = VectorVar(1) + psi0_2d = AnalyticFunction([X_2d],[cos(X_2d[0]*PI/4.-PI/2),sin(X_2d[0]*PI/4.-PI/2)]) + + id_2d = OctaSym([1, 2]) + s = OctaSym([-2, 1]) + + v_par_2d = PEIBOS(f_2d,psi0_2d,[id_2d,s,s*s,s.invert()],0.2,[-0.2,0.],True) + + figure_2d = Figure2D("Henon Map", GraphicOutput.VIBES) + figure_2d.set_window_properties([25,50],[500,500]) + figure_2d.set_axes(axis(0,[-1.4,2.2]), axis(1,[-0.4,0.3])) + + for par in v_par_2d: + figure_2d.draw_parallelepiped(par,[Color.green(),Color.green(0.5)]) + figure_2d.draw_box(par.box(), [Color.blue()]) + for vertice in par.vertices(): + figure_2d.draw_point(vertice, [Color.red(), Color.red(0.5)]) + + # 3D example of the PEIBOS algorithm + + y_3d = VectorVar(3) + f_3d = AnalyticFunction([y_3d],[sqr(y_3d[0])-sqr(y_3d[1])+y_3d[0],2*y_3d[0]*y_3d[1]+y_3d[1],y_3d[2]]) + + X_3d = VectorVar(2) + psi0_3d = AnalyticFunction([X_3d],[1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))]) + + id_3d = OctaSym([1, 2, 3]) + s1 = OctaSym([-2, 1, 3]) + s2 = OctaSym([3, 2, -1]) + + figure_3d = Figure3D("Conform") + figure_3d.draw_axes() + + figure_3d_proj = Figure2D("Conform projected", GraphicOutput.VIBES) + figure_3d_proj.set_window_properties([25,600],[500,500]) + figure_3d_proj.set_axes(axis(0,[-1.5,2.5]), axis(1,[-2,2])) + + v_par_3d = PEIBOS(f_3d,psi0_3d,[id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()],0.2,True) + + for p in v_par_3d: + figure_3d.draw_parallelepiped(p,Color.green(0.5)) + figure_3d_proj.draw_zonotope(p.proj([0, 1]) , [Color.black(),Color.green(0.2)]) + + # nD example of the PEIBOS algorithm + + y_nd = VectorVar(3) + rot_matrix_1 = Matrix([[1,0,0],[0,1/np.sqrt(2.0),-1/np.sqrt(2.0)],[0,1/np.sqrt(2.0),+1/np.sqrt(2.0)]]) + rot_matrix_2 = Matrix([[1/np.sqrt(2.0),-1/np.sqrt(2.0),0],[1/np.sqrt(2.0),1/np.sqrt(2.0),0],[0,0,1]]) + g_nd = AnalyticFunction([y_nd], [y_nd[0]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[1]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[2]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2]))]) + f_nd = AnalyticFunction([y_nd], rot_matrix_1 * rot_matrix_2 * g_nd(y_nd)) + + X_nd = VectorVar(1) + psi0_nd = AnalyticFunction([X_nd],[X_nd[0],1,1]) + + id_nd = OctaSym([1, 2, 3]) + s1_nd = OctaSym([-2, 1, 3]) + s2_nd = OctaSym([3, 2, -1]) + + figure_3d_nd = Figure3D("Cube on Sphere") + figure_3d_nd.draw_axes(0.5) + + figure_2d_nd_xy = Figure2D("XY Plane", GraphicOutput.VIBES) + figure_2d_nd_xy.set_window_properties([575,50],[500,500]) + figure_2d_nd_xy.set_axes(axis(0,[-1., 1.]), axis(1,[-1., 1.])) + + figure_2d_nd_zy = Figure2D("ZY Plane", GraphicOutput.VIBES) + figure_2d_nd_zy.set_window_properties([1125,50],[500,500]) + figure_2d_nd_zy.set_axes(axis(0,[-1., 1.]), axis(1,[-1., 1.])) + + v_par_nd = PEIBOS(f_nd,psi0_nd,[id_nd,s1_nd,s1_nd*s1_nd,s1_nd.invert(), + s2_nd,s2_nd*s1_nd,s2_nd*s1_nd*s1_nd,s2_nd*s1_nd.invert(), + s2_nd*s2_nd,s2_nd*s2_nd*s1_nd,s2_nd*s2_nd*s1_nd*s1_nd,s2_nd*s2_nd*s1_nd.invert()],0.1,True) + + for p in v_par_nd: + figure_3d_nd.draw_parallelepiped(p, Color.green(0.5)) + figure_2d_nd_xy.draw_zonotope(p.proj([0, 1]), [Color.black(), Color.green(0.2)]) + figure_2d_nd_zy.draw_zonotope(p.proj([2, 1]), [Color.black(), Color.green(0.2)]) diff --git a/python/codac/core/__init__.py b/python/codac/core/__init__.py index 2c5f7b0f1..0011571ad 100644 --- a/python/codac/core/__init__.py +++ b/python/codac/core/__init__.py @@ -68,6 +68,9 @@ def traj_eval(self,*args): def tube_eval(self,*args): return self.f.tube_eval(*args) + + def parallelepiped_eval(self,*args): + return self.f.parallelepiped_eval(*args) def diff(self,*args): return self.f.diff(*args) diff --git a/python/src/core/CMakeLists.txt b/python/src/core/CMakeLists.txt index 6df10a790..79eaccc17 100644 --- a/python/src/core/CMakeLists.txt +++ b/python/src/core/CMakeLists.txt @@ -86,6 +86,8 @@ paver/codac2_py_pave.cpp + peibos/codac2_py_peibos.cpp + separators/codac2_py_Sep.cpp separators/codac2_py_Sep.h separators/codac2_py_SepAction.cpp @@ -131,6 +133,7 @@ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/graphics/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/matrices/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/paver/ + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/peibos/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/separators/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/tools/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/trajectory/ diff --git a/python/src/core/codac2_py_core.cpp b/python/src/core/codac2_py_core.cpp index 21d36d84c..bd68495e2 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -120,6 +120,9 @@ void export_operators(py::module& m); // paver void export_pave(py::module& m); +// peibos +void export_peibos(py::module& m); + // separators py::class_ export_Sep(py::module& m); void export_SepAction(py::module& m, py::class_& pysep); @@ -273,6 +276,9 @@ PYBIND11_MODULE(_core, m) // paver export_pave(m); + // peibos + export_peibos(m); + // separators auto py_sep = export_Sep(m); export_SepAction(m,py_sep); diff --git a/python/src/core/domains/interval/codac2_py_Interval.cpp b/python/src/core/domains/interval/codac2_py_Interval.cpp index fd01bcc92..b70939aec 100644 --- a/python/src/core/domains/interval/codac2_py_Interval.cpp +++ b/python/src/core/domains/interval/codac2_py_Interval.cpp @@ -87,6 +87,12 @@ py::class_ export_Interval(py::module& m) .def("mig", &Interval::mig, DOUBLE_INTERVAL_MIG_CONST) + .def("smag", &Interval::smag, + DOUBLE_INTERVAL_SMAG_CONST) + + .def("smig", &Interval::smig, + DOUBLE_INTERVAL_SMIG_CONST) + .def("rand", &Interval::rand, DOUBLE_INTERVAL_RAND_CONST) diff --git a/python/src/core/domains/interval/codac2_py_IntervalMatrixBase.h b/python/src/core/domains/interval/codac2_py_IntervalMatrixBase.h index a3dea9e86..806f3b0c3 100644 --- a/python/src/core/domains/interval/codac2_py_IntervalMatrixBase.h +++ b/python/src/core/domains/interval/codac2_py_IntervalMatrixBase.h @@ -55,6 +55,12 @@ void export_IntervalMatrixBase(py::module& m, py::class_& pyclass) .def("mig", [](const S& x) { return x.mig(); }, MATRIXBASE_ADDONS_INTERVALMATRIXBASE_AUTO_MIG_CONST) + .def("smag", [](const S& x) { return x.smag(); }, + MATRIXBASE_ADDONS_INTERVALMATRIXBASE_AUTO_SMAG_CONST) + + .def("smig", [](const S& x) { return x.smig(); }, + MATRIXBASE_ADDONS_INTERVALMATRIXBASE_AUTO_SMIG_CONST) + .def("rand", [](const S& x) { return x.rand(); }, MATRIXBASE_ADDONS_INTERVALMATRIXBASE_AUTO_RAND_CONST) diff --git a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp index 2fae348b8..63b9ff280 100644 --- a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp +++ b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp @@ -32,18 +32,15 @@ void export_Parallelepiped(py::module& m) PARALLELEPIPED_PARALLELEPIPED_CONST_VECTOR_REF_CONST_MATRIX_REF, "z"_a, "A"_a) - .def("proj",[](const Parallelepiped& x, const std::vector& indices) - { - return x.proj(matlab::convert_indices(indices)); - }, - ZONOTOPE_PARALLELEPIPED_PROJ_CONST_VECTOR_INDEX_REF_CONST, - "indices"_a) - .def("vertices", &Parallelepiped::vertices, VECTOR_VECTOR_PARALLELEPIPED_VERTICES_CONST) .def("box", &Parallelepiped::box, INTERVALVECTOR_PARALLELEPIPED_BOX_CONST) + .def("contains", &Parallelepiped::contains, + BOOL_PARALLELEPIPED_CONTAINS_CONST_VECTOR_REF_CONST, + "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 a4b7d362c..20595f94c 100644 --- a/python/src/core/domains/zonotope/codac2_py_Zonotope.cpp +++ b/python/src/core/domains/zonotope/codac2_py_Zonotope.cpp @@ -30,6 +30,13 @@ void export_Zonotope(py::module& m) ZONOTOPE_ZONOTOPE_CONST_VECTOR_REF_CONST_MATRIX_REF, "z"_a, "A"_a) + .def("proj",[](const Zonotope& x, const std::vector& indices) + { + return x.proj(matlab::convert_indices(indices)); + }, + ZONOTOPE_ZONOTOPE_PROJ_CONST_VECTOR_INDEX_REF_CONST, + "indices"_a) + .def_readwrite("z", &Zonotope::z, VECTOR_ZONOTOPE_Z) diff --git a/python/src/core/functions/analytic/codac2_py_AnalyticFunction.h b/python/src/core/functions/analytic/codac2_py_AnalyticFunction.h index 2256aaa0e..28e01214a 100644 --- a/python/src/core/functions/analytic/codac2_py_AnalyticFunction.h +++ b/python/src/core/functions/analytic/codac2_py_AnalyticFunction.h @@ -321,12 +321,7 @@ void export_AnalyticFunction(py::module& m, const std::string& export_name) if constexpr(std::is_same_v) { - exported - - .def("parallelepiped_eval", &AnalyticFunction::parallelepiped_eval, - PARALLELEPIPED_ANALYTICFUNCTION_T_PARALLELEPIPED_EVAL_CONST_INTERVALVECTOR_REF_CONST, - "x"_a) - ; + bind_(exported, "parallelepiped_eval", parallelepiped_eval, PARALLELEPIPED_ANALYTICFUNCTION_T_PARALLELEPIPED_EVAL_CONST_ARGS_REF_VARIADIC_CONST); } exported diff --git a/python/src/core/peibos/codac2_py_peibos.cpp b/python/src/core/peibos/codac2_py_peibos.cpp new file mode 100644 index 000000000..64923ac5e --- /dev/null +++ b/python/src/core/peibos/codac2_py_peibos.cpp @@ -0,0 +1,52 @@ +/** + * Codac binding (unsupported) + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include + #include +#include +#include +#include +#include +#include "codac2_py_peibos_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_peibos(py::module& m) +{ + m.def("parallelepiped_inclusion", + [](const IntervalVector& z, const IntervalMatrix& Jf, const Matrix& Jf_tild, const py::object& psi_0, const OctaSym& sigma, const IntervalVector& X) + { + return parallelepiped_inclusion(z, Jf, Jf_tild, cast>(psi_0), sigma, X); + }, + PARALLELEPIPED_PARALLELEPIPED_INCLUSION_CONST_INTERVALVECTOR_REF_CONST_INTERVALMATRIX_REF_CONST_MATRIX_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_OCTASYM_REF_CONST_INTERVALVECTOR_REF, + "z"_a, "Jf"_a, "Jf_tild"_a, "psi_0"_a, "sigma"_a, "X"_a); + + + m.def("PEIBOS", + [](const py::object& f, const py::object& psi_0, const vector& Sigma, double epsilon, bool verbose = false) + { + return PEIBOS(cast>(f), cast>(psi_0), Sigma, epsilon, verbose); + }, + VECTOR_PARALLELEPIPED_PEIBOS_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_VECTOR_OCTASYM_REF_DOUBLE_BOOL, + "f"_a, "psi_0"_a, "Sigma"_a, "epsilon"_a, "verbose"_a = false); + + m.def("PEIBOS", + [](const py::object& f, const py::object& psi_0, const vector& Sigma, double epsilon, const Vector& offset, bool verbose = false) + { + return PEIBOS(cast>(f), cast>(psi_0), Sigma, epsilon, offset, verbose); + }, + VECTOR_PARALLELEPIPED_PEIBOS_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_VECTOR_OCTASYM_REF_DOUBLE_CONST_VECTOR_REF_BOOL, + "f"_a, "psi_0"_a, "Sigma"_a, "epsilon"_a, "offset"_a, "verbose"_a = false); +} \ No newline at end of file diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index c6d4104c6..ae0c73164 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -193,6 +193,9 @@ ${CMAKE_CURRENT_SOURCE_DIR}/paver/codac2_pave.cpp ${CMAKE_CURRENT_SOURCE_DIR}/paver/codac2_pave.h + ${CMAKE_CURRENT_SOURCE_DIR}/peibos/codac2_peibos.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/peibos/codac2_peibos.h + ${CMAKE_CURRENT_SOURCE_DIR}/proj/codac2_ProjBase.cpp ${CMAKE_CURRENT_SOURCE_DIR}/proj/codac2_ProjBase.h @@ -285,6 +288,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}/matrices/eigen/MatrixBase_addons ${CMAKE_CURRENT_SOURCE_DIR}/operators ${CMAKE_CURRENT_SOURCE_DIR}/paver + ${CMAKE_CURRENT_SOURCE_DIR}/peibos ${CMAKE_CURRENT_SOURCE_DIR}/proj ${CMAKE_CURRENT_SOURCE_DIR}/separators ${CMAKE_CURRENT_SOURCE_DIR}/tools diff --git a/src/core/domains/codac2_Domain.h b/src/core/domains/codac2_Domain.h index 61b6c6142..f7827d96e 100644 --- a/src/core/domains/codac2_Domain.h +++ b/src/core/domains/codac2_Domain.h @@ -39,6 +39,8 @@ namespace codac2 virtual V mid() const = 0; virtual V mag() const = 0; virtual V mig() const = 0; + virtual V smag() const = 0; + virtual V smig() const = 0; virtual V rad() const = 0; virtual V diam() const = 0; virtual double volume() const = 0; diff --git a/src/core/domains/interval/codac2_Interval.h b/src/core/domains/interval/codac2_Interval.h index 85006cb31..b8b1e2a7f 100644 --- a/src/core/domains/interval/codac2_Interval.h +++ b/src/core/domains/interval/codac2_Interval.h @@ -210,6 +210,25 @@ namespace codac2 */ double mig() const; + /** + * \brief Returns the signed magnitude of this + * i.e. lower bound if |lower bound|>|upper bound|, upper bound otherwise. + * + * \return the signed magnitude of the interval + */ + double smag() const; + + /** + * \brief Returns the signed mignitude of this + * + * lower bound if lower_bound > 0, + * upper bound if upper_bound < 0, + * 0 otherwise. + * + * \return the signed mignitude of the interval + */ + double smig() const; + /** * \brief Returns a random value inside the interval * diff --git a/src/core/domains/interval/codac2_Interval_impl.h b/src/core/domains/interval/codac2_Interval_impl.h index dc7b48dd0..5a8f29f95 100644 --- a/src/core/domains/interval/codac2_Interval_impl.h +++ b/src/core/domains/interval/codac2_Interval_impl.h @@ -134,6 +134,16 @@ namespace codac2 return gaol::interval::mig(); } + inline double Interval::smag() const + { + return (abs(lb()) > abs(ub())) ? lb() : ub(); + } + + inline double Interval::smig() const + { + return gaol::interval::smig(); + } + inline double Interval::rand() const { if(is_empty()) diff --git a/src/core/domains/zonotope/codac2_Parallelepiped.cpp b/src/core/domains/zonotope/codac2_Parallelepiped.cpp index 8a1a40f66..4073b57b4 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped.cpp +++ b/src/core/domains/zonotope/codac2_Parallelepiped.cpp @@ -18,23 +18,6 @@ Parallelepiped::Parallelepiped(const Vector& z_, const Matrix& A_) assert_release(A.cols() <= z.size() && "too many vectors, you are describing a zonotope"); } -Zonotope Parallelepiped::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"); - - Matrix A_cropped (indices.size(), A.cols()); - Vector z_cropped (indices.size()); - - for (size_t i = 0; i < indices.size(); ++i) - { - A_cropped.row(i) = A.row(indices[i]); - z_cropped[i] = z[indices[i]]; - } - - return Zonotope(z_cropped, A_cropped); -} - void generate_vertices(Index i, Index n, const Vector& z, const Matrix& A, vector& L_v) { if (i == n) @@ -61,4 +44,13 @@ IntervalVector Parallelepiped::box() const for(const auto& v : vertices()) box |= v; return box; +} + +bool Parallelepiped::contains(const Vector& v) const +{ + assert_release(A.rows() == A.cols() && "Matrix A must be square to check containment."); + + IntervalVector IV = Interval(-1,1)*IntervalVector::Ones(A.cols()); + + return IV.contains(A.inverse()*(v - z)); } \ No newline at end of file diff --git a/src/core/domains/zonotope/codac2_Parallelepiped.h b/src/core/domains/zonotope/codac2_Parallelepiped.h index 95960bea8..3fb791008 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped.h @@ -37,15 +37,6 @@ namespace codac2 * \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); - - /** - * \brief Projects the parallelepiped onto the subspace defined by the given indices - * - * \param indices Vector of indices of the dimensions to project onto - * - * \return A new Zonotope object representing the projection of the parallelepiped onto the specified subspace - */ - Zonotope proj(const std::vector& indices) const; /** * \brief Computes the vertices of the parallelepiped @@ -61,5 +52,14 @@ namespace codac2 */ IntervalVector box() const; + /** + * \brief Checks if a given point is contained within the parallelepiped. The matrix A has to be square and invertible. + * + * \param v The point to check (n-dimensional vector) + * + * \return true if the point is inside the parallelepiped, false otherwise + */ + bool contains(const Vector& v) const; + }; } diff --git a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h index bfa6c8420..42fbab5b5 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h @@ -11,19 +11,32 @@ #include "codac2_AnalyticFunction.h" #include "codac2_Parallelepiped.h" +#include "codac2_peibos.h" namespace codac2 { template requires std::is_base_of_v - inline Parallelepiped AnalyticFunction::parallelepiped_eval(const IntervalVector& x) const + template + inline Parallelepiped AnalyticFunction::parallelepiped_eval(const Args&... x) const { - this->check_valid_inputs(x); - assert_release(this->input_size() < this->output_size()); - assert_release(this->args().size() == 1 && "f must have only one arg"); + this->check_valid_inputs(x...); + assert_release(this->input_size() < this->output_size() && + "Parallelepiped evaluation requires more outputs than inputs."); + assert_release(this->input_size() > 0 && + "Parallelepiped evaluation requires at least one input."); - // todo + auto Y = this->eval(x.mid()...); + auto z = Y.mid(); - return Parallelepiped({{0.}},{{0.}}); + Matrix A = this->diff(x.mid()...).mid(); + + // Maximum error computation + double rho = error(Y, z, this->diff(x...), A, cart_prod(x...)); + + // Inflation of the parallelepiped + auto A_inf = inflate_flat_parallelepiped(A, cart_prod(x...).rad(), rho); + + return Parallelepiped(z, A_inf); } } \ No newline at end of file diff --git a/src/core/domains/zonotope/codac2_Zonotope.cpp b/src/core/domains/zonotope/codac2_Zonotope.cpp index cc06b7e37..d6a3538ea 100644 --- a/src/core/domains/zonotope/codac2_Zonotope.cpp +++ b/src/core/domains/zonotope/codac2_Zonotope.cpp @@ -18,3 +18,20 @@ Zonotope::Zonotope(const Vector& z_, const Matrix& A_) assert_release(z.size() == A.rows()); } +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"); + + Matrix A_cropped (indices.size(), A.cols()); + Vector z_cropped (indices.size()); + + for (size_t i = 0; i < indices.size(); ++i) + { + A_cropped.row(i) = A.row(indices[i]); + z_cropped[i] = z[indices[i]]; + } + + return Zonotope(z_cropped, A_cropped); +} + diff --git a/src/core/domains/zonotope/codac2_Zonotope.h b/src/core/domains/zonotope/codac2_Zonotope.h index 1a147b369..27a0ca212 100644 --- a/src/core/domains/zonotope/codac2_Zonotope.h +++ b/src/core/domains/zonotope/codac2_Zonotope.h @@ -37,6 +37,15 @@ namespace codac2 */ Zonotope(const Vector& z, const Matrix& A); + /** + * \brief Projects the Zonotope onto the subspace defined by the given indices + * + * \param indices Vector of indices of the dimensions to project onto + * + * \return A new Zonotope object representing the projection of the Zonotope onto the specified subspace + */ + Zonotope proj(const std::vector& indices) const; + /** * \brief Center of the zonotope */ diff --git a/src/core/functions/analytic/codac2_AnalyticFunction.h b/src/core/functions/analytic/codac2_AnalyticFunction.h index bd9908bdf..9c9b5651e 100644 --- a/src/core/functions/analytic/codac2_AnalyticFunction.h +++ b/src/core/functions/analytic/codac2_AnalyticFunction.h @@ -195,7 +195,8 @@ namespace codac2 return y; } - Parallelepiped parallelepiped_eval(const IntervalVector& x) const; + template + Parallelepiped parallelepiped_eval(const Args&... x) const; // -> is defined in codac2_Parallelepiped_eval.h file Index output_size() const diff --git a/src/core/matrices/eigen/MatrixBase_addons/codac2_MatrixBase_addons_IntervalMatrixBase.h b/src/core/matrices/eigen/MatrixBase_addons/codac2_MatrixBase_addons_IntervalMatrixBase.h index 1f2b07431..425b4bd28 100644 --- a/src/core/matrices/eigen/MatrixBase_addons/codac2_MatrixBase_addons_IntervalMatrixBase.h +++ b/src/core/matrices/eigen/MatrixBase_addons/codac2_MatrixBase_addons_IntervalMatrixBase.h @@ -145,6 +145,32 @@ inline auto mig() const degenerate_mat(mig); } +/** + * \brief Returns a matrix containing the signed magnitudes of each interval element. + * + * The signed magnitude is lower bound if |lower bound|>|upper bound|, upper bound otherwise. + * + * \return A matrix of doubles with the signed magnitudes of each interval. + */ +template + requires IsIntervalDomain +inline auto smag() const +{ + degenerate_mat(smag); +} + +/** + * \brief Returns a matrix containing the signed mignitudes of each interval element. + * + * \return A matrix of doubles with the signed mignitudes of each interval. + */ +template + requires IsIntervalDomain +inline auto smig() const +{ + degenerate_mat(smig); +} + /** * \brief Returns a matrix with random values chosen inside each interval element. * diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp new file mode 100644 index 000000000..1c554d3ee --- /dev/null +++ b/src/core/peibos/codac2_peibos.cpp @@ -0,0 +1,149 @@ +/** + * codac2_peibos.cpp + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_peibos.h" +#include "codac2_Parallelepiped_eval.h" + +using namespace std; +using namespace codac2; + +namespace codac2 +{ + double split (const IntervalVector& X, double eps, vector& boxes) + { + if (X.max_diam()<=eps) + { + boxes.push_back(X); + return X.max_diam(); + } + else + { + auto p = X.bisect_largest(0.5); + double diam1 = split(p.first,eps,boxes); + double diam2 = split (p.second,eps,boxes); + return std::max(diam1,diam2); + } + } + + double error(const IntervalVector& Y, const Vector& z, const IntervalMatrix& Jf, const Matrix& A, const IntervalVector& X) + { + auto xc = X.mid(); + + auto dX = X-xc; + + auto E = (Y - z) + (Jf - A)*dX; + + return E.norm().ub(); + } + + Matrix inflate_flat_parallelepiped(const Matrix& A, const Vector& e_vec, double rho) + { + Index m = A.cols(); + Index n = A.rows(); + + IntervalMatrix A_int (A); + + IntvFullPivLU LUdec((IntervalMatrix) A_int.transpose()); + IntervalMatrix N = LUdec.kernel(); + + IntervalMatrix A_tild (n,n); + A_tild << A, N; + + IntervalMatrix Q = inverse_enclosure(A_tild.transpose() * A_tild); + + IntervalMatrix mult = IntervalMatrix::Zero(n,n); + for (int i = 0; i < n; i++) + mult(i,i) = rho*sqrt(Q(i,i)); + + for (int i = 0; i < m; i++) + mult(i,i) += e_vec(i); + + // From here we have an IntervalMatrix A_inf, meaning that each generator is an interval vector + IntervalMatrix A_inf = A_tild * mult; + + // The initial parallelepiped is = Y*[-1,1]^n + Matrix Y = A_inf.smig(); + + // The following is similar to the previous operations. The difference is that here the Matrix Y is square and not interval + + // rho2 is the sum of the radiuses of the circle enclosing each interval generator + Interval rho2 (0.); + + for (int i = 0; i < n; i++) + rho2 += A_inf.col(i).rad().norm(); + + IntervalMatrix Q2 = inverse_enclosure(Y.transpose() * Y); + + IntervalMatrix Y2 = IntervalMatrix::Zero(n,n); + + for (int i = 0; i < n; i++) + Y2.col(i) = Y.col(i)*(1+rho2*sqrt(Q2(i,i))); + + return Y2.smag(); + } + + Parallelepiped parallelepiped_inclusion(const IntervalVector& Y, const IntervalMatrix& Jf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X) + { + // Computation of the Jacobian of g = f o sigma(psi_0) + IntervalMatrix Jg = Jf * (sigma.permutation_matrix().template cast()) * psi_0.diff(X); + + Vector z = Y.mid(); + // A is an approximation of the Jacobian of g at the center of X + Matrix A = (Jf_tild * sigma.permutation_matrix() * (psi_0.diff(X.mid()).mid())); + + // Maximum error computation + double rho = error(Y, z, Jg, A, X); + + // Inflation of the parallelepiped + auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); + + return Parallelepiped(z, A_inf); + } + + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, bool verbose) + { + return PEIBOS(f, psi_0, Sigma, epsilon, Vector::Zero(psi_0.output_size()), verbose); + } + + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, const Vector& offset, bool verbose) + { + Index m = psi_0.input_size(); + + assert_release (f.input_size() == psi_0.output_size() && "output size of psi_0 must match input size of f"); + assert_release (offset.size() == psi_0.output_size() && "offset size must match output size of psi_0"); + assert_release (m < psi_0.output_size()); + assert_release (Sigma.size() > 0 && (int) Sigma[0].size() == psi_0.output_size() && "no generator given or wrong dimension of generator (must match output size of psi_0)"); + + clock_t t_start = clock(); + + vector output; + + vector boxes; + double true_eps = split(IntervalVector::constant(m,{-1,1}), epsilon, boxes); + + for (const auto& sigma : Sigma) + { + VectorVar x(m); + AnalyticFunction g_i ({x}, f(sigma(psi_0(x))+offset)); + + for (const auto& X : boxes) + output.push_back(g_i.parallelepiped_eval(X)); + } + + if (verbose) + { + printf("\nPEIBOS statistics:\n"); + printf("------------------\n"); + printf("Real epsilon: %.4f\n", true_eps); + printf("Computation time: %.4fs\n\n", (double)(clock()-t_start)/CLOCKS_PER_SEC); + } + + return output; + } +} \ No newline at end of file diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h new file mode 100644 index 000000000..6c736ed37 --- /dev/null +++ b/src/core/peibos/codac2_peibos.h @@ -0,0 +1,94 @@ +/** + * \file codac2_peibos.h + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include "codac2_Parallelepiped.h" +#include "codac2_AnalyticFunction.h" +#include "codac2_OctaSym.h" +#include "codac2_inversion.h" +#include "codac2_OctaSym_operator.h" +#include "codac2_IntvFullPivLU.h" + + +namespace codac2 +{ + /** + * \brief Recursively split a box until the maximum diameter of each box is less than eps. Note that the resulting boxes are stored in place in the vector boxes. + * + * \param X The box to split + * \param eps The maximum diameter of each box + * \param boxes The vector to store the resulting boxes + * + * \return The maximum diameter of the resulting boxes (lower or equal to eps) + */ + double split (const IntervalVector& X, double eps, vector& boxes); + + /** + * \brief Compute the error term for the parallelepiped inclusion. The error is later used to inflate the flat parallelepiped \f$\mathbf{z}+\mathbf{A}\cdot(\left[\mathbf{x}\right]-\bar{\mathbf{x}})\f$ + * + * \param Y The box enclosing \f$\mathbf{f}(\bar{\mathbf{x}})\f$ + * \param z Any Vector with the same dimension as z. Usually we pick an approximation of \f$\mathbf{f}(\bar{\mathbf{x}})\f$ + * \param Jf The interval Jacobian matrix containing \f$\frac{d\mathbf{f}}{d\mathbf{x}}([\mathbf{x}])\f$ + * \param A Any Matrix with the same dimension as Jf. Usually we pick an approximation of \f$\frac{d\mathbf{f}}{d\mathbf{x}}(\bar{\mathbf{x}})\f$ + * \param X The box \f$[\mathbf{x}]\f$ + */ + double error(const IntervalVector& Y, const Vector& z, const IntervalMatrix& Jf, const Matrix& A, const IntervalVector& X); + + /** + * \brief Inflate the flat parallelepiped \f$\mathbf{A}\cdot e_{vec}\f$ by \f$\rho\f$. + * + * \param A The shape matrix of the flat parallelepiped (n x m matrix with m < n) + * \param e_vec The vector of scaling along each generator of the flat parallelepiped (m-dimensional vector) + * \param rho The inflation factor + * + * \return The shape matrix of the inflated parallelepiped (n x n matrix) + */ + Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); + + /** + * \brief Used in PEIBOS. Compute a parallelepiped enclosing of \f$\mathbf{g}([\mathbf{x}])\f$ where \f$\mathbf{g} = \mathbf{f}\circ \sigma \circ \psi_0\f$. + * + * \param Y The box enclosing \f$\mathbf{f}(\sigma(\psi_0(\bar{\mathbf{x}})))\f$. + * \param Jf The interval Jacobian matrix containing \f$\frac{d\mathbf{f}}{d\mathbf{y}}([\mathbf{y}])\f$ where \f$\mathbf{y} = \sigma(\psi_0(\left[\mathbf{x}\right]))\f$ + * \param Jf_tild An approximation of \f$\frac{d\mathbf{f}}{d\mathbf{y}}(\sigma(\psi_0(\bar{\mathbf{x}})))\f$ + * \param psi_0 The analytic function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param sigma The symmetry operator to construct the atlas + * \param X The box \f$[\mathbf{x}]\f$ + * + * \return A Parallelepiped enclosing \f$\mathbf{g}([\mathbf{x}])\f$ + */ + Parallelepiped parallelepiped_inclusion(const IntervalVector& Y, const IntervalMatrix& Jf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X); + + /** + * \brief Compute a set of parallelepipeds enclosing \f$\mathbf{f}(\sigma(\psi_0([-1,1]^m)))\f$ for each symmetry \f$\sigma\f$ in the set of symmetries \f$\Sigma\f$. Note that \f$\left\{\psi_0,\Sigma\right\}\f$ form a gnomonic atlas. + * + * \param f The analytic function \f$\mathbf{f}:\mathbb{R}^n\rightarrow\mathbb{R}^n\f$ + * \param psi_0 The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param Sigma The set of symmetry operators \f$\sigma\f$ to construct the atlas + * \param epsilon The maximum diameter of the boxes to split \f$[-1,1]^m\f$ before computing the parallelepiped inclusions + * \param verbose If true, print the time taken to compute the parallelepiped inclusions with other statistics + * + * \return A vector of Parallelepipeds enclosing \f$\mathbf{f}(\sigma(\psi_0([-1,1]^m)))\f$ for each symmetry \f$\sigma\f$ in the set of symmetries \f$\Sigma\f$. + */ + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, bool verbose = false); + + /** + * \brief Compute a set of parallelepipeds enclosing \f$\mathbf{f}(\sigma(\psi_0([-1,1]^m)) + offset) \f$ for each symmetry \f$\sigma\f$ in the set of symmetries \f$\Sigma\f$. Note that \f$\left\{\psi_0,\Sigma\right\}\f$ form a gnomonic atlas. + * \param f The analytic function \f$\mathbf{f}:\mathbb{R}^n\rightarrow\mathbb{R}^n\f$ + * \param psi_0 The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param Sigma The set of symmetry operators \f$\sigma\f$ to construct the atlas + * \param epsilon The maximum diameter of the boxes to split \f$[-1,1]^m\f$ before computing the parallelepiped inclusions + * \param offset The offset to add to \f$\sigma(\psi_0([-1,1]^m))\f$ (used to translate the initial manifold) + * \param verbose If true, print the time taken to compute the parallelepiped inclusions with other statistics + * + * \return A vector of Parallelepipeds enclosing \f$\mathbf{f}(\sigma(\psi_0([-1,1]^m))+ offset)\f$ for each symmetry \f$\sigma\f$ in the set of symmetries \f$\Sigma\f$. + */ + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, const Vector& offset, bool verbose = false); +} \ No newline at end of file diff --git a/src/graphics/styles/codac2_ColorMap.h b/src/graphics/styles/codac2_ColorMap.h index 51102f0bc..9f0fb0314 100644 --- a/src/graphics/styles/codac2_ColorMap.h +++ b/src/graphics/styles/codac2_ColorMap.h @@ -165,6 +165,5 @@ namespace codac2 } return cmap; } - }; } \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 90f789147..c34f808ec 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -48,6 +48,7 @@ list(APPEND SRC_TESTS # listing files without extension core/domains/interval/codac2_tests_IntervalMatrix core/domains/interval/codac2_tests_IntervalVector core/domains/zonotope/codac2_tests_Parallelepiped + core/domains/zonotope/codac2_tests_Parallelepiped_eval core/domains/tube/codac2_tests_TDomain core/domains/tube/codac2_tests_SlicedTube core/domains/tube/codac2_tests_SlicedTube_integral @@ -76,6 +77,8 @@ list(APPEND SRC_TESTS # listing files without extension core/operators/codac2_tests_operators + core/peibos/codac2_tests_peibos + core/separators/codac2_tests_SepCartProd core/separators/codac2_tests_SepCtcBoundary core/separators/codac2_tests_SepInverse diff --git a/tests/core/domains/interval/codac2_tests_Interval.cpp b/tests/core/domains/interval/codac2_tests_Interval.cpp index fcf029ab8..7c756c22c 100644 --- a/tests/core/domains/interval/codac2_tests_Interval.cpp +++ b/tests/core/domains/interval/codac2_tests_Interval.cpp @@ -103,6 +103,10 @@ TEST_CASE("Interval - tests from IBEX") CHECK(x.mid() == 1); CHECK(x.rad() == 1); CHECK(x.diam() == 2); + CHECK(x.mag() == 2); + CHECK(x.mig() == 0); + CHECK(x.smag() == 2); + CHECK(x.smig() == 0); x = Interval(-3,-1); @@ -111,6 +115,10 @@ TEST_CASE("Interval - tests from IBEX") CHECK(x.mid() == -2); CHECK(x.rad() == 1); CHECK(x.diam() == 2); + CHECK(x.mag() == 3); + CHECK(x.mig() == 1); + CHECK(x.smag() == -3); + CHECK(x.smig() == -1); x = Interval(-3,1); @@ -119,6 +127,10 @@ TEST_CASE("Interval - tests from IBEX") CHECK(x.mid() == -1); CHECK(x.rad() == 2); CHECK(x.diam() == 4); + CHECK(x.mag() == 3); + CHECK(x.mig() == 0); + CHECK(x.smag() == -3); + CHECK(x.smig() == 0); x = Interval(-oo,0); @@ -127,6 +139,10 @@ TEST_CASE("Interval - tests from IBEX") CHECK(x.mid() == -std::numeric_limits::max()); CHECK(x.rad() == oo); CHECK(x.diam() == oo); + CHECK(x.mag() == oo); + CHECK(x.mig() == 0); + CHECK(x.smag() == -oo); + CHECK(x.smig() == 0); x = Interval(-oo,oo); @@ -135,6 +151,10 @@ TEST_CASE("Interval - tests from IBEX") CHECK(x.mid() == 0); CHECK(x.rad() == oo); CHECK(x.diam() == oo); + CHECK(x.mag() == oo); + CHECK(x.mig() == 0); + CHECK(x.smag() == oo); + CHECK(x.smig() == 0); x = Interval(std::numeric_limits::max(),oo); @@ -143,6 +163,10 @@ TEST_CASE("Interval - tests from IBEX") CHECK(x.mid() == std::numeric_limits::max()); CHECK(x.rad() == oo); CHECK(x.diam() == oo); + CHECK(x.mag() == oo); + CHECK(x.mig() == std::numeric_limits::max()); + CHECK(x.smag() == oo); + CHECK(x.smig() == std::numeric_limits::max()); x = Interval(-1,1); for(size_t i = 0 ; i < 10 ; i++) diff --git a/tests/core/domains/interval/codac2_tests_Interval.py b/tests/core/domains/interval/codac2_tests_Interval.py index 3fe569e0f..460a798fa 100644 --- a/tests/core/domains/interval/codac2_tests_Interval.py +++ b/tests/core/domains/interval/codac2_tests_Interval.py @@ -97,6 +97,10 @@ def test_interval(self): self.assertTrue(x.mid() == 1) self.assertTrue(x.rad() == 1) self.assertTrue(x.diam() == 2) + self.assertTrue(x.mag() == 2) + self.assertTrue(x.mig() == 0) + self.assertTrue(x.smag() == 2) + self.assertTrue(x.smig() == 0) x = Interval(-3,-1) @@ -105,6 +109,10 @@ def test_interval(self): self.assertTrue(x.mid() == -2) self.assertTrue(x.rad() == 1) self.assertTrue(x.diam() == 2) + self.assertTrue(x.mag() == 3) + self.assertTrue(x.mig() == 1) + self.assertTrue(x.smag() == -3) + self.assertTrue(x.smig() == -1) x = Interval(-3,1) @@ -113,6 +121,10 @@ def test_interval(self): self.assertTrue(x.mid() == -1) self.assertTrue(x.rad() == 2) self.assertTrue(x.diam() == 4) + self.assertTrue(x.mag() == 3) + self.assertTrue(x.mig() == 0) + self.assertTrue(x.smag() == -3) + self.assertTrue(x.smig() == 0) x = Interval(-oo,0) @@ -121,6 +133,10 @@ def test_interval(self): self.assertTrue(x.mid() == -sys.float_info.max) self.assertTrue(x.rad() == oo) self.assertTrue(x.diam() == oo) + self.assertTrue(x.mag() == oo) + self.assertTrue(x.mig() == 0) + self.assertTrue(x.smag() == -oo) + self.assertTrue(x.smig() == 0) x = Interval(-oo,oo) @@ -129,6 +145,10 @@ def test_interval(self): self.assertTrue(x.mid() == 0) self.assertTrue(x.rad() == oo) self.assertTrue(x.diam() == oo) + self.assertTrue(x.mag() == oo) + self.assertTrue(x.mig() == 0) + self.assertTrue(x.smag() == oo) + self.assertTrue(x.smig() == 0) x = Interval(sys.float_info.max,oo) @@ -137,17 +157,21 @@ def test_interval(self): self.assertTrue(x.mid() == sys.float_info.max) self.assertTrue(x.rad() == oo) self.assertTrue(x.diam() == oo) + self.assertTrue(x.mag() == oo) + self.assertTrue(x.mig() == sys.float_info.max) + self.assertTrue(x.smag() == oo) + self.assertTrue(x.smig() == sys.float_info.max) x = Interval(-1,1) for i in range(10): self.assertTrue(x.contains(x.rand())) - x = Interval(-oo,0); + x = Interval(-oo,0) for i in range(10): self.assertTrue(x.contains(x.rand())) - x = Interval(0,oo); + x = Interval(0,oo) for i in range(10): self.assertTrue(x.contains(x.rand())) - x = Interval(-oo,oo); + x = Interval(-oo,oo) for i in range(10): self.assertTrue(x.contains(x.rand())) x = Interval.empty() diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp index cf5214f0c..138df45b3 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp @@ -1,7 +1,7 @@ /** * Codac tests * ---------------------------------------------------------------------------- - * \date 2024 + * \date 2025 * \author Maël Godard * \copyright Copyright 2024 Codac Team * \license GNU Lesser General Public License (LGPL) @@ -9,7 +9,6 @@ #include #include -#include using namespace std; using namespace codac2; @@ -18,6 +17,7 @@ TEST_CASE("Parallelepiped") { Parallelepiped p(Vector({0,2,4}), Matrix({{0.5,0.,0.},{0.,1.,0.},{0.,1.,1.}})); CHECK(p.box() == IntervalVector({{-0.5,0.5},{1.,3.},{2.,6.}})); + CHECK(p.contains(Vector({0.1,2.1,4.1}))); Zonotope z = p.proj({2,1,0}); CHECK(z.z == Vector({4,2,0})); diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py index 196ed764c..888a9c85e 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py @@ -2,7 +2,7 @@ # Codac tests # ---------------------------------------------------------------------------- -# \date 2024 +# \date 2025 # \author Maël Godard # \copyright Copyright 2024 Codac Team # \license GNU Lesser General Public License (LGPL) @@ -18,6 +18,7 @@ def test_parallelepiped(self): p = Parallelepiped(Vector([0,2,4]), Matrix([[0.5,0,0],[0,1,0],[0,1,1]])) self.assertTrue(p.box() == IntervalVector([[-0.5,0.5],[1,3],[2,6]])) + self.assertTrue(p.contains(Vector([0.1,2.1,4.1]))) z = p.proj([2,1,0]) self.assertTrue(z.z == Vector([4,2,0])) diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp new file mode 100644 index 000000000..5205b8122 --- /dev/null +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp @@ -0,0 +1,60 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("Parallelepiped_eval") +{ + ScalarVar x1; + ScalarVar x2; + VectorVar x(2); + + AnalyticFunction f1({x1},{x1,sqr(x1)}); + + AnalyticFunction f2({x1,x2}, {x1, x2, sqr(x1)+sqr(x2)}); + AnalyticFunction f3({x}, {x[0], x[1], sqr(x[0])+sqr(x[1])}); + + auto p1 = f1.parallelepiped_eval(Interval(-0.1,0.1)); + + CHECK(Approx(p1.z,1e-6) == Vector({0,0})); + CHECK(Approx(p1.A,1e-6) == Matrix({{0.12,0},{0,0.02}})); + + auto p = f2.parallelepiped_eval(Interval(-0.1,0.1), Interval(-0.1,0.1)); + + CHECK(Approx(p.z,1e-6) == Vector({0,0,0})); + CHECK(Approx(p.A,1e-6) == Matrix({{0.14,0,0},{0,0.14,0},{0,0,0.04}})); + + + double dx = 0.4; + double x0 = -2.0; + while (x0<2.0) + { + Interval X0(x0, x0+dx); + double y0=-2.0; + while (y0<2.0) + { + Interval Y0(y0, y0+dx); + + 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.A,1e-6) == p3.A); + y0+=dx; + } + x0+=dx; + } + +} + diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py new file mode 100644 index 000000000..d2d232765 --- /dev/null +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py @@ -0,0 +1,56 @@ +#!/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 TestParallelepipedEval(unittest.TestCase): + + def test_parallelepiped_eval(self): + x1 = ScalarVar() + x2 = ScalarVar() + x = VectorVar(2) + + f1 = AnalyticFunction([x1], [x1,sqr(x1)]) + + f2 = AnalyticFunction([x1,x2], [x1, x2, sqr(x1)+sqr(x2)]) + f3 = AnalyticFunction([x], [x[0], x[1], sqr(x[0])+sqr(x[1])]) + + p1 = f1.parallelepiped_eval(Interval(-0.1,0.1)) + + self.assertTrue(Approx(p1.z,1e-6)==Vector([0.0,0.0])) + self.assertTrue(Approx(p1.A,1e-6)==Matrix([[0.12,0.0],[0.0,0.02]])) + + p = f2.parallelepiped_eval(Interval(-0.1,0.1), Interval(-0.1,0.1)) + + self.assertTrue(Approx(p.z,1e-6)==Vector([0,0,0])) + self.assertTrue(Approx(p.A,1e-6)==Matrix([[0.14,0,0],[0,0.14,0],[0,0,0.04]])) + + + dx = 0.4 + x0 = -2.0 + while x0<2.0: + X0 = Interval(x0,x0+dx) + y0 = -2.0 + while y0<2.0: + Y0 = Interval(y0,y0+dx) + + 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.A,1e-6)==p3.A) + y0 += dx + x0 += dx + + +if __name__ == '__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 new file mode 100644 index 000000000..d51015c38 --- /dev/null +++ b/tests/core/peibos/codac2_tests_peibos.cpp @@ -0,0 +1,85 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("Peibos") +{ + + // 2D checks of the PEIBOS algorithm + + VectorVar y_2d(2); + AnalyticFunction f_2d({y_2d},{sqr(y_2d[0])-sqr(y_2d[1])+y_2d[0], + 2.*y_2d[0]*y_2d[1] + y_2d[1]}); + + VectorVar X_2d(1); + AnalyticFunction psi0_2d ({X_2d},{cos(X_2d[0]*PI/4.),sin(X_2d[0]*PI/4.)}); + + OctaSym id_2d ({1,2}); + OctaSym s ({-2,1}); + + auto v_par_2d = PEIBOS(f_2d, psi0_2d, {id_2d,s,s*s,s.invert()}, 0.25, {-0.2,0.}); + + Vector b0 ({-0.5,0.}); + Vector b1 ({0.,1.45}); + Vector b2 ({-1.165,0.}); + + int count_b0 = 0; int count_b1 = 0; int count_b2 = 0; + + for (const auto& p : v_par_2d) + { + if (p.contains(b0)) + count_b0++; + if (p.contains(b1)) + count_b1++; + if (p.contains(b2)) + count_b2++; + } + + CHECK(count_b0 == 0); + CHECK(count_b1 == 1); + CHECK(count_b2 == 2); + + // 3D checks of the PEIBOS algorithm + + VectorVar y_3d(3); + AnalyticFunction f_3d({y_3d},{y_3d[0],y_3d[1],y_3d[2]}); + + VectorVar X_3d(2); + AnalyticFunction psi0_3d ({X_3d},{1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))}); + + OctaSym id_3d ({1,2,3}); + OctaSym s1 ({-2,1,3}); + OctaSym s2 ({3,2,-1}); + + auto v_par_3d = PEIBOS(f_3d, psi0_3d, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 2.); + + 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.})); + + double a = 4.35066; + + CHECK(Approx(v_par_3d[0].A,1e-5) == Matrix({{0.,0.,a},{a+1,0.,0.},{0.,a+1,0.}})); + CHECK(Approx(v_par_3d[1].A,1e-5) == Matrix({{-(a+1),0.,0.},{0.,0.,a},{0.,a+1,0.}})); + CHECK(Approx(v_par_3d[2].A,1e-5) == Matrix({{0.,0.,a},{-(a+1),0.,0.},{0.,a+1,0.}})); + CHECK(Approx(v_par_3d[3].A,1e-5) == Matrix({{a+1,0.,0.},{0.,0.,a},{0.,a+1,0.}})); + CHECK(Approx(v_par_3d[4].A,1e-5) == Matrix({{0.,a+1,0.},{a+1,0.,0.},{0.,0.,a}})); + CHECK(Approx(v_par_3d[5].A,1e-5) == Matrix({{0.,-(a+1),0.},{a+1,0.,0.},{0.,0.,a}})); +} \ No newline at end of file diff --git a/tests/core/peibos/codac2_tests_peibos.py b/tests/core/peibos/codac2_tests_peibos.py new file mode 100644 index 000000000..3e6f6587a --- /dev/null +++ b/tests/core/peibos/codac2_tests_peibos.py @@ -0,0 +1,83 @@ +#!/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 TestPeibos(unittest.TestCase): + + def test_peibos(self): + + # 2D checks of the PEIBOS algorithm + + y_2d = VectorVar(2) + f_2d = AnalyticFunction([y_2d],[sqr(y_2d[0])-sqr(y_2d[1])+y_2d[0],2.*y_2d[0]*y_2d[1] + y_2d[1]]) + + X_2d = VectorVar(1) + psi0_2d = AnalyticFunction([X_2d],[cos(X_2d[0]*PI/4.),sin(X_2d[0]*PI/4.)]) + + id_2d = OctaSym([1, 2]) + s = OctaSym([-2, 1]) + + v_par_2d = PEIBOS(f_2d,psi0_2d,[id_2d,s,s*s,s.invert()],0.25,[-0.2,0.]) + + b0 = Vector([-0.5,0.0]) + b1 = Vector([0.0,1.45]) + b2 = Vector([-1.165,0.0]) + + count_b0, count_b1, count_b2 = 0, 0, 0 + + for p in v_par_2d: + if p.contains(b0): + count_b0 += 1 + if p.contains(b1): + count_b1 += 1 + if p.contains(b2): + count_b2 += 1 + + self.assertTrue(count_b0 == 0) + self.assertTrue(count_b1 == 1) + self.assertTrue(count_b2 == 2) + + # 3D checks of the PEIBOS algorithm + + y_3d = VectorVar(3) + f_3d = AnalyticFunction([y_3d],[y_3d[0],y_3d[1],y_3d[2]]) + + X_3d = VectorVar(2) + psi0_3d = AnalyticFunction([X_3d],[1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))]) + + id_3d = OctaSym([1, 2, 3]) + s1 = OctaSym([-2, 1, 3]) + s2 = OctaSym([3, 2, -1]) + + v_par_3d = PEIBOS(f_3d,psi0_3d,[id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()],2.0) + + 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.])) + + a = 4.35066 + + self.assertTrue(Approx(v_par_3d[0].A,1e-5) == Matrix([[0.,0.,a],[a+1,0.,0.],[0.,a+1,0.]])) + self.assertTrue(Approx(v_par_3d[1].A,1e-5) == Matrix([[-(a+1),0.,0.],[0.,0.,a],[0.,a+1,0.]])) + self.assertTrue(Approx(v_par_3d[2].A,1e-5) == Matrix([[0.,0.,a],[-(a+1),0.,0.],[0.,a+1,0.]])) + self.assertTrue(Approx(v_par_3d[3].A,1e-5) == Matrix([[a+1,0.,0.],[0.,0.,a],[0.,a+1,0.]])) + self.assertTrue(Approx(v_par_3d[4].A,1e-5) == Matrix([[0.,a+1,0.],[a+1,0.,0.],[0.,0.,a]])) + self.assertTrue(Approx(v_par_3d[5].A,1e-5) == Matrix([[0.,-(a+1),0.],[a+1,0.,0.],[0.,0.,a]])) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file