Codac2 polytope

CMakeLists.txt

@@ -210,4 +210,4 @@
   set(CPACK_DEBIAN_PACKAGE_HOMEPAGE ${CODAC_URL})
   # todo: finish deb package
 
-  include(CPack)
+  include(CPack)

doc/manual/manual/geometry/facets.rst

@@ -0,0 +1,33 @@
+.. _sec-polytope:
+
+Facet classes
+=============
+
+  Main author: `Damien Massé <https://labsticc.fr/fr/annuaire/masse-damien>`_
+
+Facets
+------
+
+.. doxygentypedef:: codac2::Facet
+  :project: codac
+
+.. doxygennamespace:: codac2::Facet_
+  :project: codac
+  :members:
+
+.. doxygenclass:: codac2::FacetBase
+  :project: codac
+  :members:
+
+.. doxygenclass:: codac2::FacetRhs
+  :project: codac
+  :members:
+
+
+Collection of facets
+--------------------
+
+.. doxygenclass:: codac2::CollectFacets
+  :project: codac
+  :members:
+

doc/manual/manual/geometry/index.rst

@@ -3,7 +3,7 @@
 Geometry
 ========
 
-Codac provides a set of utility functions for basic 2d geometric calculations.
+Codac provides a set of utility functions for basic geometric calculations.
 
 Several basic geometric types are available under the form of classes representing enclosures of these types. For instance, a 2d point if represented by an ``IntervalVector`` enclosing it, a segment between two 2d points is implemented by the ``Segment`` class, the endpoints of which are ``IntervalVector`` objects. Furthermore, ``Polygon`` and ``ConvexPolygon`` are also available. Operations between these structure will reliably meet the uncertainties associated with the enclosure of points (vertices), as well as floating-point calculations.
 
@@ -16,6 +16,7 @@ Because computations are based on interval arithmetic, all these functions provi
   segment.rst
   polygon.rst
   zonotope.rst
+  polytope.rst
 
 Related interval enumerations
 -----------------------------

doc/manual/manual/geometry/polytope.rst

@@ -0,0 +1,60 @@
+.. _sec-polytope:
+
+Polytope
+========
+
+  Main author: `Damien Massé <https://labsticc.fr/fr/annuaire/masse-damien>`_
+
+The Polytope class represent convex polytopes (or, more generally, 
+convex polyhedra).
+Polytopes are internally represented as intersections of linear constraints
+(and a bounding box), in a class name CollectFacets.
+A double-description algorithm is used to compute generators for operations
+requiring them (union, projections...). Due to the imprecise nature of
+floating-point computations, the algorithm is designed such that the
+convex hull of the "generators" encloses the polytope. However, no guarantee
+is given that each vertice of the polytope is associated to a generator.
+As a result,
+going back to the hyperplane representation may return a larger polytope.
+
+.. toctree::
+  :maxdepth: 1
+
+  facets.rst
+  polytope_class.rst
+
+
+Building a polytope
+-------------------
+Polytopes can be defined from a set of facets (linear inequalities 
+or equalities), or a set of vertices. 
+
+.. tabs::
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src_polytope.cpp
+      :language: c++
+      :start-after: [polytope-1-beg]
+      :end-before: [polytope-1-end]
+      :dedent: 2
+
+Output of a polytope
+--------------------
+Polytopes are displayed as a set of linear constraints and a bounding box.
+To display the set of vertices, one can use the 
+vertices() method.
+
+.. tabs::
+
+  .. group-tab:: C++
+
+    .. literalinclude:: src_polytope.cpp
+      :language: c++
+      :start-after: [polytope-2-beg]
+      :end-before: [polytope-2-end]
+      :dedent: 2
+
+Operations on polytopes
+-----------------------
+

doc/manual/manual/geometry/polytope_class.rst

@@ -0,0 +1,14 @@
+.. _sec-polytope:
+
+The Polytope class
+==================
+
+  Main author: `Damien Massé <https://labsticc.fr/fr/annuaire/masse-damien>`_
+
+Polytope
+--------
+
+.. doxygenclass:: codac2::Polytope
+  :project: codac
+  :members:
+

doc/manual/manual/geometry/src_polytope.cpp

