Adding peibos

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)
 
 

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]
   }

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]
 

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
 ------

examples/07_centered_2D/main.cpp

@@ -35,7 +35,10 @@ int main(){
        Vector inflationbox = cent.rad() + T.rad()*df.rad();
       Matrix v (2,3);
       v << (T.rad()*df.mid()), Vector({ inflationbox[0], 0.0 }), Vector({ 0.0, inflationbox[1] });
-       fig4.draw_zonotope({cent.mid(),v},{Color::red(),Color::yellow(0.1)});
+       fig4.draw_zonotope({cent.mid(),v},{{Color::red(),Color::yellow(0.1)},"zonotopes"});
+
+      Parallelepiped p = f1.parallelepiped_eval(T);
+      fig4.draw_parallelepiped(p,{{Color::black(),Color::green(0.1)},"parallelepipeds"});
      }
      time = time+dt;
   }

examples/11_peibos/CMakeLists.txt

@@ -0,0 +1,34 @@
+# ==================================================================
+#  codac / basics example - cmake configuration file
+# ==================================================================
+
+  cmake_minimum_required(VERSION 3.5)
+  project(codac_example LANGUAGES CXX)
+
+  set(CMAKE_CXX_STANDARD 20)
+  set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Adding Codac
+
+  # In case you installed Codac in a local directory, you need 
+  # to specify its path with the CMAKE_PREFIX_PATH option.
+  # set(CMAKE_PREFIX_PATH "~/codac/build_install")
+
+  find_package(CODAC REQUIRED)
+  message(STATUS "Found Codac version ${CODAC_VERSION}")
+
+# Initializating Ibex
+
+  ibex_init_common()
+
+# Compilation
+
+  if(FAST_RELEASE)
+    add_compile_definitions(FAST_RELEASE)
+    message(STATUS "You are running Codac in fast release mode. (option -DCMAKE_BUILD_TYPE=Release is required)")
+  endif()
+
+  add_executable(${PROJECT_NAME} main.cpp)
+  target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS})
+  target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS})
+  target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES})

examples/11_peibos/main.cpp

