From 92b49ce39684058c7e67746d0aa7955c64604b44 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 30 Sep 2025 09:25:35 +0200 Subject: [PATCH 01/24] [doc] minor update --- doc/manual/manual/visualization/colors.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 ------ From 7ff5ae4505e9254e83d3af4cf8270a09bf2eeb04 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 30 Sep 2025 12:16:47 +0200 Subject: [PATCH 02/24] [peibos] first clean implementation (C++ only) --- examples/07_centered_2D/main.cpp | 9 +- examples/10_peibos/CMakeLists.txt | 34 ++++++ examples/10_peibos/main.cpp | 108 ++++++++++++++++++ src/core/CMakeLists.txt | 4 + src/core/peibos/codac2_peibos.cpp | 176 ++++++++++++++++++++++++++++++ src/core/peibos/codac2_peibos.h | 30 +++++ 6 files changed, 360 insertions(+), 1 deletion(-) create mode 100644 examples/10_peibos/CMakeLists.txt create mode 100644 examples/10_peibos/main.cpp create mode 100644 src/core/peibos/codac2_peibos.cpp create mode 100644 src/core/peibos/codac2_peibos.h diff --git a/examples/07_centered_2D/main.cpp b/examples/07_centered_2D/main.cpp index af464ecd1..166dc3ce6 100644 --- a/examples/07_centered_2D/main.cpp +++ b/examples/07_centered_2D/main.cpp @@ -17,10 +17,14 @@ int main(){ // we use Fermat's spiral AnalyticFunction f1 ({t},{a*sqrt(t)*cos(t),a*sqrt(t)*sin(t)}); + VectorVar t_vec(1); + AnalyticFunction f2 ({t_vec},{a*sqrt(t_vec[0])*cos(t_vec[0]),a*sqrt(t_vec[0])*sin(t_vec[0])}); + double dt=0.2; double time=0.0; while (time<100.0) { Interval T(time,time+dt); + IntervalVector Tvec ({T}); // using eval_ would be easier, but it's not allowed :( IntervalVector df = f1.diff(T); if (df.is_empty() || df.is_unbounded()) { @@ -35,7 +39,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 = parallelepiped_inclusion(f2,Tvec); + fig4.draw_parallelepiped(p,{{Color::black(),Color::green(0.1)},"parallelepipeds"}); } time = time+dt; } diff --git a/examples/10_peibos/CMakeLists.txt b/examples/10_peibos/CMakeLists.txt new file mode 100644 index 000000000..cf9ef35ce --- /dev/null +++ b/examples/10_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/10_peibos/main.cpp b/examples/10_peibos/main.cpp new file mode 100644 index 000000000..47c63ab0f --- /dev/null +++ b/examples/10_peibos/main.cpp @@ -0,0 +1,108 @@ +#include +#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)); + Zonotope z = p.proj({0,1}); + figure_3d_proj.draw_zonotope(z , {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} }); + 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]))}); + 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}); + OctaSym s3_nd ({1,-2,-3}); + + + 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)); + Zonotope z_xy = p.proj({0,1}); + Zonotope z_zy = p.proj({2,1}); + figure_2d_nd_xy.draw_zonotope(z_xy, {Color::black(),Color::green(0.2)}); + figure_2d_nd_zy.draw_zonotope(z_zy, {Color::black(),Color::green(0.2)}); + } +} \ No newline at end of file diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 97c741682..fec66ad14 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -190,6 +190,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 @@ -282,6 +285,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/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp new file mode 100644 index 000000000..07c62399a --- /dev/null +++ b/src/core/peibos/codac2_peibos.cpp @@ -0,0 +1,176 @@ +/** + * 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" + +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 IntervalMatrix& JJg, const IntervalMatrix& JJg_punc, const IntervalVector& X) + { + auto xc = X.mid(); + + auto dX=X-xc; + + auto E = (JJg - JJg_punc)*dX; + Interval N(0.); + + for (int i = 0; i < E.rows(); i++) + N += sqr(E[i]); + + return sqrt(N).ub(); + } + + double error(const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) + { + auto xc = X.mid(); + + IntervalMatrix JJg_punc=JJf_punc*(symmetry.permutation_matrix().template cast())*psi_0.diff(xc); + + IntervalMatrix JJg=JJf*(symmetry.permutation_matrix().template cast())*psi_0.diff(X); + + return error(JJg, JJg_punc, X); + } + + Matrix inflate_flat_parallelepiped(const Matrix& A, const Vector& e_vec, double rho) + { + Index m = A.cols(); + Index n = A.rows(); + + Eigen::FullPivLU lu_decomp(A.transpose()); + Eigen::MatrixXd N = lu_decomp.kernel(); + + Matrix A_tild (n,n); + A_tild << A, N; + + Matrix Q = (A_tild.transpose() * A_tild).inverse(); + + // Matrix mult (n, n); + Matrix mult = Matrix::Zero(n,n); + for (int i = 0; i < n; i++) + mult(i,i) = rho*std::sqrt(Q(i,i)); + + for (int i = 0; i < m; i++) + mult(i,i) += e_vec(i); + + return A_tild*mult; + } + + Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X) + { + assert_release (f.input_size() < f.output_size() && "input size must be strictly lower than output size"); + assert_release (X.size() == f.input_size() && "dimension of X must match input size of f"); + + VectorVar psi_0(f.input_size()); + + // Maximum error computation + double rho = error(f.diff(X), f.diff(X.mid()).mid(), X); + + Vector x_rad (X.size()); + for (int i = 0; i < X.size(); i++) + x_rad(i) = X(i).rad(); + + Matrix A = f.diff(X.mid()).mid(); + + // Inflation of the parallelepiped + auto A_inf = inflate_flat_parallelepiped(A, x_rad, rho); + + return Parallelepiped(f.eval(X.mid()).mid(), A_inf); + } + + Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X, double true_eps) + { + // Maximum error computation + double rho = error( JJf, JJf_punc, psi_0, symmetry, X); + + auto Jz = (JJf_punc * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X.mid())).mid(); + + // Inflation of the parallelepiped + + auto A = inflate_flat_parallelepiped(Jz, 0.5*true_eps*Vector::ones(X.size()), rho); + + return Parallelepiped(z, A); + } + + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose) + { + return PEIBOS(f, psi_0, symmetries, epsilon, Vector::Zero(psi_0.output_size()), verbose); + } + + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, 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 (symmetries.size() > 0 && (int) symmetries[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(Interval(-1.,1.)*IntervalVector::Ones(m), epsilon, boxes); + + for (const auto& symmetry : symmetries) + { + for (const auto& X : boxes) + { + IntervalVector Y = symmetry(psi_0.eval(X)) + offset; + + auto JJf=f.diff(Y); + + auto xc = X.mid(); + auto yc = (symmetry(psi_0.eval(xc)) + offset).mid(); + + auto JJf_punc=f.diff(yc).mid(); + + // Center of the parallelepiped + auto z = f.eval(yc).mid(); + + auto p = parallelepiped_inclusion(z, JJf, JJf_punc, psi_0, symmetry, X, true_eps); + + output.push_back(p); + + } + } + + if (verbose) + { + printf("\nPEIBOS statistics:\n"); + printf("------------------\n"); + printf("Number of symmetries: %ld\n", symmetries.size()); + 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..87a4c7c09 --- /dev/null +++ b/src/core/peibos/codac2_peibos.h @@ -0,0 +1,30 @@ +/** + * \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" + +namespace codac2 +{ + double split (const IntervalVector& X, double eps, vector& boxes); + + double error(const IntervalMatrix& JJg, const IntervalMatrix& JJg_punc, const IntervalVector& X); + double error(const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); + + Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); + + Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X); + Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X, double true_eps); + + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose = false); + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false); +} \ No newline at end of file From 62d09e4a648079a13e777945d7f704b0dbaf1d81 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 30 Sep 2025 14:05:30 +0200 Subject: [PATCH 03/24] [peibos] simplifying code --- examples/10_peibos/main.cpp | 3 +++ src/core/peibos/codac2_peibos.cpp | 18 ++++++------------ src/core/peibos/codac2_peibos.h | 1 - 3 files changed, 9 insertions(+), 13 deletions(-) diff --git a/examples/10_peibos/main.cpp b/examples/10_peibos/main.cpp index 47c63ab0f..7c08f8717 100644 --- a/examples/10_peibos/main.cpp +++ b/examples/10_peibos/main.cpp @@ -68,7 +68,10 @@ int main() 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); diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 07c62399a..1960e7d35 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -45,17 +45,6 @@ namespace codac2 return sqrt(N).ub(); } - double error(const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) - { - auto xc = X.mid(); - - IntervalMatrix JJg_punc=JJf_punc*(symmetry.permutation_matrix().template cast())*psi_0.diff(xc); - - IntervalMatrix JJg=JJf*(symmetry.permutation_matrix().template cast())*psi_0.diff(X); - - return error(JJg, JJg_punc, X); - } - Matrix inflate_flat_parallelepiped(const Matrix& A, const Vector& e_vec, double rho) { Index m = A.cols(); @@ -104,8 +93,12 @@ namespace codac2 Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X, double true_eps) { + // Computation of the Jacobian of g = f o symmetry(psi_0) + IntervalMatrix JJg=JJf*(symmetry.permutation_matrix().template cast())*psi_0.diff(X); + IntervalMatrix JJg_punc=JJf_punc*(symmetry.permutation_matrix().template cast())*psi_0.diff(X.mid()); + // Maximum error computation - double rho = error( JJf, JJf_punc, psi_0, symmetry, X); + double rho = error(JJg, JJg_punc, X); auto Jz = (JJf_punc * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X.mid())).mid(); @@ -139,6 +132,7 @@ namespace codac2 for (const auto& symmetry : symmetries) { + // Simon : would it be possible to create a function g = f o symmetry(psi_0) ? for (const auto& X : boxes) { IntervalVector Y = symmetry(psi_0.eval(X)) + offset; diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index 87a4c7c09..c2cb0c6be 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -18,7 +18,6 @@ namespace codac2 double split (const IntervalVector& X, double eps, vector& boxes); double error(const IntervalMatrix& JJg, const IntervalMatrix& JJg_punc, const IntervalVector& X); - double error(const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); From d56d986dbf244be08faffb89fa8c91856cd3edd1 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 30 Sep 2025 14:28:55 +0200 Subject: [PATCH 04/24] [peibos] minor modification of example --- examples/10_peibos/main.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/examples/10_peibos/main.cpp b/examples/10_peibos/main.cpp index 7c08f8717..9c3944de9 100644 --- a/examples/10_peibos/main.cpp +++ b/examples/10_peibos/main.cpp @@ -55,8 +55,7 @@ int main() for (const auto& p : v_par_3d) { figure3d.draw_parallelepiped(p, Color::green(0.5)); - Zonotope z = p.proj({0,1}); - figure_3d_proj.draw_zonotope(z , {Color::black(),Color::green(0.2)}); + figure_3d_proj.draw_zonotope(p.proj({0,1}) , {Color::black(),Color::green(0.2)}); } // nD example of the PEIBOS algorithm @@ -103,9 +102,7 @@ int main() for (const auto& p : v_par_nd) { figure_3d_nd.draw_parallelepiped(p, Color::green(0.5)); - Zonotope z_xy = p.proj({0,1}); - Zonotope z_zy = p.proj({2,1}); - figure_2d_nd_xy.draw_zonotope(z_xy, {Color::black(),Color::green(0.2)}); - figure_2d_nd_zy.draw_zonotope(z_zy, {Color::black(),Color::green(0.2)}); + 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 From 30efde23eefbab61c68140cf5a2ea7d6951eb827 Mon Sep 17 00:00:00 2001 From: godardma Date: Wed, 1 Oct 2025 14:06:51 +0200 Subject: [PATCH 05/24] [graphics] added colormap rainbow_05 --- src/graphics/styles/codac2_ColorMap.h | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/graphics/styles/codac2_ColorMap.h b/src/graphics/styles/codac2_ColorMap.h index 51102f0bc..62371175b 100644 --- a/src/graphics/styles/codac2_ColorMap.h +++ b/src/graphics/styles/codac2_ColorMap.h @@ -166,5 +166,21 @@ namespace codac2 return cmap; } + /** + * \brief Rainbow color map with 50% saturation and 50% transparency + * + * Rainbow color map with 50% saturation and 50% transparency + */ + static ColorMap rainbow_05() + { + ColorMap cmap( Model::HSV ); + int i = 0; + for(int h = 300 ; h > 0 ; h-=10) + { + cmap[i]=Color({(float)h,50.,100.,50.},Model::HSV); + i++; + } + return cmap; + } }; } \ No newline at end of file From 150e9929ba90b81b0a69a65511520b81d4d01f75 Mon Sep 17 00:00:00 2001 From: godardma Date: Thu, 2 Oct 2025 12:40:18 +0200 Subject: [PATCH 06/24] [peibos] code simplification (thanks Simon) --- examples/10_peibos/main.cpp | 1 - src/core/peibos/codac2_peibos.cpp | 24 ++++-------------------- src/core/peibos/codac2_peibos.h | 2 ++ 3 files changed, 6 insertions(+), 21 deletions(-) diff --git a/examples/10_peibos/main.cpp b/examples/10_peibos/main.cpp index 9c3944de9..524844c45 100644 --- a/examples/10_peibos/main.cpp +++ b/examples/10_peibos/main.cpp @@ -1,5 +1,4 @@ #include -#include using namespace std; using namespace codac2; diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 1960e7d35..24ea881df 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -103,7 +103,6 @@ namespace codac2 auto Jz = (JJf_punc * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X.mid())).mid(); // Inflation of the parallelepiped - auto A = inflate_flat_parallelepiped(Jz, 0.5*true_eps*Vector::ones(X.size()), rho); return Parallelepiped(z, A); @@ -132,26 +131,11 @@ namespace codac2 for (const auto& symmetry : symmetries) { - // Simon : would it be possible to create a function g = f o symmetry(psi_0) ? - for (const auto& X : boxes) - { - IntervalVector Y = symmetry(psi_0.eval(X)) + offset; - - auto JJf=f.diff(Y); - - auto xc = X.mid(); - auto yc = (symmetry(psi_0.eval(xc)) + offset).mid(); - - auto JJf_punc=f.diff(yc).mid(); - - // Center of the parallelepiped - auto z = f.eval(yc).mid(); - - auto p = parallelepiped_inclusion(z, JJf, JJf_punc, psi_0, symmetry, X, true_eps); - - output.push_back(p); + VectorVar x(m); + AnalyticFunction g_i ({x}, f(symmetry(psi_0(x))+offset)); - } + for (const auto& X : boxes) + output.push_back( parallelepiped_inclusion(g_i, X)); } if (verbose) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index c2cb0c6be..a4c2f47b0 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -12,6 +12,8 @@ #include "codac2_Parallelepiped.h" #include "codac2_AnalyticFunction.h" #include "codac2_OctaSym.h" +#include + namespace codac2 { From 5c01aaf97e0af049f480f8914be63acb1ee20357 Mon Sep 17 00:00:00 2001 From: godardma Date: Thu, 2 Oct 2025 12:41:18 +0200 Subject: [PATCH 07/24] [peibos] minor modification --- src/core/peibos/codac2_peibos.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 24ea881df..cd34a4013 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -148,7 +148,5 @@ namespace codac2 } return output; - } - } \ No newline at end of file From 5c291b6c3fb2b045da0d3155ad7b09104e9c63c7 Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 3 Oct 2025 10:11:24 +0200 Subject: [PATCH 08/24] [peibos] python bindings --- examples/10_peibos/main.cpp | 2 - examples/10_peibos/main.py | 88 +++++++++++++++++++++ python/src/core/CMakeLists.txt | 3 + python/src/core/codac2_py_core.cpp | 6 ++ python/src/core/peibos/codac2_py_peibos.cpp | 60 ++++++++++++++ src/core/peibos/codac2_peibos.cpp | 10 ++- src/core/peibos/codac2_peibos.h | 2 +- 7 files changed, 165 insertions(+), 6 deletions(-) create mode 100644 examples/10_peibos/main.py create mode 100644 python/src/core/peibos/codac2_py_peibos.cpp diff --git a/examples/10_peibos/main.cpp b/examples/10_peibos/main.cpp index 524844c45..fbbf1ec51 100644 --- a/examples/10_peibos/main.cpp +++ b/examples/10_peibos/main.cpp @@ -78,8 +78,6 @@ int main() OctaSym id_nd ({1,2,3}); OctaSym s1_nd ({-2,1,3}); OctaSym s2_nd ({3,2,-1}); - OctaSym s3_nd ({1,-2,-3}); - Figure3D figure_3d_nd ("Cube on Sphere"); figure_3d_nd.draw_axes(0.5); diff --git a/examples/10_peibos/main.py b/examples/10_peibos/main.py new file mode 100644 index 000000000..cc871bbe9 --- /dev/null +++ b/examples/10_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.]) + + 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) + + 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) + + 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/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 c0f307992..2a09742b0 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -119,6 +119,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); @@ -272,6 +275,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/peibos/codac2_py_peibos.cpp b/python/src/core/peibos/codac2_py_peibos.cpp new file mode 100644 index 000000000..dd874d740 --- /dev/null +++ b/python/src/core/peibos/codac2_py_peibos.cpp @@ -0,0 +1,60 @@ +/** + * 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 py::object& f, const IntervalVector& X) + { + return parallelepiped_inclusion(cast>(f), X); + }, + PARALLELEPIPED_PARALLELEPIPED_INCLUSION_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_INTERVALVECTOR_REF, + "f"_a, "X"_a); + + m.def("parallelepiped_inclusion", + [](const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const py::object& psi_0, const OctaSym& symmetry, const IntervalVector& X) + { + return parallelepiped_inclusion(z, JJf, JJf_punc, cast>(psi_0), symmetry, X); + }, + PARALLELEPIPED_PARALLELEPIPED_INCLUSION_CONST_VECTOR_REF_CONST_INTERVALMATRIX_REF_CONST_INTERVALMATRIX_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_OCTASYM_REF_CONST_INTERVALVECTOR_REF, + "z"_a, "JJf"_a, "JJf_punc"_a, "psi_0"_a, "symmetry"_a, "X"_a); + + + m.def("PEIBOS", + [](const py::object& f, const py::object& psi_0, const vector& symmetries, double epsilon, bool verbose = false) + { + return PEIBOS(cast>(f), cast>(psi_0), symmetries, epsilon, verbose); + }, + VECTOR_PARALLELEPIPED_PEIBOS_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_VECTOR_OCTASYM_REF_DOUBLE_BOOL, + "f"_a, "psi_0"_a, "symmetries"_a, "epsilon"_a, "verbose"_a = false); + + m.def("PEIBOS", + [](const py::object& f, const py::object& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false) + { + return PEIBOS(cast>(f), cast>(psi_0), symmetries, 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, "symmetries"_a, "epsilon"_a, "offset"_a, "verbose"_a = false); +} \ No newline at end of file diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index cd34a4013..6a04081eb 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -91,7 +91,7 @@ namespace codac2 return Parallelepiped(f.eval(X.mid()).mid(), A_inf); } - Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X, double true_eps) + Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) { // Computation of the Jacobian of g = f o symmetry(psi_0) IntervalMatrix JJg=JJf*(symmetry.permutation_matrix().template cast())*psi_0.diff(X); @@ -102,8 +102,12 @@ namespace codac2 auto Jz = (JJf_punc * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X.mid())).mid(); + Vector x_rad (X.size()); + for (int i = 0; i < X.size(); i++) + x_rad(i) = X(i).rad(); + // Inflation of the parallelepiped - auto A = inflate_flat_parallelepiped(Jz, 0.5*true_eps*Vector::ones(X.size()), rho); + auto A = inflate_flat_parallelepiped(Jz, x_rad, rho); return Parallelepiped(z, A); } @@ -135,7 +139,7 @@ namespace codac2 AnalyticFunction g_i ({x}, f(symmetry(psi_0(x))+offset)); for (const auto& X : boxes) - output.push_back( parallelepiped_inclusion(g_i, X)); + output.push_back(parallelepiped_inclusion(g_i, X)); } if (verbose) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index a4c2f47b0..9dfda8922 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -24,7 +24,7 @@ namespace codac2 Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X); - Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X, double true_eps); + Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose = false); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false); From 469bcca36c83f1bdc5f9b720888619c9e6ac9263 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 7 Oct 2025 16:00:25 +0200 Subject: [PATCH 09/24] [peibos] minor modifications following PR comments (Thanks Simon) --- src/core/peibos/codac2_peibos.cpp | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 6a04081eb..b43ee41cb 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -74,14 +74,10 @@ namespace codac2 assert_release (f.input_size() < f.output_size() && "input size must be strictly lower than output size"); assert_release (X.size() == f.input_size() && "dimension of X must match input size of f"); - VectorVar psi_0(f.input_size()); - // Maximum error computation double rho = error(f.diff(X), f.diff(X.mid()).mid(), X); - Vector x_rad (X.size()); - for (int i = 0; i < X.size(); i++) - x_rad(i) = X(i).rad(); + Vector x_rad = X.rad(); Matrix A = f.diff(X.mid()).mid(); @@ -102,9 +98,7 @@ namespace codac2 auto Jz = (JJf_punc * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X.mid())).mid(); - Vector x_rad (X.size()); - for (int i = 0; i < X.size(); i++) - x_rad(i) = X(i).rad(); + Vector x_rad = X.rad(); // Inflation of the parallelepiped auto A = inflate_flat_parallelepiped(Jz, x_rad, rho); @@ -131,7 +125,7 @@ namespace codac2 vector output; vector boxes; - double true_eps = split(Interval(-1.,1.)*IntervalVector::Ones(m), epsilon, boxes); + double true_eps = split(IntervalVector::constant(m,{-1,1}), epsilon, boxes); for (const auto& symmetry : symmetries) { From 1c45b83807c6cb7f17ace8674792259b4182acd6 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 7 Oct 2025 17:30:19 +0200 Subject: [PATCH 10/24] [peibos] modification to preserve the guarantee --- python/src/core/peibos/codac2_py_peibos.cpp | 4 +-- src/core/peibos/codac2_peibos.cpp | 33 +++++++++++---------- src/core/peibos/codac2_peibos.h | 4 +-- 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/python/src/core/peibos/codac2_py_peibos.cpp b/python/src/core/peibos/codac2_py_peibos.cpp index dd874d740..946285a41 100644 --- a/python/src/core/peibos/codac2_py_peibos.cpp +++ b/python/src/core/peibos/codac2_py_peibos.cpp @@ -34,11 +34,11 @@ void export_peibos(py::module& m) "f"_a, "X"_a); m.def("parallelepiped_inclusion", - [](const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const py::object& psi_0, const OctaSym& symmetry, const IntervalVector& X) + [](const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& JJf_punc, const py::object& psi_0, const OctaSym& symmetry, const IntervalVector& X) { return parallelepiped_inclusion(z, JJf, JJf_punc, cast>(psi_0), symmetry, X); }, - PARALLELEPIPED_PARALLELEPIPED_INCLUSION_CONST_VECTOR_REF_CONST_INTERVALMATRIX_REF_CONST_INTERVALMATRIX_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_OCTASYM_REF_CONST_INTERVALVECTOR_REF, + PARALLELEPIPED_PARALLELEPIPED_INCLUSION_CONST_INTERVALVECTOR_REF_CONST_INTERVALMATRIX_REF_CONST_MATRIX_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_OCTASYM_REF_CONST_INTERVALVECTOR_REF, "z"_a, "JJf"_a, "JJf_punc"_a, "psi_0"_a, "symmetry"_a, "X"_a); diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index b43ee41cb..109c21630 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -30,7 +30,7 @@ namespace codac2 } } - double error(const IntervalMatrix& JJg, const IntervalMatrix& JJg_punc, const IntervalVector& X) + double error(const IntervalMatrix& JJg, const Matrix& JJg_punc, const IntervalVector& X) { auto xc = X.mid(); @@ -74,36 +74,39 @@ namespace codac2 assert_release (f.input_size() < f.output_size() && "input size must be strictly lower than output size"); assert_release (X.size() == f.input_size() && "dimension of X must match input size of f"); - // Maximum error computation - double rho = error(f.diff(X), f.diff(X.mid()).mid(), X); - - Vector x_rad = X.rad(); + auto z = f.eval(X.mid()); Matrix A = f.diff(X.mid()).mid(); + // Maximum error computation + double rho = error(f.diff(X), A, X); + + // We need to add the radius of z to rho to account for the fact that we have a (small) box enclosing f(x_bar) and not f(x_bar) itself + rho += z.rad().maxCoeff(); + // Inflation of the parallelepiped - auto A_inf = inflate_flat_parallelepiped(A, x_rad, rho); + auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); - return Parallelepiped(f.eval(X.mid()).mid(), A_inf); + return Parallelepiped(z.mid(), A_inf); } - Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) + Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) { // Computation of the Jacobian of g = f o symmetry(psi_0) - IntervalMatrix JJg=JJf*(symmetry.permutation_matrix().template cast())*psi_0.diff(X); - IntervalMatrix JJg_punc=JJf_punc*(symmetry.permutation_matrix().template cast())*psi_0.diff(X.mid()); + IntervalMatrix JJg = JJf * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X); + + Matrix JJg_punc = (JJf_punc * symmetry.permutation_matrix() * (psi_0.diff(X.mid()).mid())); // Maximum error computation double rho = error(JJg, JJg_punc, X); - auto Jz = (JJf_punc * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X.mid())).mid(); - - Vector x_rad = X.rad(); + // We need to add the radius of z to rho to account for the fact that we have a (small) box enclosing z and not z itself + rho += z.rad().maxCoeff(); // Inflation of the parallelepiped - auto A = inflate_flat_parallelepiped(Jz, x_rad, rho); + auto A = inflate_flat_parallelepiped(JJg_punc, X.rad(), rho); - return Parallelepiped(z, A); + return Parallelepiped(z.mid(), A); } vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index 9dfda8922..a5e337812 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -19,12 +19,12 @@ namespace codac2 { double split (const IntervalVector& X, double eps, vector& boxes); - double error(const IntervalMatrix& JJg, const IntervalMatrix& JJg_punc, const IntervalVector& X); + double error(const IntervalMatrix& JJg, const Matrix& JJg_punc, const IntervalVector& X); Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X); - Parallelepiped parallelepiped_inclusion(const Vector& z, const IntervalMatrix& JJf, const IntervalMatrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); + Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose = false); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false); From bd4901d08c3bc6f7e752690d785a2a7cf30686f7 Mon Sep 17 00:00:00 2001 From: godardma Date: Tue, 7 Oct 2025 21:04:47 +0200 Subject: [PATCH 11/24] [peibos] modification to preserve the guarantee --- src/core/peibos/codac2_peibos.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 109c21630..4fa9b886f 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -37,12 +37,8 @@ namespace codac2 auto dX=X-xc; auto E = (JJg - JJg_punc)*dX; - Interval N(0.); - for (int i = 0; i < E.rows(); i++) - N += sqr(E[i]); - - return sqrt(N).ub(); + return E.norm().ub(); } Matrix inflate_flat_parallelepiped(const Matrix& A, const Vector& e_vec, double rho) @@ -82,7 +78,7 @@ namespace codac2 double rho = error(f.diff(X), A, X); // We need to add the radius of z to rho to account for the fact that we have a (small) box enclosing f(x_bar) and not f(x_bar) itself - rho += z.rad().maxCoeff(); + rho += z.norm().diam(); // Inflation of the parallelepiped auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); @@ -101,7 +97,7 @@ namespace codac2 double rho = error(JJg, JJg_punc, X); // We need to add the radius of z to rho to account for the fact that we have a (small) box enclosing z and not z itself - rho += z.rad().maxCoeff(); + rho += z.norm().diam(); // Inflation of the parallelepiped auto A = inflate_flat_parallelepiped(JJg_punc, X.rad(), rho); From a6b984a081b8631b31c359f4b1cbbb26ad3fa7d2 Mon Sep 17 00:00:00 2001 From: godardma Date: Wed, 8 Oct 2025 14:29:19 +0200 Subject: [PATCH 12/24] [peibos] modification to preserve the guarantee (last one I hope) --- python/src/core/peibos/codac2_py_peibos.cpp | 6 ++-- src/core/peibos/codac2_peibos.cpp | 31 +++++++++------------ src/core/peibos/codac2_peibos.h | 4 +-- 3 files changed, 18 insertions(+), 23 deletions(-) diff --git a/python/src/core/peibos/codac2_py_peibos.cpp b/python/src/core/peibos/codac2_py_peibos.cpp index 946285a41..f0acc8580 100644 --- a/python/src/core/peibos/codac2_py_peibos.cpp +++ b/python/src/core/peibos/codac2_py_peibos.cpp @@ -34,12 +34,12 @@ void export_peibos(py::module& m) "f"_a, "X"_a); m.def("parallelepiped_inclusion", - [](const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& JJf_punc, const py::object& psi_0, const OctaSym& symmetry, const IntervalVector& X) + [](const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_punc, const py::object& psi_0, const OctaSym& symmetry, const IntervalVector& X) { - return parallelepiped_inclusion(z, JJf, JJf_punc, cast>(psi_0), symmetry, X); + return parallelepiped_inclusion(z, JJf, Jf_punc, cast>(psi_0), symmetry, 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, "JJf"_a, "JJf_punc"_a, "psi_0"_a, "symmetry"_a, "X"_a); + "z"_a, "JJf"_a, "Jf_punc"_a, "psi_0"_a, "symmetry"_a, "X"_a); m.def("PEIBOS", diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 4fa9b886f..d4ee09049 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -30,13 +30,15 @@ namespace codac2 } } - double error(const IntervalMatrix& JJg, const Matrix& JJg_punc, const IntervalVector& X) + double error(const IntervalVector& z, const IntervalMatrix& JJg, const Matrix& B, const IntervalVector& X) { auto xc = X.mid(); - auto dX=X-xc; + auto dX = X-xc; - auto E = (JJg - JJg_punc)*dX; + auto a = z.mid(); + + auto E = (z - a) + (JJg - B)*dX; return E.norm().ub(); } @@ -54,7 +56,6 @@ namespace codac2 Matrix Q = (A_tild.transpose() * A_tild).inverse(); - // Matrix mult (n, n); Matrix mult = Matrix::Zero(n,n); for (int i = 0; i < n; i++) mult(i,i) = rho*std::sqrt(Q(i,i)); @@ -72,35 +73,29 @@ namespace codac2 auto z = f.eval(X.mid()); - Matrix A = f.diff(X.mid()).mid(); + Matrix B = f.diff(X.mid()).mid(); // Maximum error computation - double rho = error(f.diff(X), A, X); - - // We need to add the radius of z to rho to account for the fact that we have a (small) box enclosing f(x_bar) and not f(x_bar) itself - rho += z.norm().diam(); + double rho = error(z, f.diff(X), B, X); // Inflation of the parallelepiped - auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); + auto A = inflate_flat_parallelepiped(B, X.rad(), rho); - return Parallelepiped(z.mid(), A_inf); + return Parallelepiped(z.mid(), A); } - Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) + Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) { // Computation of the Jacobian of g = f o symmetry(psi_0) IntervalMatrix JJg = JJf * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X); - Matrix JJg_punc = (JJf_punc * symmetry.permutation_matrix() * (psi_0.diff(X.mid()).mid())); + Matrix B = (Jf_punc * symmetry.permutation_matrix() * (psi_0.diff(X.mid()).mid())); // Maximum error computation - double rho = error(JJg, JJg_punc, X); - - // We need to add the radius of z to rho to account for the fact that we have a (small) box enclosing z and not z itself - rho += z.norm().diam(); + double rho = error(z, JJg, B, X); // Inflation of the parallelepiped - auto A = inflate_flat_parallelepiped(JJg_punc, X.rad(), rho); + auto A = inflate_flat_parallelepiped(B, X.rad(), rho); return Parallelepiped(z.mid(), A); } diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index a5e337812..0123ce6b0 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -19,12 +19,12 @@ namespace codac2 { double split (const IntervalVector& X, double eps, vector& boxes); - double error(const IntervalMatrix& JJg, const Matrix& JJg_punc, const IntervalVector& X); + double error(const IntervalVector& z, const IntervalMatrix& JJg, const Matrix& B, const IntervalVector& X); Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X); - Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& JJf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); + Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose = false); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false); From 6ca22003e75be24825788d245c78ac0b6153bae3 Mon Sep 17 00:00:00 2001 From: godardma Date: Wed, 8 Oct 2025 16:21:57 +0200 Subject: [PATCH 13/24] [peibos] added doc --- src/core/peibos/codac2_peibos.cpp | 9 +++--- src/core/peibos/codac2_peibos.h | 49 ++++++++++++++++++++++++++++++- 2 files changed, 53 insertions(+), 5 deletions(-) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index d4ee09049..e4488bb0d 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -84,12 +84,13 @@ namespace codac2 return Parallelepiped(z.mid(), A); } - Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X) + Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X) { - // Computation of the Jacobian of g = f o symmetry(psi_0) - IntervalMatrix JJg = JJf * (symmetry.permutation_matrix().template cast()) * psi_0.diff(X); + // Computation of the Jacobian of g = f o sigma(psi_0) + IntervalMatrix JJg = JJf * (sigma.permutation_matrix().template cast()) * psi_0.diff(X); - Matrix B = (Jf_punc * symmetry.permutation_matrix() * (psi_0.diff(X.mid()).mid())); + // B is an approximation of the Jacobian of g at the center of X + Matrix B = (Jf_tild * sigma.permutation_matrix() * (psi_0.diff(X.mid()).mid())); // Maximum error computation double rho = error(z, JJg, B, X); diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index 0123ce6b0..0a629326e 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -17,14 +17,61 @@ 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{a}+\mathbf{B}\cdot(\left[\mathbf{x}\right]-\bar{\mathbf{x}})\f$ + * + * \param z The box enclosing \f$\mathbf{f}(\bar{\mathbf{x}})\f$. \f$\mathbf{a}\f$ is its center + * \param JJg The interval Jacobian matrix containing \f$\frac{d\mathbf{f}}{d\mathbf{x}}([\mathbf{x}])\f$ + * \param B Any Matrix with the same dimension as JJg. 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& z, const IntervalMatrix& JJg, const Matrix& B, 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 Compute a parallelepiped enclosing of \f$\mathbf{f}([\mathbf{x}])\f$ when \f$\mathbf{f}:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ is an analytic function with \f$m < n\f$. + * + * \param f The analytic function \f$\mathbf{f}\f$ + * \param X The box \f$[\mathbf{x}]\f$ + * + * \return A Parallelepiped enclosing \f$\mathbf{f}([\mathbf{x}])\f$ + */ Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X); - Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_punc, const AnalyticFunction& psi_0, const OctaSym& symmetry, const IntervalVector& X); + + /** + * \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 z The box enclosing \f$\mathbf{f}(\sigma(\psi_0(\bar{\mathbf{x}})))\f$. + * \param JJf 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& z, const IntervalMatrix& JJf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose = false); vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false); From f916d147c08fe5927e494932ea8ed04802d44b16 Mon Sep 17 00:00:00 2001 From: godardma Date: Thu, 9 Oct 2025 11:18:53 +0200 Subject: [PATCH 14/24] [peibos] finished Doxygen doc + renaming variables --- python/src/core/peibos/codac2_py_peibos.cpp | 18 ++++----- src/core/peibos/codac2_peibos.cpp | 45 ++++++++++----------- src/core/peibos/codac2_peibos.h | 42 ++++++++++++++----- 3 files changed, 64 insertions(+), 41 deletions(-) diff --git a/python/src/core/peibos/codac2_py_peibos.cpp b/python/src/core/peibos/codac2_py_peibos.cpp index f0acc8580..ee3087cf9 100644 --- a/python/src/core/peibos/codac2_py_peibos.cpp +++ b/python/src/core/peibos/codac2_py_peibos.cpp @@ -34,27 +34,27 @@ void export_peibos(py::module& m) "f"_a, "X"_a); m.def("parallelepiped_inclusion", - [](const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_punc, const py::object& psi_0, const OctaSym& symmetry, const IntervalVector& X) + [](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, JJf, Jf_punc, cast>(psi_0), symmetry, 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, "JJf"_a, "Jf_punc"_a, "psi_0"_a, "symmetry"_a, "X"_a); + "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& symmetries, double epsilon, bool verbose = false) + [](const py::object& f, const py::object& psi_0, const vector& Sigma, double epsilon, bool verbose = false) { - return PEIBOS(cast>(f), cast>(psi_0), symmetries, epsilon, verbose); + 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, "symmetries"_a, "epsilon"_a, "verbose"_a = false); + "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& symmetries, double epsilon, const Vector& offset, bool verbose = false) + [](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), symmetries, epsilon, offset, verbose); + 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, "symmetries"_a, "epsilon"_a, "offset"_a, "verbose"_a = false); + "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/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index e4488bb0d..ab4aa680c 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -30,15 +30,13 @@ namespace codac2 } } - double error(const IntervalVector& z, const IntervalMatrix& JJg, const Matrix& B, const IntervalVector& X) + 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 a = z.mid(); - - auto E = (z - a) + (JJg - B)*dX; + auto E = (Y - z) + (Jf - A)*dX; return E.norm().ub(); } @@ -71,49 +69,51 @@ namespace codac2 assert_release (f.input_size() < f.output_size() && "input size must be strictly lower than output size"); assert_release (X.size() == f.input_size() && "dimension of X must match input size of f"); - auto z = f.eval(X.mid()); + auto Y = f.eval(X.mid()); + auto z = Y.mid(); - Matrix B = f.diff(X.mid()).mid(); + Matrix A = f.diff(X.mid()).mid(); // Maximum error computation - double rho = error(z, f.diff(X), B, X); + double rho = error(Y, z, f.diff(X), A, X); // Inflation of the parallelepiped - auto A = inflate_flat_parallelepiped(B, X.rad(), rho); + auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); - return Parallelepiped(z.mid(), A); + return Parallelepiped(z, A_inf); } - Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X) + 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 JJg = JJf * (sigma.permutation_matrix().template cast()) * psi_0.diff(X); + IntervalMatrix Jg = Jf * (sigma.permutation_matrix().template cast()) * psi_0.diff(X); - // B is an approximation of the Jacobian of g at the center of X - Matrix B = (Jf_tild * sigma.permutation_matrix() * (psi_0.diff(X.mid()).mid())); + 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(z, JJg, B, X); + double rho = error(Y, z, Jg, A, X); // Inflation of the parallelepiped - auto A = inflate_flat_parallelepiped(B, X.rad(), rho); + auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); - return Parallelepiped(z.mid(), A); + return Parallelepiped(z, A_inf); } - vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose) + vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& Sigma, double epsilon, bool verbose) { - return PEIBOS(f, psi_0, symmetries, epsilon, Vector::Zero(psi_0.output_size()), 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& symmetries, double epsilon, const Vector& offset, bool 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 (symmetries.size() > 0 && (int) symmetries[0].size() == psi_0.output_size() && "no generator given or wrong dimension of generator (must match output size of psi_0)"); + 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(); @@ -122,10 +122,10 @@ namespace codac2 vector boxes; double true_eps = split(IntervalVector::constant(m,{-1,1}), epsilon, boxes); - for (const auto& symmetry : symmetries) + for (const auto& sigma : Sigma) { VectorVar x(m); - AnalyticFunction g_i ({x}, f(symmetry(psi_0(x))+offset)); + AnalyticFunction g_i ({x}, f(sigma(psi_0(x))+offset)); for (const auto& X : boxes) output.push_back(parallelepiped_inclusion(g_i, X)); @@ -135,7 +135,6 @@ namespace codac2 { printf("\nPEIBOS statistics:\n"); printf("------------------\n"); - printf("Number of symmetries: %ld\n", symmetries.size()); printf("Real epsilon: %.4f\n", true_eps); printf("Computation time: %.4fs\n\n", (double)(clock()-t_start)/CLOCKS_PER_SEC); } diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index 0a629326e..17f71f90b 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -31,12 +31,13 @@ namespace codac2 /** * \brief Compute the error term for the parallelepiped inclusion. The error is later used to inflate the flat parallelepiped \f$\mathbf{a}+\mathbf{B}\cdot(\left[\mathbf{x}\right]-\bar{\mathbf{x}})\f$ * - * \param z The box enclosing \f$\mathbf{f}(\bar{\mathbf{x}})\f$. \f$\mathbf{a}\f$ is its center - * \param JJg The interval Jacobian matrix containing \f$\frac{d\mathbf{f}}{d\mathbf{x}}([\mathbf{x}])\f$ - * \param B Any Matrix with the same dimension as JJg. Usually we pick an approximation of \f$\frac{d\mathbf{f}}{d\mathbf{x}}(\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& z, const IntervalMatrix& JJg, const Matrix& B, const IntervalVector& X); + 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$. @@ -62,8 +63,8 @@ namespace codac2 /** * \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 z The box enclosing \f$\mathbf{f}(\sigma(\psi_0(\bar{\mathbf{x}})))\f$. - * \param JJf 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 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 @@ -71,8 +72,31 @@ namespace codac2 * * \return A Parallelepiped enclosing \f$\mathbf{g}([\mathbf{x}])\f$ */ - Parallelepiped parallelepiped_inclusion(const IntervalVector& z, const IntervalMatrix& JJf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X); + Parallelepiped parallelepiped_inclusion(const IntervalVector& Y, const IntervalMatrix& Jf, const Matrix& Jf_tild, const AnalyticFunction& psi_0, const OctaSym& sigma, const IntervalVector& X); - vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, bool verbose = false); - vector PEIBOS(const AnalyticFunction& f, const AnalyticFunction& psi_0, const vector& symmetries, double epsilon, const Vector& offset, bool verbose = false); + /** + * \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 From 4b0e2d91bbe9ddaf6788246e032bf6048f90867d Mon Sep 17 00:00:00 2001 From: godardma Date: Thu, 9 Oct 2025 11:31:07 +0200 Subject: [PATCH 15/24] [peibos] correcting doc --- src/core/peibos/codac2_peibos.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index 17f71f90b..a069878fa 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -29,7 +29,7 @@ namespace codac2 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{a}+\mathbf{B}\cdot(\left[\mathbf{x}\right]-\bar{\mathbf{x}})\f$ + * \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$ From 33e0de51d37306529dd2df0ed9e130b335900405 Mon Sep 17 00:00:00 2001 From: godardma Date: Thu, 9 Oct 2025 15:23:09 +0200 Subject: [PATCH 16/24] [peibos] modification to preserve the guarantee (it was not the last one) --- src/core/peibos/codac2_peibos.cpp | 9 +++++---- src/core/peibos/codac2_peibos.h | 3 ++- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index ab4aa680c..113abfb91 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -52,16 +52,17 @@ namespace codac2 Matrix A_tild (n,n); A_tild << A, N; - Matrix Q = (A_tild.transpose() * A_tild).inverse(); + IntervalMatrix A_tild_int (A_tild); + IntervalMatrix Q = inverse_enclosure(A_tild_int.transpose() * A_tild_int); - Matrix mult = Matrix::Zero(n,n); + IntervalMatrix mult = Matrix::Zero(n,n); for (int i = 0; i < n; i++) - mult(i,i) = rho*std::sqrt(Q(i,i)); + mult(i,i) = rho*sqrt(Q(i,i)); for (int i = 0; i < m; i++) mult(i,i) += e_vec(i); - return A_tild*mult; + return A_tild * mult.ub(); } Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index a069878fa..feb9cf10c 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -12,7 +12,8 @@ #include "codac2_Parallelepiped.h" #include "codac2_AnalyticFunction.h" #include "codac2_OctaSym.h" -#include +#include "codac2_inversion.h" +#include "codac2_OctaSym_operator.h" namespace codac2 From a8d32cd068931f891e66439b680ec399a0e01877 Mon Sep 17 00:00:00 2001 From: godardma Date: Thu, 9 Oct 2025 20:37:24 +0200 Subject: [PATCH 17/24] [peibos] added smag and used it --- .../manual/intervals/Interval_class.rst | 2 ++ doc/manual/manual/intervals/src.cpp | 2 ++ doc/manual/manual/intervals/src.py | 2 ++ .../domains/interval/codac2_py_Interval.cpp | 6 ++++ .../interval/codac2_py_IntervalMatrixBase.h | 6 ++++ src/core/domains/codac2_Domain.h | 2 ++ src/core/domains/interval/codac2_Interval.h | 19 ++++++++++++ .../domains/interval/codac2_Interval_impl.h | 10 +++++++ ...ac2_MatrixBase_addons_IntervalMatrixBase.h | 26 ++++++++++++++++ src/core/peibos/codac2_peibos.cpp | 29 +++++++++++++----- src/core/peibos/codac2_peibos.h | 1 + .../interval/codac2_tests_Interval.cpp | 24 +++++++++++++++ .../domains/interval/codac2_tests_Interval.py | 30 +++++++++++++++++-- 13 files changed, 149 insertions(+), 10 deletions(-) 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/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/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..2ce9a85f9 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/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 index 113abfb91..041aa74f7 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -46,23 +46,38 @@ namespace codac2 Index m = A.cols(); Index n = A.rows(); - Eigen::FullPivLU lu_decomp(A.transpose()); - Eigen::MatrixXd N = lu_decomp.kernel(); + IntervalMatrix A_int (A); - Matrix A_tild (n,n); + IntvFullPivLU LUdec((IntervalMatrix) A_int.transpose()); + IntervalMatrix N = LUdec.kernel(); + + IntervalMatrix A_tild (n,n); A_tild << A, N; - IntervalMatrix A_tild_int (A_tild); - IntervalMatrix Q = inverse_enclosure(A_tild_int.transpose() * A_tild_int); + IntervalMatrix Q = inverse_enclosure(A_tild.transpose() * A_tild); - IntervalMatrix mult = Matrix::Zero(n,n); + 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 vector is an interval vector + IntervalMatrix A_inf = A_tild * mult; + + // We inflate each generator's IntervalVector by the diameter of the other generators + IntervalMatrix A_add = A_inf.diam()*Interval(-1,1); + + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) + if (i != j) + A_inf.col(i) += A_add.col(j); + + // We return the matrix of the signed magnitude of A_inf. + // Thanks to the previous inflation step, we are certain to contain every combination of generators in the resulting parallelepiped - return A_tild * mult.ub(); + return A_inf.smag(); } Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index feb9cf10c..cb3968d45 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -14,6 +14,7 @@ #include "codac2_OctaSym.h" #include "codac2_inversion.h" #include "codac2_OctaSym_operator.h" +#include "codac2_IntvFullPivLU.h" namespace codac2 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() From 9cf3885a91c3d2a05269fef339a38258f750c455 Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 10 Oct 2025 23:26:43 +0200 Subject: [PATCH 18/24] [peibos] new version for guarantee --- src/core/domains/interval/codac2_Interval.h | 4 +-- src/core/peibos/codac2_peibos.cpp | 31 +++++++++++++-------- 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/core/domains/interval/codac2_Interval.h b/src/core/domains/interval/codac2_Interval.h index 2ce9a85f9..b8b1e2a7f 100644 --- a/src/core/domains/interval/codac2_Interval.h +++ b/src/core/domains/interval/codac2_Interval.h @@ -221,8 +221,8 @@ namespace codac2 /** * \brief Returns the signed mignitude of this * - * lower bound if lower_bound > 0 - * upper bound if upper_bound < 0 + * lower bound if lower_bound > 0, + * upper bound if upper_bound < 0, * 0 otherwise. * * \return the signed mignitude of the interval diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 041aa74f7..55980bfa1 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -63,21 +63,30 @@ namespace codac2 for (int i = 0; i < m; i++) mult(i,i) += e_vec(i); - // From here we have an IntervalMatrix A_inf, meaning that each vector is an interval vector + // From here we have an IntervalMatrix A_inf, meaning that each generator is an interval vector IntervalMatrix A_inf = A_tild * mult; - // We inflate each generator's IntervalVector by the diameter of the other generators - IntervalMatrix A_add = A_inf.diam()*Interval(-1,1); + // 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)-A_inf.col(i).mid()).ub().norm(); + + IntervalMatrix Q2 = inverse_enclosure(Y.transpose() * Y); + + IntervalMatrix mult2 = IntervalMatrix::Zero(n,n); for (int i = 0; i < n; i++) - for (int j = 0; j < n; j++) - if (i != j) - A_inf.col(i) += A_add.col(j); + mult2(i,i) = 1+rho2*sqrt(Q2(i,i)); - // We return the matrix of the signed magnitude of A_inf. - // Thanks to the previous inflation step, we are certain to contain every combination of generators in the resulting parallelepiped + IntervalMatrix Y2 = Y.template cast() * mult2; - return A_inf.smag(); + return (Y2).smag(); } Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X) @@ -94,7 +103,7 @@ namespace codac2 double rho = error(Y, z, f.diff(X), A, X); // Inflation of the parallelepiped - auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); + auto A_inf = inflate_flat_parallelepiped(A, (X-X.mid()).ub(), rho); return Parallelepiped(z, A_inf); } @@ -112,7 +121,7 @@ namespace codac2 double rho = error(Y, z, Jg, A, X); // Inflation of the parallelepiped - auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); + auto A_inf = inflate_flat_parallelepiped(A, (X-X.mid()).ub(), rho); return Parallelepiped(z, A_inf); } From 5f5fdb4618c54dfb95fd353a2b7fb2612caf4dea Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 10 Oct 2025 23:32:55 +0200 Subject: [PATCH 19/24] [peibos] new version for guarantee --- src/core/peibos/codac2_peibos.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/core/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 55980bfa1..0bb16b083 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -79,14 +79,12 @@ namespace codac2 IntervalMatrix Q2 = inverse_enclosure(Y.transpose() * Y); - IntervalMatrix mult2 = IntervalMatrix::Zero(n,n); + IntervalMatrix Y2 = IntervalMatrix::Zero(n,n); for (int i = 0; i < n; i++) - mult2(i,i) = 1+rho2*sqrt(Q2(i,i)); - - IntervalMatrix Y2 = Y.template cast() * mult2; + Y2.col(i) = Y.col(i)*(1+rho2*sqrt(Q2(i,i))); - return (Y2).smag(); + return Y2.smag(); } Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X) From e3601821b32849fbf9f272dd01ad70d304892727 Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 17 Oct 2025 13:21:03 +0200 Subject: [PATCH 20/24] [graphics] removed colormap rainbow_05 --- src/graphics/styles/codac2_ColorMap.h | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/src/graphics/styles/codac2_ColorMap.h b/src/graphics/styles/codac2_ColorMap.h index 62371175b..9f0fb0314 100644 --- a/src/graphics/styles/codac2_ColorMap.h +++ b/src/graphics/styles/codac2_ColorMap.h @@ -165,22 +165,5 @@ namespace codac2 } return cmap; } - - /** - * \brief Rainbow color map with 50% saturation and 50% transparency - * - * Rainbow color map with 50% saturation and 50% transparency - */ - static ColorMap rainbow_05() - { - ColorMap cmap( Model::HSV ); - int i = 0; - for(int h = 300 ; h > 0 ; h-=10) - { - cmap[i]=Color({(float)h,50.,100.,50.},Model::HSV); - i++; - } - return cmap; - } }; } \ No newline at end of file From 53e14e26ec3fd677d9bfa721480e6a0815dd55e2 Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 17 Oct 2025 15:28:17 +0200 Subject: [PATCH 21/24] [functions] completed parallelepiped_eval with variadic argument --- examples/07_centered_2D/main.cpp | 6 +---- .../{10_peibos => 11_peibos}/CMakeLists.txt | 0 examples/{10_peibos => 11_peibos}/main.cpp | 0 examples/{10_peibos => 11_peibos}/main.py | 6 ++--- python/codac/core/__init__.py | 3 +++ .../analytic/codac2_py_AnalyticFunction.h | 7 +---- python/src/core/peibos/codac2_py_peibos.cpp | 8 ------ .../zonotope/codac2_Parallelepiped_eval.h | 20 ++++++++++---- .../analytic/codac2_AnalyticFunction.h | 3 ++- src/core/peibos/codac2_peibos.cpp | 26 +++---------------- src/core/peibos/codac2_peibos.h | 10 ------- 11 files changed, 29 insertions(+), 60 deletions(-) rename examples/{10_peibos => 11_peibos}/CMakeLists.txt (100%) rename examples/{10_peibos => 11_peibos}/main.cpp (100%) rename examples/{10_peibos => 11_peibos}/main.py (97%) diff --git a/examples/07_centered_2D/main.cpp b/examples/07_centered_2D/main.cpp index 166dc3ce6..48ca64520 100644 --- a/examples/07_centered_2D/main.cpp +++ b/examples/07_centered_2D/main.cpp @@ -17,14 +17,10 @@ int main(){ // we use Fermat's spiral AnalyticFunction f1 ({t},{a*sqrt(t)*cos(t),a*sqrt(t)*sin(t)}); - VectorVar t_vec(1); - AnalyticFunction f2 ({t_vec},{a*sqrt(t_vec[0])*cos(t_vec[0]),a*sqrt(t_vec[0])*sin(t_vec[0])}); - double dt=0.2; double time=0.0; while (time<100.0) { Interval T(time,time+dt); - IntervalVector Tvec ({T}); // using eval_ would be easier, but it's not allowed :( IntervalVector df = f1.diff(T); if (df.is_empty() || df.is_unbounded()) { @@ -41,7 +37,7 @@ int main(){ 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)},"zonotopes"}); - Parallelepiped p = parallelepiped_inclusion(f2,Tvec); + Parallelepiped p = f1.parallelepiped_eval(T); fig4.draw_parallelepiped(p,{{Color::black(),Color::green(0.1)},"parallelepipeds"}); } time = time+dt; diff --git a/examples/10_peibos/CMakeLists.txt b/examples/11_peibos/CMakeLists.txt similarity index 100% rename from examples/10_peibos/CMakeLists.txt rename to examples/11_peibos/CMakeLists.txt diff --git a/examples/10_peibos/main.cpp b/examples/11_peibos/main.cpp similarity index 100% rename from examples/10_peibos/main.cpp rename to examples/11_peibos/main.cpp diff --git a/examples/10_peibos/main.py b/examples/11_peibos/main.py similarity index 97% rename from examples/10_peibos/main.py rename to examples/11_peibos/main.py index cc871bbe9..6e4fd52b4 100644 --- a/examples/10_peibos/main.py +++ b/examples/11_peibos/main.py @@ -15,7 +15,7 @@ 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.]) + 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]) @@ -46,7 +46,7 @@ figure_3d_proj.set_window_properties([25,600],[500,500]) figure_3d_proj.set_axes(axis(0,[-1.5,2.5]), axis(1,[-2,2])) - v_par_3d = PEIBOS(f_3d,psi0_3d,[id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()],0.2) + 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)) @@ -80,7 +80,7 @@ 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) + 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)) 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/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 index ee3087cf9..64923ac5e 100644 --- a/python/src/core/peibos/codac2_py_peibos.cpp +++ b/python/src/core/peibos/codac2_py_peibos.cpp @@ -25,14 +25,6 @@ using namespace pybind11::literals; void export_peibos(py::module& m) { - m.def("parallelepiped_inclusion", - [](const py::object& f, const IntervalVector& X) - { - return parallelepiped_inclusion(cast>(f), X); - }, - PARALLELEPIPED_PARALLELEPIPED_INCLUSION_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_INTERVALVECTOR_REF, - "f"_a, "X"_a); - 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) { diff --git a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h index bfa6c8420..2f888cab6 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h @@ -11,19 +11,29 @@ #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); + 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"); - // 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/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/peibos/codac2_peibos.cpp b/src/core/peibos/codac2_peibos.cpp index 0bb16b083..1c554d3ee 100644 --- a/src/core/peibos/codac2_peibos.cpp +++ b/src/core/peibos/codac2_peibos.cpp @@ -8,6 +8,7 @@ */ #include "codac2_peibos.h" +#include "codac2_Parallelepiped_eval.h" using namespace std; using namespace codac2; @@ -75,7 +76,7 @@ namespace codac2 Interval rho2 (0.); for (int i = 0; i < n; i++) - rho2 += (A_inf.col(i)-A_inf.col(i).mid()).ub().norm(); + rho2 += A_inf.col(i).rad().norm(); IntervalMatrix Q2 = inverse_enclosure(Y.transpose() * Y); @@ -87,25 +88,6 @@ namespace codac2 return Y2.smag(); } - Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X) - { - assert_release (f.input_size() < f.output_size() && "input size must be strictly lower than output size"); - assert_release (X.size() == f.input_size() && "dimension of X must match input size of f"); - - auto Y = f.eval(X.mid()); - auto z = Y.mid(); - - Matrix A = f.diff(X.mid()).mid(); - - // Maximum error computation - double rho = error(Y, z, f.diff(X), A, X); - - // Inflation of the parallelepiped - auto A_inf = inflate_flat_parallelepiped(A, (X-X.mid()).ub(), rho); - - return Parallelepiped(z, A_inf); - } - 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) @@ -119,7 +101,7 @@ namespace codac2 double rho = error(Y, z, Jg, A, X); // Inflation of the parallelepiped - auto A_inf = inflate_flat_parallelepiped(A, (X-X.mid()).ub(), rho); + auto A_inf = inflate_flat_parallelepiped(A, X.rad(), rho); return Parallelepiped(z, A_inf); } @@ -151,7 +133,7 @@ namespace codac2 AnalyticFunction g_i ({x}, f(sigma(psi_0(x))+offset)); for (const auto& X : boxes) - output.push_back(parallelepiped_inclusion(g_i, X)); + output.push_back(g_i.parallelepiped_eval(X)); } if (verbose) diff --git a/src/core/peibos/codac2_peibos.h b/src/core/peibos/codac2_peibos.h index cb3968d45..6c736ed37 100644 --- a/src/core/peibos/codac2_peibos.h +++ b/src/core/peibos/codac2_peibos.h @@ -52,16 +52,6 @@ namespace codac2 */ Matrix inflate_flat_parallelepiped (const Matrix& A, const Vector& e_vec, double rho); - /** - * \brief Compute a parallelepiped enclosing of \f$\mathbf{f}([\mathbf{x}])\f$ when \f$\mathbf{f}:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ is an analytic function with \f$m < n\f$. - * - * \param f The analytic function \f$\mathbf{f}\f$ - * \param X The box \f$[\mathbf{x}]\f$ - * - * \return A Parallelepiped enclosing \f$\mathbf{f}([\mathbf{x}])\f$ - */ - Parallelepiped parallelepiped_inclusion(const AnalyticFunction& f, const IntervalVector& X); - /** * \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$. * From 2cbb21d414383c88630998c4d4f2c414f1dd746e Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 17 Oct 2025 16:06:16 +0200 Subject: [PATCH 22/24] [zonotope] proj is now a method of Zonotope --- .../zonotope/codac2_py_Parallelepiped.cpp | 7 ------- .../domains/zonotope/codac2_py_Zonotope.cpp | 7 +++++++ .../domains/zonotope/codac2_Parallelepiped.cpp | 17 ----------------- .../domains/zonotope/codac2_Parallelepiped.h | 9 --------- src/core/domains/zonotope/codac2_Zonotope.cpp | 17 +++++++++++++++++ src/core/domains/zonotope/codac2_Zonotope.h | 9 +++++++++ 6 files changed, 33 insertions(+), 33 deletions(-) diff --git a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp index 2fae348b8..98474dc3e 100644 --- a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp +++ b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp @@ -32,13 +32,6 @@ 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) 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/src/core/domains/zonotope/codac2_Parallelepiped.cpp b/src/core/domains/zonotope/codac2_Parallelepiped.cpp index 8a1a40f66..34a9213ad 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) diff --git a/src/core/domains/zonotope/codac2_Parallelepiped.h b/src/core/domains/zonotope/codac2_Parallelepiped.h index 95960bea8..9761e26a9 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 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 */ From 7235e7beb01d641ba5bfdebea3a562d53a7a5570 Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 17 Oct 2025 17:36:35 +0200 Subject: [PATCH 23/24] [functions] added tests for parallelepiped_eval --- .../zonotope/codac2_Parallelepiped_eval.h | 5 +- tests/CMakeLists.txt | 1 + .../zonotope/codac2_tests_Parallelepiped.cpp | 3 +- .../zonotope/codac2_tests_Parallelepiped.py | 2 +- .../codac2_tests_Parallelepiped_eval.cpp | 60 +++++++++++++++++++ .../codac2_tests_Parallelepiped_eval.py | 56 +++++++++++++++++ 6 files changed, 123 insertions(+), 4 deletions(-) create mode 100644 tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp create mode 100644 tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py diff --git a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h index 2f888cab6..2f36e0d31 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h @@ -21,7 +21,10 @@ namespace codac2 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->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."); auto Y = this->eval(x.mid()...); auto z = Y.mid(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 90f789147..95c3e1d11 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 diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp index cf5214f0c..f8d5d4dfc 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; diff --git a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py index 196ed764c..72fa047a3 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) 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 From 7765633d8a467dd30994c606a833bd88cf449e5d Mon Sep 17 00:00:00 2001 From: godardma Date: Fri, 17 Oct 2025 18:52:06 +0200 Subject: [PATCH 24/24] [peibos] added tests --- .../zonotope/codac2_py_Parallelepiped.cpp | 4 + .../zonotope/codac2_Parallelepiped.cpp | 9 ++ .../domains/zonotope/codac2_Parallelepiped.h | 9 ++ .../zonotope/codac2_Parallelepiped_eval.h | 2 +- tests/CMakeLists.txt | 2 + .../zonotope/codac2_tests_Parallelepiped.cpp | 1 + .../zonotope/codac2_tests_Parallelepiped.py | 1 + tests/core/peibos/codac2_tests_peibos.cpp | 85 +++++++++++++++++++ tests/core/peibos/codac2_tests_peibos.py | 83 ++++++++++++++++++ 9 files changed, 195 insertions(+), 1 deletion(-) create mode 100644 tests/core/peibos/codac2_tests_peibos.cpp create mode 100644 tests/core/peibos/codac2_tests_peibos.py diff --git a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp index 98474dc3e..63b9ff280 100644 --- a/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp +++ b/python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp @@ -38,5 +38,9 @@ void export_Parallelepiped(py::module& m) .def("box", &Parallelepiped::box, INTERVALVECTOR_PARALLELEPIPED_BOX_CONST) + .def("contains", &Parallelepiped::contains, + BOOL_PARALLELEPIPED_CONTAINS_CONST_VECTOR_REF_CONST, + "x"_a) + ; } diff --git a/src/core/domains/zonotope/codac2_Parallelepiped.cpp b/src/core/domains/zonotope/codac2_Parallelepiped.cpp index 34a9213ad..4073b57b4 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped.cpp +++ b/src/core/domains/zonotope/codac2_Parallelepiped.cpp @@ -44,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 9761e26a9..3fb791008 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped.h @@ -52,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 2f36e0d31..42fbab5b5 100644 --- a/src/core/domains/zonotope/codac2_Parallelepiped_eval.h +++ b/src/core/domains/zonotope/codac2_Parallelepiped_eval.h @@ -24,7 +24,7 @@ namespace codac2 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."); + "Parallelepiped evaluation requires at least one input."); auto Y = this->eval(x.mid()...); auto z = Y.mid(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 95c3e1d11..c34f808ec 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -77,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/zonotope/codac2_tests_Parallelepiped.cpp b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp index f8d5d4dfc..138df45b3 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.cpp @@ -17,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 72fa047a3..888a9c85e 100644 --- a/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py +++ b/tests/core/domains/zonotope/codac2_tests_Parallelepiped.py @@ -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/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