@@ -0,0 +1,83 @@
+// Author : Damien Massé
+// Graphical illustration of the polytope test
+
+#include <codac>
+#include <vector>
+#include <memory>
+
+using namespace std;
+using namespace codac2;
+
+void draw_polytope(Figure3D &fig, const Polytope &P, const StyleProperties &st) {
+   std::vector<std::vector<Vector>> facets3D=P.vertices_3Dfacets();
+   Vector center = Vector::zero(3);
+   Matrix transfo = Matrix::Identity(3,3);
+   for (const std::vector<Vector> &vec : facets3D) {
+      fig.draw_polygon(center,transfo,vec,st);
+   }
+}
+
+int main(int argc, char *argv[])
+{
+//   std::cout << std::scientific << std::setprecision(20);
+
+   Figure3D fig("manual_polytope");
+
+   fig.draw_axes();
+
+   // [polytope-1-beg]
+   // definition using vertices
+   std::vector<Vector> vertices 
+       { {0,0,0}, {2,1,0}, {-1,1,0}, {0,1.5,2}, {-1.5,0,-1.5}, {0,3,0} };
+   Polytope p1(vertices); // vertices {0,0,0}, {2,1,0}, etc.
+   // definition using facets
+   std::vector<std::pair<Row,double>> facets 
+       { { {1,1,0.75}, 3 },        // x+y+0.75z <= 3
+         { {-1,-0.4,0.6} ,0.6 },   // -x-0.4y+0.6z <= 0.6
+         { {-1,0.5,0.375} ,1.5 },  // ...
+         { {-1,0.5,0} ,1.5 },
+         { {0.5,-1,0.75} ,0 },
+         { {-0.75,-1,0.75} ,0 },
+         { {0.5,-1,-0.5} ,0 },
+         { {1.0/3.0,1.0/3.0,-1} ,1 } };
+   // the first argument is a bounding box, here the whole space
+   Polytope p2(IntervalVector(3), facets);
+   // p1 and p2 are almost the same polytope
+   // [polytope-1-end]
+
+
+   // [polytope-2-beg]
+   std::cout << p1 << std::endl;
+   /* output, each facet is a row (sequences of double) and a right-hand-side
+    Polytope(bbox [ [-1.5, 2] ; [0, 3] ; [-1.5, 2] ]) : 
+    Collectfacets : 8 facets
+    1 :    1    1 0.75<=3
+    2 :        -1 -0.399999       0.6<=0.600001
+    3 :    -1   0.5 0.375<=1.5
+    4 :          -1         0.5 1.77636e-16<=1.5
+    5 :  0.5   -1 0.75<=0
+    6 :    -0.75       -1 0.750001<=5.92119e-16
+    7 :  0.5   -1 -0.5<=0
+    8 : 0.333334 0.333334       -1<=1
+     end Collectfacets
+    EndPolytope
+   */
+   std::vector<IntervalVector> vp2 = p2.vertices(); /* hull of p2 */
+   for (auto &v : vp2) std::cout << v << std::endl;
+   /* output : "vertices" are boxes which enclose the polytope, but each
+      box may not enclose a "real vertice". E.g. the vertice (0,0,0) is
+      not inside the second box displayed.
+    [ [1.99999, 2.00001] ; [0.999999, 1.00001] ; [-2.4743e-14, -1.70641e-15] ]
+    [ [-7.76463e-15, 1.69539e-14] ; [-3.15401e-14, -2.88916e-15] ; [-2.52802e-14, -2.40763e-15] ]
+    [ [-1.50001, -1.5] ; [-1.16016e-14, -9.33785e-15] ; [-1.50001, -1.5] ]
+    [ [-7.11673e-15, 1.609e-14] ; [1.49999, 1.50001] ; [1.99999, 2.00001] ]
+    [ [-2.47539e-15, 9.90153e-16] ; [2.99999, 3.00001] ; [-2.31036e-15, 3.30051e-15] ]
+    [ [-1.00001, -0.999999] ; [1, 1.00001] ; [1.83361e-15, 3.20883e-15] ]
+   */
+   // [polytope-2-end]
+
+   fig.draw_polytope(p1,StyleProperties(Color::dark_red(0.8),"p1"));
+   fig.draw_polytope(p2,StyleProperties(Color::dark_blue(0.8),"p2"));
+
+   return 0;
+}

doc/manual/manual/visualization/3d_visualization.rst

@@ -38,6 +38,7 @@ Geometric shapes
   - Box
   - Parallelepiped
   - Zonotope
+  - Polytope
   - Arrow
   - Parametric surface
   - Sphere
@@ -85,6 +86,9 @@ generally a sequence of adjacent triangles sharing a same vertex.
 .. doxygenfunction:: codac2::Figure3D::draw_zonotope(const Zonotope&, const StyleProperties&)
   :project: codac
 
+.. doxygenfunction:: codac2::Figure3D::draw_polytope(const Polytope&, const StyleProperties&)
+  :project: codac
+
 .. doxygenfunction:: codac2::Figure3D::draw_arrow(const Vector&, const Matrix& A, const StyleProperties&)
   :project: codac
 

doc/manual/manual/visualization/functions.rst

@@ -65,6 +65,7 @@ Geometric shapes
   - Polygone
   - Parallelepiped
   - Zonotope
+  - Polytope
   - Pie
   - Ellipse
   - Ellipsoid