@@ -0,0 +1,105 @@
+#include <codac>
+
+using namespace std;
+using namespace codac2;
+
+int main()
+{
+  // 2D example of the PEIBOS algorithm
+  VectorVar y_2d(2);
+  double a = 1.4; double b = 0.3;
+  AnalyticFunction f_2d({y_2d},{y_2d[1]+1-a*sqr(y_2d[0]),b*y_2d[0]});
+
+  VectorVar X_2d(1);
+  AnalyticFunction psi0_2d ({X_2d},{cos(X_2d[0]*PI/4.-PI/2),sin(X_2d[0]*PI/4.-PI/2)});
+
+  OctaSym id_2d ({1,2});
+  OctaSym s ({-2,1});
+
+  auto v_par_2d = PEIBOS(f_2d, psi0_2d, {id_2d,s,s*s,s.invert()}, 0.2, {-0.2,0.}, true);
+
+  Figure2D figure_2d ("Henon Map", GraphicOutput::VIBES);
+  figure_2d.set_window_properties({25,50},{500,500});
+  figure_2d.set_axes({0,{-1.4,2.2}}, {1,{-0.4,0.3}});
+
+  for (const auto& p : v_par_2d)
+  {
+
+    figure_2d.draw_parallelepiped(p, {Color::green(),Color::green(0.5)});
+    figure_2d.draw_box(p.box(), {Color::blue()});
+    for (const auto& vertice : p.vertices())
+      figure_2d.draw_point(vertice, {Color::red(),Color::red(0.5)});
+  } 
+
+  // 3D example of the PEIBOS algorithm
+  VectorVar y_3d(3);
+  AnalyticFunction f_3d({y_3d},{sqr(y_3d[0])-sqr(y_3d[1])+y_3d[0],2*y_3d[0]*y_3d[1]+y_3d[1],y_3d[2]});
+
+  VectorVar X_3d(2);
+  AnalyticFunction psi0_3d ({X_3d},{1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))});
+
+  OctaSym id_3d ({1,2,3});
+  OctaSym s1 ({-2,1,3});
+  OctaSym s2 ({3,2,-1});
+
+  Figure3D figure3d ("Conform");
+  figure3d.draw_axes();
+
+  Figure2D figure_3d_proj ("Conform projected", GraphicOutput::VIBES);
+  figure_3d_proj.set_window_properties({25,600},{500,500});
+  figure_3d_proj.set_axes({0,{-1.5,2.5}}, {1,{-2,2}});
+
+  auto v_par_3d = PEIBOS(f_3d, psi0_3d, {id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()}, 0.2, true);  
+
+  for (const auto& p : v_par_3d)
+  {
+    figure3d.draw_parallelepiped(p, Color::green(0.5));
+    figure_3d_proj.draw_zonotope(p.proj({0,1}) , {Color::black(),Color::green(0.2)});
+  }
+
+  // nD example of the PEIBOS algorithm
+
+  VectorVar y_nd(3);
+  Matrix rot_matrix_1 ({ {1,0,0},
+                          {0,1/std::sqrt(2.0),-1/std::sqrt(2.0)},
+                          {0,1/std::sqrt(2.0),+1/std::sqrt(2.0)} });
+  Matrix rot_matrix_2 ({ {1/std::sqrt(2.0),-1/std::sqrt(2.0),0},
+                          {1/std::sqrt(2.0),+1/std::sqrt(2.0),0},
+                          {0,0,1} });
+
+  // Function mapping cube to sphere
+  AnalyticFunction g_nd ({y_nd}, {y_nd[0]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[1]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[2]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2]))});
+  // Apply two rotations
+  AnalyticFunction f_nd ({y_nd}, rot_matrix_1 * rot_matrix_2 * g_nd(y_nd));
+
+  VectorVar X_nd(1);
+  AnalyticFunction psi0_nd ({X_nd},{X_nd[0],1,1});
+
+  OctaSym id_nd ({1,2,3});
+  OctaSym s1_nd ({-2,1,3});
+  OctaSym s2_nd ({3,2,-1});
+
+  Figure3D figure_3d_nd ("Cube on Sphere");
+  figure_3d_nd.draw_axes(0.5);
+
+  Figure2D figure_2d_nd_xy ("XY Plane", GraphicOutput::VIBES);
+  figure_2d_nd_xy.set_window_properties({575,50},{500,500});
+  figure_2d_nd_xy.set_axes(axis(0,{-1.,1.}), axis(1,{-1.,1.}));
+
+  Figure2D figure_2d_nd_zy ("ZY Plane", GraphicOutput::VIBES);
+  figure_2d_nd_zy.set_window_properties({1125,50},{500,500});
+  figure_2d_nd_zy.set_axes(axis(0,{-1.,1.}), axis(1,{-1.,1.}));
+
+  // 12 symmetries are needed : id, s1, s1^2, s1^-1, the four same multiplied by s2, and the four same multiplied by s2^2
+  auto v_par_nd = PEIBOS(f_nd, psi0_nd, {id_nd,s1_nd,s1_nd*s1_nd,s1_nd.invert(),
+                                          s2_nd,s2_nd*s1_nd,s2_nd*s1_nd*s1_nd,s2_nd*s1_nd.invert(),
+                                          s2_nd*s2_nd,s2_nd*s2_nd*s1_nd,s2_nd*s2_nd*s1_nd*s1_nd,s2_nd*s2_nd*s1_nd.invert()}, 
+                                          0.1, true);
+
+  for (const auto& p : v_par_nd)
+  {
+    figure_3d_nd.draw_parallelepiped(p, Color::green(0.5));
+    figure_2d_nd_xy.draw_zonotope(p.proj({0,1}), {Color::black(),Color::green(0.2)});
+    figure_2d_nd_zy.draw_zonotope(p.proj({2,1}), {Color::black(),Color::green(0.2)});
+  }
+}

examples/11_peibos/main.py

@@ -0,0 +1,88 @@
+from codac import *
+import numpy as np
+
+if __name__=="__main__":
+
+  # 2D example of the PEIBOS algorithm
+
+  y_2d = VectorVar(2)
+  a,b = 1.4,0.3
+  f_2d = AnalyticFunction([y_2d],[y_2d[1]+1-a*sqr(y_2d[0]),b*y_2d[0]])
+
+  X_2d = VectorVar(1)
+  psi0_2d = AnalyticFunction([X_2d],[cos(X_2d[0]*PI/4.-PI/2),sin(X_2d[0]*PI/4.-PI/2)])
+
+  id_2d = OctaSym([1, 2])
+  s = OctaSym([-2, 1])
+
+  v_par_2d = PEIBOS(f_2d,psi0_2d,[id_2d,s,s*s,s.invert()],0.2,[-0.2,0.],True)
+
+  figure_2d = Figure2D("Henon Map", GraphicOutput.VIBES)
+  figure_2d.set_window_properties([25,50],[500,500])
+  figure_2d.set_axes(axis(0,[-1.4,2.2]), axis(1,[-0.4,0.3]))
+
+  for par in v_par_2d:
+    figure_2d.draw_parallelepiped(par,[Color.green(),Color.green(0.5)])
+    figure_2d.draw_box(par.box(), [Color.blue()])
+    for vertice in par.vertices():
+      figure_2d.draw_point(vertice, [Color.red(), Color.red(0.5)])
+
+  # 3D example of the PEIBOS algorithm
+
+  y_3d = VectorVar(3)
+  f_3d = AnalyticFunction([y_3d],[sqr(y_3d[0])-sqr(y_3d[1])+y_3d[0],2*y_3d[0]*y_3d[1]+y_3d[1],y_3d[2]])
+
+  X_3d = VectorVar(2)
+  psi0_3d = AnalyticFunction([X_3d],[1/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[0]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1])),X_3d[1]/sqrt(1+sqr(X_3d[0])+sqr(X_3d[1]))])
+
+  id_3d = OctaSym([1, 2, 3])
+  s1 = OctaSym([-2, 1, 3])
+  s2 = OctaSym([3, 2, -1])
+
+  figure_3d = Figure3D("Conform")
+  figure_3d.draw_axes()
+
+  figure_3d_proj = Figure2D("Conform projected", GraphicOutput.VIBES)
+  figure_3d_proj.set_window_properties([25,600],[500,500])
+  figure_3d_proj.set_axes(axis(0,[-1.5,2.5]), axis(1,[-2,2]))
+
+  v_par_3d = PEIBOS(f_3d,psi0_3d,[id_3d,s1,s1*s1,s1.invert(),s2,s2.invert()],0.2,True)
+
+  for p in v_par_3d:
+    figure_3d.draw_parallelepiped(p,Color.green(0.5))
+    figure_3d_proj.draw_zonotope(p.proj([0, 1]) , [Color.black(),Color.green(0.2)])
+
+  # nD example of the PEIBOS algorithm
+
+  y_nd = VectorVar(3)
+  rot_matrix_1 = Matrix([[1,0,0],[0,1/np.sqrt(2.0),-1/np.sqrt(2.0)],[0,1/np.sqrt(2.0),+1/np.sqrt(2.0)]])
+  rot_matrix_2 = Matrix([[1/np.sqrt(2.0),-1/np.sqrt(2.0),0],[1/np.sqrt(2.0),1/np.sqrt(2.0),0],[0,0,1]])
+  g_nd = AnalyticFunction([y_nd], [y_nd[0]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[1]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2])), y_nd[2]/sqrt(sqr(y_nd[0])+sqr(y_nd[1])+sqr(y_nd[2]))])
+  f_nd = AnalyticFunction([y_nd], rot_matrix_1 * rot_matrix_2 * g_nd(y_nd))
+
+  X_nd = VectorVar(1)
+  psi0_nd = AnalyticFunction([X_nd],[X_nd[0],1,1])
+
+  id_nd = OctaSym([1, 2, 3])
+  s1_nd = OctaSym([-2, 1, 3])
+  s2_nd = OctaSym([3, 2, -1])
+
+  figure_3d_nd = Figure3D("Cube on Sphere")
+  figure_3d_nd.draw_axes(0.5)
+
+  figure_2d_nd_xy = Figure2D("XY Plane", GraphicOutput.VIBES)
+  figure_2d_nd_xy.set_window_properties([575,50],[500,500])
+  figure_2d_nd_xy.set_axes(axis(0,[-1., 1.]), axis(1,[-1., 1.]))
+
+  figure_2d_nd_zy = Figure2D("ZY Plane", GraphicOutput.VIBES)
+  figure_2d_nd_zy.set_window_properties([1125,50],[500,500])
+  figure_2d_nd_zy.set_axes(axis(0,[-1., 1.]), axis(1,[-1., 1.]))
+
+  v_par_nd = PEIBOS(f_nd,psi0_nd,[id_nd,s1_nd,s1_nd*s1_nd,s1_nd.invert(),
+                                  s2_nd,s2_nd*s1_nd,s2_nd*s1_nd*s1_nd,s2_nd*s1_nd.invert(),
+                                  s2_nd*s2_nd,s2_nd*s2_nd*s1_nd,s2_nd*s2_nd*s1_nd*s1_nd,s2_nd*s2_nd*s1_nd.invert()],0.1,True)
+
+  for p in v_par_nd:
+    figure_3d_nd.draw_parallelepiped(p, Color.green(0.5))
+    figure_2d_nd_xy.draw_zonotope(p.proj([0, 1]), [Color.black(), Color.green(0.2)])
+    figure_2d_nd_zy.draw_zonotope(p.proj([2, 1]), [Color.black(), Color.green(0.2)])

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)

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/