@@ -122,6 +123,9 @@ Geometric shapes
 .. doxygenfunction:: codac2::Figure2D::draw_zonotope(const Zonotope&, const StyleProperties&)
   :project: codac
 
+.. doxygenfunction:: codac2::Figure2D::draw_polytope(const Polytope&, const StyleProperties&)
+  :project: codac
+
 .. doxygenfunction:: codac2::Figure2D::draw_pie(const Vector&, const Interval&, const Interval&, const StyleProperties&)
   :project: codac
 

examples/07_centered_2D/main.cpp

@@ -39,6 +39,8 @@ int main(){
 
       Parallelepiped p = f1.parallelepiped_eval(T);
       fig4.draw_parallelepiped(p,{{Color::black(),Color::green(0.1)},"parallelepipeds"});
+      Polytope p2 = f1.polytope_eval(T);
+      fig4.draw_polytope(p2,{{Color::black(),Color::red(0.1)},"polytopes"});
      }
      time = time+dt;
   }

examples/08_centered_3D/main.cpp

@@ -5,6 +5,7 @@
 using namespace std;
 using namespace codac2;
 
+
 int main()
 {
   VectorVar x(2);
@@ -47,6 +48,29 @@ int main()
       }
       phi=phi+dphi;
   }
+  phi=0.0;
+  while (phi<2*PI) {
+      psi=0.0;
+      while (psi<2*PI) {
+        IntervalVector T { {phi,phi+dphi}, {psi,psi+dpsi} };
+        Polytope P = f.polytope_eval(T);
+        fig_zon.draw_polytope(P, { Color::blue(0.5), "polytope" });
+        psi=psi+dpsi;
+      }
+      phi=phi+dphi;
+  }
+
+  phi=0.0;
+  while (phi<2*PI) {
+      psi=0.0;
+      while (psi<2*PI) {
+        IntervalVector T { {phi,phi+dphi}, {psi,psi+dpsi} };
+        IntervalVector Q = f.eval(T);
+        fig_zon.draw_box(Q, { Color::purple(0.5), "box" });
+        psi=psi+dpsi;
+      }
+      phi=phi+dphi;
+  }
 
   fig_zon.draw_surface({0,0,0}, Matrix::Identity(3,3),
 		 Interval(0,2*PI),
@@ -58,6 +82,6 @@ int main()
           2.0*sin(phi)*cos(phi)/(1+(sin(phi)*sin(phi)))
       }; 
     } ,
-		{ Color::green(0.6), "surface" });
+		{ Color::green(0.8), "surface" });
 
 }

examples/polytope_examples/CMakeLists.txt

@@ -0,0 +1,44 @@
+# ==================================================================
+#  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 IBEX
+
+# In case you installed IBEX in a local directory, you need 
+# to specify its path with the CMAKE_PREFIX_PATH option.
+# set(CMAKE_PREFIX_PATH "~/ibex-lib/build_install")
+
+find_package(IBEX REQUIRED)
+ibex_init_common() # IBEX should have installed this function
+message(STATUS "Found IBEX version ${IBEX_VERSION}")
+
+# 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}")
+
+# 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()
+
+set(PROGS main-F2V main-3Dgraphics main-V2F main-3Dhull main-3Dtest main-manual main-contractor)
+
+foreach(PROG ${PROGS})
+  add_executable(${PROG} ${PROG}.cpp)
+  target_compile_options(${PROG} PUBLIC ${CODAC_CXX_FLAGS})
+  target_include_directories(${PROG} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS})
+  target_link_libraries(${PROG} PUBLIC ${CODAC_LIBRARIES} Ibex::ibex)
+endforeach()

examples/polytope_examples/bla.ine

@@ -0,0 +1,17 @@
+H-representation
+begin
+ 12 7 real
+  0  0  1  1  0  0 -1
+  0  1  0  1  0 -1  0
+  0  1  1  0 -1  0  0
+  0  1  0 -1  0  1  0
+  0  1 -1  0  1  0  0
+  0  0  0  0  1  1 -1
+  0  0  1 -1  0  0  1
+  0  0  0  0  1 -1  1
+  0 -1  1  0  1  0  0
+  0  0  0  0 -1  1  1
+  0  0 -1  1  0  0  1
+  0 -1  0  1  0  1  0
+end
+

examples/polytope_examples/bug45.ine

@@ -0,0 +1,14 @@
+* This is to check whether issue #45 is fixed.
+* Only the origin satisfies all equations/inequalities.
+* The correct output of "cddexec --redcheck <examples/bug45.ine"
+* is in bug45res.ine.
+H-representation
+linearity 2 1 2
+begin
+ 5 4 real
+ 0.0 0.0 1.0 1.0
+ -1.0 0.0 1.0 1.0
+ -1.0 2.0 -1.0 0.0
+ -1.0 4.0 -2.0 0.0
+ 1 0 0 0
+end