python/src/core/codac2_py_core.cpp

@@ -120,6 +120,9 @@ void export_operators(py::module& m);
 // paver
 void export_pave(py::module& m);
 
+// peibos
+void export_peibos(py::module& m);
+
 // separators
 py::class_<SepBase,pySep> export_Sep(py::module& m);
 void export_SepAction(py::module& m, py::class_<SepBase,pySep>& pysep);
@@ -273,6 +276,9 @@ PYBIND11_MODULE(_core, m)
   // paver
   export_pave(m);
 
+  // peibos
+  export_peibos(m);
+
   // separators
   auto py_sep = export_Sep(m);
   export_SepAction(m,py_sep);

python/src/core/domains/interval/codac2_py_Interval.cpp

@@ -87,6 +87,12 @@ py::class_<Interval> 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)
 

python/src/core/domains/interval/codac2_py_IntervalMatrixBase.h

@@ -55,6 +55,12 @@ void export_IntervalMatrixBase(py::module& m, py::class_<S>& 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)
 

python/src/core/domains/zonotope/codac2_py_Parallelepiped.cpp

@@ -32,18 +32,15 @@ void export_Parallelepiped(py::module& m)
       PARALLELEPIPED_PARALLELEPIPED_CONST_VECTOR_REF_CONST_MATRIX_REF,
       "z"_a, "A"_a)
 
-    .def("proj",[](const Parallelepiped& x, const std::vector<Index_type>& indices)
-        {
-          return x.proj(matlab::convert_indices(indices));
-        },
-      ZONOTOPE_PARALLELEPIPED_PROJ_CONST_VECTOR_INDEX_REF_CONST,
-      "indices"_a)
-
     .def("vertices", &Parallelepiped::vertices,
       VECTOR_VECTOR_PARALLELEPIPED_VERTICES_CONST)
 
     .def("box", &Parallelepiped::box,
       INTERVALVECTOR_PARALLELEPIPED_BOX_CONST)
 
+    .def("contains", &Parallelepiped::contains,
+      BOOL_PARALLELEPIPED_CONTAINS_CONST_VECTOR_REF_CONST,
+      "x"_a)
+
   ;
 }

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<Index_type>& 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)