diff --git a/CMakeLists.txt b/CMakeLists.txt index 7eae26a62..cf4943e4d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -210,4 +210,4 @@ set(CPACK_DEBIAN_PACKAGE_HOMEPAGE ${CODAC_URL}) # todo: finish deb package - include(CPack) \ No newline at end of file + include(CPack) diff --git a/doc/manual/manual/geometry/facets.rst b/doc/manual/manual/geometry/facets.rst new file mode 100644 index 000000000..3bb553190 --- /dev/null +++ b/doc/manual/manual/geometry/facets.rst @@ -0,0 +1,33 @@ +.. _sec-polytope: + +Facet classes +============= + + Main author: `Damien Massé `_ + +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: + diff --git a/doc/manual/manual/geometry/index.rst b/doc/manual/manual/geometry/index.rst index 8e998195c..eb328afbc 100644 --- a/doc/manual/manual/geometry/index.rst +++ b/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 ----------------------------- diff --git a/doc/manual/manual/geometry/polytope.rst b/doc/manual/manual/geometry/polytope.rst new file mode 100644 index 000000000..c899dfaa0 --- /dev/null +++ b/doc/manual/manual/geometry/polytope.rst @@ -0,0 +1,60 @@ +.. _sec-polytope: + +Polytope +======== + + Main author: `Damien Massé `_ + +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 +----------------------- + diff --git a/doc/manual/manual/geometry/polytope_class.rst b/doc/manual/manual/geometry/polytope_class.rst new file mode 100644 index 000000000..ac0ea74b5 --- /dev/null +++ b/doc/manual/manual/geometry/polytope_class.rst @@ -0,0 +1,14 @@ +.. _sec-polytope: + +The Polytope class +================== + + Main author: `Damien Massé `_ + +Polytope +-------- + +.. doxygenclass:: codac2::Polytope + :project: codac + :members: + diff --git a/doc/manual/manual/geometry/src_polytope.cpp b/doc/manual/manual/geometry/src_polytope.cpp new file mode 100644 index 000000000..1230c47e9 --- /dev/null +++ b/doc/manual/manual/geometry/src_polytope.cpp @@ -0,0 +1,83 @@ +// Author : Damien Massé +// Graphical illustration of the polytope test + +#include +#include +#include + +using namespace std; +using namespace codac2; + +void draw_polytope(Figure3D &fig, const Polytope &P, const StyleProperties &st) { + std::vector> facets3D=P.vertices_3Dfacets(); + Vector center = Vector::zero(3); + Matrix transfo = Matrix::Identity(3,3); + for (const std::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 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> 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 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; +} diff --git a/doc/manual/manual/visualization/3d_visualization.rst b/doc/manual/manual/visualization/3d_visualization.rst index d9d632159..bd9399b50 100644 --- a/doc/manual/manual/visualization/3d_visualization.rst +++ b/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 diff --git a/doc/manual/manual/visualization/functions.rst b/doc/manual/manual/visualization/functions.rst index 90216a8b3..841086745 100644 --- a/doc/manual/manual/visualization/functions.rst +++ b/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 diff --git a/examples/07_centered_2D/main.cpp b/examples/07_centered_2D/main.cpp index 48ca64520..45c2d0b60 100644 --- a/examples/07_centered_2D/main.cpp +++ b/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; } diff --git a/examples/08_centered_3D/main.cpp b/examples/08_centered_3D/main.cpp index 3b2400beb..bad719ef9 100644 --- a/examples/08_centered_3D/main.cpp +++ b/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" }); } diff --git a/examples/polytope_examples/CMakeLists.txt b/examples/polytope_examples/CMakeLists.txt new file mode 100644 index 000000000..a3113d999 --- /dev/null +++ b/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() diff --git a/examples/polytope_examples/bla.ine b/examples/polytope_examples/bla.ine new file mode 100644 index 000000000..987633e1d --- /dev/null +++ b/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 + diff --git a/examples/polytope_examples/bug45.ine b/examples/polytope_examples/bug45.ine new file mode 100644 index 000000000..f6aa4318c --- /dev/null +++ b/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 +#include +#include + +using namespace std; +using namespace codac2; + +void draw_polytope(Figure3D &fig, const Polytope &P, const StyleProperties &st) { + std::vector> facets3D=P.vertices_3Dfacets(); + Vector center = Vector::zero(3); + Matrix transfo = Matrix::Identity(3,3); + for (const std::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("3D polytopes"); + + /* box */ + IntervalVector box1 { { 1, 3 }, { 1, 2 }, { 2, 2.5 } }; + Polytope pbox1(box1); + std::cout << "Box : " << pbox1 << "\n"; + draw_polytope(fig,pbox1,StyleProperties(Color::dark_green(),"box")); + + /* flat box */ + IntervalVector box2 { { 1, 3 }, { -1, -1 }, { 2, 2.5 } }; + Polytope pbox2(box2); + std::cout << "Flat box : " << pbox2 << "\n"; + draw_polytope(fig,pbox2,StyleProperties(Color::dark_red(),"flat box")); + + /* parallelepiped */ + Polytope paral(Parallelepiped({-1,1,1}, {{ 0.5,0.6,0.0 }, + { -0.3,0.0,0.4 }, + { 0.0, -0.6, 0.8 }})); + std::cout << "Parallelotop : " << paral << "\n"; + draw_polytope(fig,paral,StyleProperties(Color::dark_blue(0.6),"parallelepiped")); + + /* zonotope */ + Polytope zono(Zonotope({-4,1,1}, {{ 0.5,0.6,0.0,0.4,0.0,1.1 }, + { -0.3,0.0,0.4,-0.5,0.3,-0.3 }, + { -0.2,0.0,0.3,-0.4,-0.2,-0.2 }})); + std::cout << "Zonotope : " << zono << "\n"; + draw_polytope(fig,zono,StyleProperties(Color::dark_yellow(0.6),"zonotope")); + + /* intersection hyperplane */ + Polytope iplan = zono.meet_with_hyperplane(2,1.5); + std::cout << "Intersection with hyperplan (iplan) : " << iplan << "\n"; + iplan.minimize_constraints(); + std::cout << "iplan (minimized) : " << iplan << "\n"; + draw_polytope(fig,iplan,StyleProperties(Color::dark_purple(),"iplan")); + + + /* intersection */ + Polytope zono2 = zono; + zono2.meet_with_polytope(paral); + + std::cout << "Intersection : " << zono2 << "\n"; + draw_polytope(fig,zono2,StyleProperties(Color::dark_gray(),"intersection")); + zono2.minimize_constraints(); + std::cout << "Intersection (minimized) : " << zono2 << "\n"; + + /* union */ + Polytope uni = Polytope::union_of_polytopes({ pbox1, pbox2, paral, zono }); + std::cout << "Union : " << uni << "\n"; + draw_polytope(fig,uni,StyleProperties(Color::dark_orange(0.3),"union")); + + /* homothety */ + uni.homothety(IntervalVector({ 0,0,-3 }), 2.0); + std::cout << "Homot : " << uni << "\n"; + draw_polytope(fig,uni,StyleProperties(Color::dark_red(0.3),"homot")); + + return 0; +} diff --git a/examples/polytope_examples/main-3Dhull.cpp b/examples/polytope_examples/main-3Dhull.cpp new file mode 100644 index 000000000..3514e8db5 --- /dev/null +++ b/examples/polytope_examples/main-3Dhull.cpp @@ -0,0 +1,56 @@ +// Author : Damien Massé +// Generate random points on a ball of radius 10 and draw the convex hull + +#include +#include +#include +#include +#include + +using namespace std; +using namespace codac2; + +void draw_polytope(Figure3D &fig, const Polytope &P, const StyleProperties &st) { + std::vector> facets3D=P.vertices_3Dfacets(); + Vector center = Vector::zero(3); + Matrix transfo = Matrix::Identity(3,3); + for (const std::vector &vec : facets3D) { + fig.draw_polygon(center,transfo,vec,st); + } +} + + +std::random_device rd; // Will be used to obtain a seed for the random number engine +std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd() +std::uniform_real_distribution<> dis(-1.0,1.0); + +double radius=10.0; + +Vector random_point() { + double a=dis(gen)*radius; + double b=dis(gen)*std::sqrt(radius*radius-a*a); + double c=(dis(gen)>0.0 ? 1.0 : -1.0) * std::sqrt(radius*radius-a*a-b*b); + return { a,b,c }; +} + +int main(int argc, char *argv[]) +{ + Figure3D fig("3D hull"); + + std::cout << std::scientific << std::setprecision(9); + + int nbvertices=20; + if (argc>1) nbvertices= std::atoi(argv[1]); + std::vector vertices; + Matrix scale = 0.1*Matrix::Identity(3,3); + for (int i=0;i +#include +#include + +using namespace std; +using namespace codac2; + +int main(int argc, char *argv[]) +{ +// std::cout << std::scientific << std::setprecision(20); + + Figure3D fig("test_polytope"); + + fig.draw_axes(); +/* polytope with vertices + (0,0,0) , (2,0,1) , (1,0.1,0) , (0,-0.1,1), (1,2,0.5) + inequalities + -x + 20y + 82z <= 80 x - 20y - 2z <= 0 + -20.5x + 10y +z <= 0 -x + 10y - 38z <= 0 + 3.9x + y - 3.8z <= 4 x -10y - 2z <= 0 */ + + std::vector> facets + { { {-1,20,82}, 80 }, + { {-20.5,10,1} ,0 }, + { {3.9,1,-3.8} ,4 }, + { {1,-20,-2} ,0 }, + { {-1,10,-38} ,0 }, + { {1,-10,-2} ,0 } }; + + Polytope p(IntervalVector({{-4,4},{-4,4},{-4,4}}), facets, true); + fig.draw_polytope(p,StyleProperties(Color::dark_red(0.8), + "initial polytope")); + + std::vector vertices = p.vertices(); + + fig.draw_box(IntervalVector({{0.2,0.4},{0.0,0.3},{0.3,0.5}}), + StyleProperties(Color::black(),"box included")); + + Polytope q = p; + q.inflate(0.5); + for (IntervalVector &vt : vertices) { + fig.draw_box(vt+IntervalVector::constant(3,Interval(-0.5,0.5)), + StyleProperties(Color::orange(0.6), + "box_inflate")); + } + fig.draw_polytope(q,StyleProperties(Color::red(0.4), + "inflate")); + Polytope r = p; + r.inflate_ball(0.5); + for (IntervalVector &vt : vertices) { + fig.draw_sphere(vt.mid(),0.5*Matrix::Identity(3,3), + StyleProperties(Color::orange(0.6), + "ball_inflate")); + } + fig.draw_polytope(r,StyleProperties(Color::purple(0.4), + "inflate_ball")); + + Polytope s = p; + s.unflat(2,0.5); + for (IntervalVector &vt : vertices) { + fig.draw_box(vt+IntervalVector({{-0.01,0.01},{-0.01,0.01},{-0.5,0.5}}), + StyleProperties(Color::black(), + "box_unflat")); + } + fig.draw_polytope(s,StyleProperties(Color::yellow(0.4), + "unflat")); + + Polytope t = p; + t.homothety(IntervalVector({0,0,-0.5}),2); + fig.draw_polytope(t,StyleProperties(Color::dark_blue(0.8), + "homothety")); + + Polytope u = t; + u.meet_with_polytope(p); + fig.draw_polytope(u,StyleProperties(Color::dark_gray(), + "intersection")); + + Polytope v = Polytope::union_of_polytopes({ p,t }); + fig.draw_polytope(v,StyleProperties(Color::cyan(0.4), + "union")); + + IntervalMatrix M { { cos(PI/3) , sin(PI/3) , 0 }, + { -sin(PI/3), cos(PI/3) , 0 }, + { 0 , 0 , 1 } }; + IntervalVector P { 1-cos(PI/3) , sin(PI/3) , 0 }; + Polytope w = p.reverse_affine_transform(M,P, + IntervalVector({{-4,4},{-4,4},{-4,4}})); + + fig.draw_polytope(w,StyleProperties(Color::blue(0.8), + "transformation")); + + return 0; +} diff --git a/examples/polytope_examples/main-F2V.cpp b/examples/polytope_examples/main-F2V.cpp new file mode 100644 index 000000000..b0f4088fc --- /dev/null +++ b/examples/polytope_examples/main-F2V.cpp @@ -0,0 +1,40 @@ +// Author : Damien Massé +// Transformation facets to vertex from a .ine file + +#include +#include +#include +#include + +using namespace std; +using namespace codac2; + +int main(int argc, char *argv[]) +{ + std::cout << std::scientific << std::setprecision(9); + + const char *namefile = "cross12.ine"; + if (argc>1) namefile= argv[1]; + Polytope pol = Polytope::from_ineFile(namefile); + std::cout << pol << std::endl; + + std::vector + lvect = pol.vertices(); + std::cout << lvect.size() << " vertices" << std::endl; + for (IntervalVector a : lvect) + std::cout << a << std::endl; + +#if 0 + IntervalVector box = IntervalVector::constant((*fc)[0]->first.row.size(), + Interval(-10,10)); + DDbuildF2V build2(box.size(),box,fc,false); + std::cout << "Initial build : " << build2 << std::endl; + for (int i=0;inbfcts();i++) { + std::cout << "add facet " << (i+1) << " : " << (*(*fc)[i]) << std::endl; + build2.add_facet(i); +// std::cout << build2; + } + std::cout << build2; +#endif + return 0; +} diff --git a/examples/polytope_examples/main-V2F.cpp b/examples/polytope_examples/main-V2F.cpp new file mode 100644 index 000000000..395ba76f9 --- /dev/null +++ b/examples/polytope_examples/main-V2F.cpp @@ -0,0 +1,29 @@ +// Author : Damien Massé +// Transformation vertices -> facets from a .ext file + +#include +#include +#include +#include + +using namespace std; +using namespace codac2; + + +int main(int argc, char *argv[]) +{ + std::cout << std::scientific << std::setprecision(9); + + const char *namefile = "irbox20-4.ext"; + if (argc>1) namefile= argv[1]; + Polytope pol = Polytope::from_extFile(namefile); + + std::cout << pol << std::endl; + + std::vector lst = pol.vertices(); + std::cout << "nb vertices : " << lst.size() << std::endl; + for (int i=0;i +#include +#include + +using namespace std; +using namespace codac2; + +int main(int argc, char *argv[]) +{ +// std::cout << std::scientific << std::setprecision(20); + + // definition using facets + std::vector> facets + { { {1,1}, 3.1 }, + { {-1,0.6} ,1.6 }, + { {0.5,0.2} ,1.4 }, + { {0.5,0.75} ,1.75 }, + { {-0.75,-1} ,1.33 }, + { {0.5,-1} ,1 }, + { {-1.0/3.0,-1} ,1 } }; + // the first argument is a bounding box, here the whole space + Polytope p({{-2,3.0},{-1,2.1}}, facets); + + SepPolytope sep(p); + +// auto pv1 = pave({{-3.0,3.0},{-3.0,3.0}},sep,0.04); + + DefaultFigure::pave({{-2.0,3.0},{-3.0,3.0}},sep,0.04); + + DefaultFigure::draw_polytope(p); + + return 0; +} diff --git a/examples/polytope_examples/main-manual.cpp b/examples/polytope_examples/main-manual.cpp new file mode 100644 index 000000000..9daad7669 --- /dev/null +++ b/examples/polytope_examples/main-manual.cpp @@ -0,0 +1,46 @@ +// Author : Damien Massé +// Graphical illustration of the polytope test + +#include +#include +#include + +using namespace std; +using namespace codac2; + +int main(int argc, char *argv[]) +{ +// std::cout << std::scientific << std::setprecision(20); + + Figure3D fig("manual_polytope"); + + fig.draw_axes(); + + // definition using vertices + std::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); + // definition using facets + std::vector> facets + { { {1,1,0.75}, 3 }, + { {-1,-0.4,0.6} ,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::Constant(3,Interval()), facets); + + // p1 and p2 are ``almost'' the same polytope + + fig.draw_polytope(p1,StyleProperties(Color::dark_red(0.8),"p1")); + std::cout << p1 << std::endl; + + std::vector vp2 = p2.vertices(); + for (auto &v : vp2) std::cout << v << std::endl; + fig.draw_polytope(p2,StyleProperties(Color::dark_blue(0.8),"p2")); + + return 0; +} diff --git a/examples/polytope_examples/sampleh8.ine b/examples/polytope_examples/sampleh8.ine new file mode 100644 index 000000000..7a8698d86 --- /dev/null +++ b/examples/polytope_examples/sampleh8.ine @@ -0,0 +1,110 @@ +* test file for redundancy removal +* the redundant rows are +* 10 11 12 14 21 23 26 32 39 40 41 44 46 +* 50 51 54 56 57 59 62 63 64 69 76 77 78 +* 79 83 84 85 87 88 91 94 97 +* LP type = 1 Seed = 123 +H-representation +begin + 100 10 integer + 0 1 0 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 0 + 0 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 0 1 + -10000 651 693 84 697 637 340 368 824 663 + -10000 725 742 387 219 751 430 202 745 356 + -10000 377 674 979 167 815 988 412 676 475 + -10000 710 275 949 284 629 1 422 974 510 + -10000 692 945 725 488 271 430 724 225 726 + -10000 465 258 450 343 87 168 161 103 919 + -10000 86 79 656 493 832 514 791 506 29 + -10000 63 630 874 918 877 272 992 119 480 + -10000 598 926 42 378 288 66 927 919 99 + -10000 256 354 106 979 641 160 395 225 837 + -10000 202 388 900 471 160 751 300 731 818 + -10000 342 502 825 563 639 261 194 984 990 + -10000 266 406 364 216 448 675 145 694 866 + -10000 362 983 732 378 134 902 946 877 205 + -10000 926 125 949 888 234 630 275 707 67 + -10000 634 81 192 768 652 822 311 961 895 + -10000 983 597 743 314 696 585 367 396 826 + -10000 511 545 539 97 111 996 477 35 372 + -10000 474 103 152 753 159 120 929 161 563 + -10000 549 793 307 456 444 184 149 792 894 + -10000 839 488 917 192 168 788 959 245 25 + -10000 750 165 338 182 392 381 962 117 713 + -10000 738 827 943 507 914 814 951 663 815 + -10000 493 339 225 351 450 788 992 167 792 + -10000 174 773 247 247 180 517 445 599 596 + -10000 303 26 967 39 535 4 7 335 217 + -10000 772 173 189 291 668 191 610 677 544 + -10000 848 642 40 125 865 100 259 534 648 + -10000 501 622 398 624 118 416 30 17 236 + -10000 218 602 697 892 322 314 361 573 985 + -10000 958 856 608 492 563 478 311 614 740 + -10000 582 913 938 949 715 338 39 726 998 + -10000 521 805 708 221 624 316 24 127 322 + -10000 491 189 412 774 418 200 193 633 315 + -10000 144 679 383 447 989 939 441 631 482 + -10000 940 241 153 215 149 457 254 207 125 + -10000 80 873 207 904 684 600 940 431 825 + -10000 75 100 353 637 432 377 940 758 164 + -10000 627 721 915 710 8 786 96 17 576 + -10000 247 104 607 432 540 164 597 282 317 + -10000 553 787 881 942 152 318 44 509 518 + -10000 181 772 942 993 916 900 713 694 812 + -10000 625 26 638 261 385 196 676 717 572 + -10000 11 911 296 761 658 631 757 817 400 + -10000 964 989 879 491 490 751 490 97 405 + -10000 28 461 343 452 480 504 411 242 892 + -10000 768 746 347 352 724 487 185 767 287 + -10000 924 228 958 554 352 460 765 742 430 + -10000 737 986 499 993 721 29 821 416 813 + -10000 483 458 249 781 584 321 779 209 570 + -10000 222 166 975 53 765 256 859 469 164 + -10000 82 696 181 166 375 145 415 496 261 + -10000 540 309 526 918 439 489 846 417 775 + -10000 699 145 327 873 467 866 711 393 603 + -10000 445 668 729 224 654 970 279 789 680 + -10000 28 659 365 306 186 520 427 176 191 + -10000 919 58 859 436 105 35 560 716 14 + -10000 641 238 340 524 205 255 272 14 472 + -10000 817 22 904 782 573 847 649 100 280 + -10000 925 277 930 495 736 960 276 726 657 + -10000 336 54 357 979 941 765 593 630 286 + -10000 65 113 104 915 222 158 580 520 188 + -10000 42 617 904 414 588 857 416 94 988 + -10000 716 49 849 67 607 99 934 184 217 + -10000 458 32 500 552 108 980 255 998 935 + -10000 860 936 115 268 231 281 948 209 766 + -10000 293 330 940 570 857 421 574 481 364 + -10000 959 949 372 757 306 757 752 312 137 + -10000 769 303 737 260 187 695 946 723 433 + -10000 980 564 605 363 776 893 558 616 603 + -10000 945 23 632 820 54 677 795 616 625 + -10000 629 692 585 33 674 703 304 489 949 + -10000 52 903 723 746 692 232 917 933 211 + -10000 367 566 364 318 607 818 424 234 605 + -10000 915 305 354 829 469 447 925 188 799 + -10000 387 596 946 888 751 614 539 295 299 + -10000 969 902 10 739 289 923 808 28 151 + -10000 537 152 546 442 182 301 258 265 762 + -10000 671 288 343 199 694 391 457 947 250 + -10000 415 734 423 435 696 111 53 227 784 + -10000 800 992 720 76 40 419 608 182 683 + -10000 619 627 278 526 453 960 772 992 86 + -10000 836 623 0 21 371 4 676 964 658 + -10000 258 677 64 968 62 219 587 308 299 + -10000 224 823 258 449 575 526 908 585 405 + -10000 593 660 663 341 36 79 50 242 293 + -10000 80 761 270 460 82 133 874 559 41 + -10000 169 539 257 827 222 743 873 929 523 + -10000 403 167 259 377 76 160 749 448 902 + -10000 290 5 276 509 342 909 193 443 552 + -10000 928 326 757 405 598 765 143 833 150 +end diff --git a/python/src/core/CMakeLists.txt b/python/src/core/CMakeLists.txt index 1b76181f1..5ed12cc16 100644 --- a/python/src/core/CMakeLists.txt +++ b/python/src/core/CMakeLists.txt @@ -32,6 +32,7 @@ contractors/codac2_py_CtcPointCloud.cpp contractors/codac2_py_CtcPolar.cpp contractors/codac2_py_CtcPolygon.cpp + contractors/codac2_py_CtcPolytope.cpp contractors/codac2_py_CtcProj.cpp contractors/codac2_py_CtcSegment.cpp contractors/codac2_py_CtcUnion.cpp @@ -52,6 +53,8 @@ domains/paving/codac2_py_Subpaving.cpp domains/zonotope/codac2_py_Parallelepiped.cpp domains/zonotope/codac2_py_Zonotope.cpp + domains/polytope/codac2_py_Facet.cpp + domains/polytope/codac2_py_Polytope.cpp domains/tube/codac2_py_Slice.h domains/tube/codac2_py_SlicedTube.h domains/tube/codac2_py_TDomain.cpp @@ -83,6 +86,7 @@ matrices/codac2_py_MatrixBase.h matrices/codac2_py_MatrixBlock.h matrices/codac2_py_Vector.cpp + matrices/codac2_py_Row.cpp matrices/codac2_py_VectorBase.h operators/codac2_py_operators.cpp @@ -104,6 +108,7 @@ separators/codac2_py_SepInverse.cpp separators/codac2_py_SepNot.cpp separators/codac2_py_SepPolygon.cpp + separators/codac2_py_SepPolytope.cpp separators/codac2_py_SepProj.cpp separators/codac2_py_SepTransform.cpp separators/codac2_py_SepUnion.cpp diff --git a/python/src/core/codac2_py_core.cpp b/python/src/core/codac2_py_core.cpp index c8cd3144a..62cc417a9 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -52,6 +52,7 @@ void export_CtcNot(py::module& m, py::class_,pyCtcInterv void export_CtcPointCloud(py::module& m, py::class_,pyCtcIntervalVector>& ctc); void export_CtcPolar(py::module& m, py::class_,pyCtcIntervalVector>& ctc); void export_CtcPolygon(py::module& m, py::class_,pyCtcIntervalVector>& ctc); +void export_CtcPolytope(py::module& m, py::class_,pyCtcIntervalVector>& ctc); void export_CtcProj(py::module& m, py::class_,pyCtcIntervalVector>& ctc); void export_CtcSegment(py::module& m, py::class_,pyCtcIntervalVector>& ctc); void export_CtcUnion(py::module& m, py::class_,pyCtcIntervalVector>& ctc); @@ -72,6 +73,8 @@ void export_PavingNode(py::module& m); void export_Subpaving(py::module& m); void export_Zonotope(py::module& m); void export_Parallelepiped(py::module& m); +void export_Facet(py::module& m); +void export_Polytope(py::module& m); void export_TDomain(py::module& m); void export_TimePropag(py::module& m); void export_TSlice(py::module& m); @@ -138,6 +141,7 @@ void export_SepInter(py::module& m, py::class_& sep); void export_SepInverse(py::module& m, py::class_& sep); void export_SepNot(py::module& m, py::class_& sep); void export_SepPolygon(py::module& m, py::class_& sep); +void export_SepPolytope(py::module& m, py::class_& sep); void export_SepProj(py::module& m, py::class_& sep); void export_SepTransform(py::module& m, py::class_& sep); void export_SepUnion(py::module& m, py::class_& sep); @@ -192,6 +196,7 @@ PYBIND11_MODULE(_core, m) export_CtcPointCloud(m, py_ctc_iv); export_CtcPolar(m, py_ctc_iv); export_CtcPolygon(m, py_ctc_iv); + export_CtcPolytope(m, py_ctc_iv); export_CtcProj(m, py_ctc_iv); export_CtcSegment(m, py_ctc_iv); export_CtcUnion(m, py_ctc_iv); @@ -200,7 +205,8 @@ PYBIND11_MODULE(_core, m) // matrices export_cart_prod(m); - py::class_ exported_row_class(m, "Row", DOC_TO_BE_DEFINED); +// py::class_ exported_row_class(m, "Row", DOC_TO_BE_DEFINED); + export_Row(m); auto py_V = export_Vector(m); auto py_M = export_Matrix(m); auto py_B = export_EigenBlock(m, "MatrixBlock"); @@ -246,6 +252,9 @@ PYBIND11_MODULE(_core, m) export_Zonotope(m); export_Parallelepiped(m); + export_Facet(m); + export_Polytope(m); + // function py::enum_(m, "EvalMode") .value("NATURAL", EvalMode::NATURAL) @@ -299,6 +308,7 @@ PYBIND11_MODULE(_core, m) export_SepInverse(m,py_sep); export_SepNot(m,py_sep); export_SepPolygon(m,py_sep); + export_SepPolytope(m,py_sep); export_SepProj(m,py_sep); export_SepTransform(m,py_sep); export_SepUnion(m,py_sep); diff --git a/python/src/core/contractors/codac2_py_CtcPolytope.cpp b/python/src/core/contractors/codac2_py_CtcPolytope.cpp new file mode 100644 index 000000000..f70030e47 --- /dev/null +++ b/python/src/core/contractors/codac2_py_CtcPolytope.cpp @@ -0,0 +1,38 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Simon Rohou, Damien Massé + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include +#include "codac2_py_Ctc.h" +#include "codac2_py_CtcPolytope_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_CtcPolytope(py::module& m, py::class_,pyCtcIntervalVector>& pyctc) +{ + py::class_ exported(m, "CtcPolytope", pyctc, CTCPOLYTOPE_MAIN); + exported + + .def(py::init(), + CTCPOLYTOPE_CTCPOLYTOPE_CONST_POLYTOPE_REF, + "p"_a) + + .def(CONTRACT_BOX_METHOD(CtcPolytope, + VOID_CTCPOLYTOPE_CONTRACT_INTERVALVECTOR_REF_CONST)) + + ; + + py::implicitly_convertible(); +} diff --git a/python/src/core/domains/polytope/codac2_py_Facet.cpp b/python/src/core/domains/polytope/codac2_py_Facet.cpp new file mode 100644 index 000000000..4aebb2837 --- /dev/null +++ b/python/src/core/domains/polytope/codac2_py_Facet.cpp @@ -0,0 +1,69 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include "codac2_py_Facet_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): +#include "codac2_py_cast.h" + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_Facet(py::module& m) +{ + py::class_ + exported1(m, "FacetBase", FACETBASE_MAIN); + exported1 + + .def(py::init(), + FACETBASE_FACETBASE_CONST_ROW_REF, "row"_a) + + .def(py::init(), + FACETBASE_FACETBASE_CONST_FACETBASE_REF,"fc"_a) + + .def("get_row",&FacetBase::get_row, + CONST_ROW_REF_FACETBASE_GET_ROW_CONST) + + .def("is_null",&FacetBase::is_null, + BOOL_FACETBASE_IS_NULL_CONST) + + .def("is_coord",&FacetBase::is_coord, + BOOL_FACETBASE_IS_COORD_CONST) + + .def("gt_dim",&FacetBase::gt_dim, + INDEX_FACETBASE_GT_DIM_CONST) + + ; + + py::class_ + exported2(m, "FacetRhs", FACETRHS_MAIN); + exported2 + + .def(py::init(), + FACETRHS_FACETRHS_DOUBLE_BOOL_INDEX, "rhs"_a, "eqcst"_a, "Id"_a) + + .def("get_rhs",&FacetRhs::get_rhs, + DOUBLE_FACETRHS_GET_RHS_CONST) + + .def("is_eqcst",&FacetRhs::is_eqcst, + BOOL_FACETRHS_IS_EQCST_CONST) + + .def("get_Id",&FacetRhs::get_Id, + INDEX_FACETRHS_GET_ID_CONST) + + ; + + m.def("Facet_make",[](const Row &row, double rhs, bool eqcst) + { return Facet_::make(row,rhs,eqcst); }); + +} diff --git a/python/src/core/domains/polytope/codac2_py_Polytope.cpp b/python/src/core/domains/polytope/codac2_py_Polytope.cpp new file mode 100644 index 000000000..a86799d5f --- /dev/null +++ b/python/src/core/domains/polytope/codac2_py_Polytope.cpp @@ -0,0 +1,165 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include "codac2_py_Polytope_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): +#include "codac2_py_cast.h" + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_Polytope(py::module& m) +{ + py::class_ + exported(m, "Polytope", POLYTOPE_MAIN); + exported + + .def(py::init<>(), + POLYTOPE_POLYTOPE) + + .def(py::init(), + POLYTOPE_POLYTOPE_INDEX,"dim"_a) + + .def(py::init(), + POLYTOPE_POLYTOPE_INDEX_BOOL, + "dim"_a,"empty"_a) + + .def(py::init(), + POLYTOPE_POLYTOPE_CONST_INTERVALVECTOR_REF, + "box"_a) + + .def(py::init(), + POLYTOPE_POLYTOPE_CONST_POLYTOPE_REF, + "P"_a) + + .def(py::init &>(), + POLYTOPE_POLYTOPE_CONST_VECTOR_VECTOR_REF, + "vertices"_a) + + .def(py::init &>(), + POLYTOPE_POLYTOPE_CONST_VECTOR_INTERVALVECTOR_REF, + "vertices"_a) + + .def(py::init &>(), + POLYTOPE_POLYTOPE_CONST_VECTOR_INTERVALVECTOR_REF, + "vertices"_a) + + .def(py::init(), + POLYTOPE_POLYTOPE_CONST_PARALLELEPIPED_REF, + "par"_a) + + .def(py::init(), + POLYTOPE_POLYTOPE_CONST_ZONOTOPE_REF, + "zon"_a) + + .def(py::init> &, bool>(), + POLYTOPE_POLYTOPE_CONST_INTERVALVECTOR_REF_CONST_VECTOR_PAIR_ROWDOUBLE_REF_BOOL, + "box"_a,"facets"_a,"minimize"_a) + + .def("copy",[](Polytope &Q) { Polytope P=Q; return P; }, + POLYTOPE_REF_POLYTOPE_OPERATORAFF_CONST_POLYTOPE_REF) + + .def("dim",&Polytope::dim, + INDEX_POLYTOPE_DIM_CONST) + + .def("dim",&Polytope::size, + INDEX_POLYTOPE_SIZE_CONST) + + .def("nb_facets",&Polytope::nb_facets, + INDEX_POLYTOPE_NB_FACETS_CONST) + + .def("nb_eq_facets",&Polytope::nb_eq_facets, + INDEX_POLYTOPE_NB_EQ_FACETS_CONST) + + .def("is_empty",[](const Polytope &P) { return P.is_empty(true); }, + BOOL_POLYTOPE_IS_EMPTY_BOOL_CONST) + + .def("is_flat",&Polytope::is_flat, + BOOL_POLYTOPE_IS_FLAT_CONST) + + .def("bound_row",&Polytope::bound_row, + BOOL_POLYTOPE_IS_FLAT_CONST, "r"_a) + + .def("contains",&Polytope::contains, + BOOLINTERVAL_POLYTOPE_CONTAINS_CONST_INTERVALVECTOR_REF_CONST, "p"_a) + + .def("is_subset",[](const Polytope &Q, const Polytope &P) + { return Q.is_subset(P); }, + VIRTUAL_BOOL_POLYTOPE_IS_SUBSET_CONST_POLYTOPE_REF_CONST, "P"_a) + + .def("box",[](const Polytope &Q) + { return Q.box(true); }, + CONST_INTERVALVECTOR_REF_POLYTOPE_BOX_BOOL_CONST) + + .def("meet_with_polytope",&Polytope::meet_with_polytope, + INT_POLYTOPE_MEET_WITH_POLYTOPE_CONST_POLYTOPE_REF,"P"_a) + + .def("inflate",[](Polytope &Q, double rad) { return Q.inflate(rad); }, + POLYTOPE_REF_POLYTOPE_INFLATE_DOUBLE,"rad"_a) + + .def("inflate",[](Polytope &Q, IntervalVector &box) + { return Q.inflate(box); }, + POLYTOPE_REF_POLYTOPE_INFLATE_CONST_INTERVALVECTOR_REF, "box"_a) + + .def("inflate_ball",[](Polytope &Q, double rad) + { return Q.inflate_ball(rad); }, + POLYTOPE_REF_POLYTOPE_INFLATE_BALL_DOUBLE,"rad"_a) + + .def("unflat",&Polytope::unflat, + POLYTOPE_REF_POLYTOPE_UNFLAT_INDEX_DOUBLE,"dm"_a,"rad"_a) + + .def("homothety",&Polytope::homothety, + POLYTOPE_REF_POLYTOPE_HOMOTHETY_CONST_INTERVALVECTOR_REF_DOUBLE, + "c"_a,"delta"_a) + + .def("minimize_constraints",&Polytope::minimize_constraints, + VIRTUAL_VOID_POLYTOPE_MINIMIZE_CONSTRAINTS_CONST) + + .def("reverse_affine_transform",&Polytope::reverse_affine_transform, + POLYTOPE_POLYTOPE_REVERSE_AFFINE_TRANSFORM_CONST_INTERVALMATRIX_REF_CONST_INTERVALVECTOR_REF_CONST_INTERVALVECTOR_REF_CONST, + "M"_a,"P"_a,"bbox"_a) + + .def("affine_transform",&Polytope::affine_transform, + POLYTOPE_POLYTOPE_AFFINE_TRANSFORM_CONST_INTERVALMATRIX_REF_CONST_INTERVALVECTOR_REF_CONST_INTERVALVECTOR_REF_BOOL_CONST, + "M"_a,"P"_a,"B"_a,"use_direct"_a) + + .def("direct_affine_transform",&Polytope::direct_affine_transform, + POLYTOPE_POLYTOPE_DIRECT_AFFINE_TRANSFORM_CONST_INTERVALMATRIX_REF_CONST_INTERVALVECTOR_REF_CONST, + "M"_a,"P"_a) + + .def("vertices",&Polytope::vertices, + VECTOR_INTERVALVECTOR_POLYTOPE_VERTICES_CONST) + + .def("vertices_3Dfacets",&Polytope::vertices_3Dfacets, + VECTOR_VECTOR_VECTOR_POLYTOPE_VERTICES_3DFACETS_CONST) + + .def_static("from_ineFile",&Polytope::from_ineFile, + STATIC_POLYTOPE_POLYTOPE_FROM_INEFILE_CONST_CHAR_PTR,"filename"_a) + + .def_static("from_extFile",&Polytope::from_extFile, + STATIC_POLYTOPE_POLYTOPE_FROM_EXTFILE_CONST_CHAR_PTR,"filename"_a) + + .def_static("union_of_polytopes", [](const std::vector &lst) + { return Polytope::union_of_polytopes(lst); }, + STATIC_POLYTOPE_POLYTOPE_UNION_OF_POLYTOPES_CONST_VECTOR_POLYTOPE_REF, + "lst"_a) + + .def("__repr__", [](const Polytope& P) { + std::ostringstream stream; + stream << P; + return string(stream.str()); + }, + OSTREAM_REF_OPERATOROUT_OSTREAM_REF_CONST_POLYTOPE_REF) + ; +} diff --git a/python/src/core/matrices/codac2_py_Row.cpp b/python/src/core/matrices/codac2_py_Row.cpp new file mode 100644 index 000000000..481a01206 --- /dev/null +++ b/python/src/core/matrices/codac2_py_Row.cpp @@ -0,0 +1,100 @@ +/** + * Row binding, from Vector bindings + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Simon Rohou, Damien Massé + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include +#include + +#include "codac2_py_doc.h" +#include "codac2_py_Matrix_addons_Base_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +//#include "codac2_py_Matrix_addons_IntervalMatrix_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_Matrix_addons_IntervalMatrixBase_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_Matrix_addons_IntervalVector_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +//#include "codac2_py_Matrix_addons_Matrix_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_Matrix_addons_MatrixBase_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_Matrix_addons_Vector_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_Matrix_addons_VectorBase_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_MatrixBase_addons_Base_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +//#include "codac2_py_MatrixBase_addons_IntervalMatrix_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_MatrixBase_addons_IntervalMatrixBase_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_MatrixBase_addons_IntervalVector_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +//#include "codac2_py_MatrixBase_addons_Matrix_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +//#include "codac2_py_MatrixBase_addons_MatrixBase_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_MatrixBase_addons_Vector_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_MatrixBase_addons_VectorBase_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_matrices_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py) +#include "codac2_py_Row_docs.h" + +#include "codac2_py_VectorBase.h" + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +py::class_ export_Row(py::module& m) +{ + py::class_ exported_row_class(m, "Row", DOC_TO_BE_DEFINED); + export_MatrixBase(m, exported_row_class); + export_VectorBase(m, exported_row_class); + + exported_row_class + + .def(py::init( + [](Index_type n) + { + matlab::test_integer(n); + return std::make_unique((int)n); + }), + DOC_TO_BE_DEFINED, + "n"_a) + + .def(py::init(), + "x"_a) + + .def(py::init( // this constructor must be the last one to be declared + [](const std::vector& v_) + { + auto v = std::make_unique(v_.size()); + for(size_t i = 0 ; i < v_.size() ; i++) + (*v)[i] = v_[i]; + return v; + }), + MATRIX_ADDONS_VECTOR_MATRIX_INITIALIZER_LIST_DOUBLE, + "v"_a) + + .def("min_coeff_index", [](const Row& x) + { + return matlab::output_index(x.min_coeff_index()); + }, + MATRIXBASE_ADDONS_VECTOR_INDEX_MIN_COEFF_INDEX_CONST) + + .def("max_coeff_index", [](const Row& x) + { + return matlab::output_index(x.max_coeff_index()); + }, + MATRIXBASE_ADDONS_VECTOR_INDEX_MAX_COEFF_INDEX_CONST) + + .def("__repr__", [](const Row& x) + { + std::ostringstream s; + s << x; + return string(s.str()); + }, + DOC_TO_BE_DEFINED) + + ; + + py::implicitly_convertible(); + + return exported_row_class; +} diff --git a/python/src/core/matrices/codac2_py_inversion.cpp b/python/src/core/matrices/codac2_py_inversion.cpp index 8c513e7e2..c81d84041 100644 --- a/python/src/core/matrices/codac2_py_inversion.cpp +++ b/python/src/core/matrices/codac2_py_inversion.cpp @@ -36,6 +36,22 @@ void export_inversion(py::module& m) INTERVALMATRIX_INVERSE_ENCLOSURE_CONST_EIGEN_MATRIXBASE_OTHERDERIVED_REF, "A"_a) + .def("left_inverse_enclosure", [](const Matrix& A) { return left_inverse_enclosure(A); }, + INTERVALMATRIX_LEFT_INVERSE_ENCLOSURE_CONST_EIGEN_MATRIXBASE_OTHERDERIVED_REF, + "A"_a) + + .def("left_inverse_enclosure", [](const IntervalMatrix& A) { return left_inverse_enclosure(A); }, + INTERVALMATRIX_LEFT_INVERSE_ENCLOSURE_CONST_EIGEN_MATRIXBASE_OTHERDERIVED_REF, + "A"_a) + + .def("right_inverse_enclosure", [](const Matrix& A) { return right_inverse_enclosure(A); }, + INTERVALMATRIX_RIGHT_INVERSE_ENCLOSURE_CONST_EIGEN_MATRIXBASE_OTHERDERIVED_REF, + "A"_a) + + .def("right_inverse_enclosure", [](const IntervalMatrix& A) { return right_inverse_enclosure(A); }, + INTERVALMATRIX_RIGHT_INVERSE_ENCLOSURE_CONST_EIGEN_MATRIXBASE_OTHERDERIVED_REF, + "A"_a) + .def("inverse_correction", [](const IntervalMatrix& A, const IntervalMatrix& B, bool left_inv = true) { diff --git a/python/src/core/separators/codac2_py_SepPolytope.cpp b/python/src/core/separators/codac2_py_SepPolytope.cpp new file mode 100644 index 000000000..8aa17f5ce --- /dev/null +++ b/python/src/core/separators/codac2_py_SepPolytope.cpp @@ -0,0 +1,37 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Simon Rohou, Damien Massé + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include "codac2_py_Sep.h" +#include "codac2_py_SepPolytope_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_SepPolytope(py::module& m, py::class_& pysep) +{ + py::class_ exported(m, "SepPolytope", pysep, SEPPOLYTOPE_MAIN); + exported + + .def(py::init(), + SEPPOLYTOPE_SEPPOLYTOPE_CONST_POLYTOPE_REF, + "p"_a) + + .def("separate", &SepPolytope::separate, + BOXPAIR_SEPPOLYTOPE_SEPARATE_CONST_INTERVALVECTOR_REF_CONST, + "x"_a) + ; + + py::implicitly_convertible(); +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 34511403e..68af50d6c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -115,6 +115,7 @@ endif() + file(APPEND ${CODAC_CMAKE_CONFIG_FILE} " set(CODAC_LIBRARIES \${CODAC_LIBRARIES} \${CODAC_GRAPHICS_LIBRARY} \${CODAC_CORE_LIBRARY}) ") @@ -133,4 +134,4 @@ get_filename_component(header_name ${header_path} NAME) file(APPEND ${CODAC_MAIN_HEADER} "#include <${header_name}>\n") endforeach() - install(FILES ${CODAC_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/) \ No newline at end of file + install(FILES ${CODAC_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/) diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 2c6a4172e..e87b86b50 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -45,6 +45,8 @@ ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcPolar.h ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcPolygon.cpp ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcPolygon.h + ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcPolytope.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcPolytope.h ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcProj.cpp ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcProj.h ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcSegment.cpp @@ -77,6 +79,14 @@ ${CMAKE_CURRENT_SOURCE_DIR}/domains/zonotope/codac2_Parallelepiped.cpp ${CMAKE_CURRENT_SOURCE_DIR}/domains/zonotope/codac2_Zonotope.h ${CMAKE_CURRENT_SOURCE_DIR}/domains/zonotope/codac2_Zonotope.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Facet.h + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Facet.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Polytope_util.h + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Polytope_util.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Polytope_dd.h + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Polytope_dd.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Polytope.h + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope/codac2_Polytope.cpp ${CMAKE_CURRENT_SOURCE_DIR}/domains/tube/codac2_Slice.h ${CMAKE_CURRENT_SOURCE_DIR}/domains/tube/codac2_SliceBase.cpp ${CMAKE_CURRENT_SOURCE_DIR}/domains/tube/codac2_SliceBase.h @@ -228,6 +238,8 @@ ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepNot.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolygon.cpp ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolygon.h + ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolytope.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepPolytope.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepProj.cpp ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepProj.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepTransform.cpp @@ -284,6 +296,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}/domains/interval/eigen ${CMAKE_CURRENT_SOURCE_DIR}/domains/paving ${CMAKE_CURRENT_SOURCE_DIR}/domains/zonotope + ${CMAKE_CURRENT_SOURCE_DIR}/domains/polytope ${CMAKE_CURRENT_SOURCE_DIR}/domains/tube ${CMAKE_CURRENT_SOURCE_DIR}/functions ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic diff --git a/src/core/contractors/codac2_CtcPolytope.cpp b/src/core/contractors/codac2_CtcPolytope.cpp new file mode 100644 index 000000000..51484eb40 --- /dev/null +++ b/src/core/contractors/codac2_CtcPolytope.cpp @@ -0,0 +1,31 @@ +/** + * codac2_CtcPolytope.cpp + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_CtcPolytope.h" + +using namespace std; +using namespace codac2; + + +CtcPolytope::CtcPolytope(const Polytope& p) + : Ctc(p.size()), _p(p) +{ + assert(p.size()>0); +} + +void CtcPolytope::contract(IntervalVector& x) const +{ + assert(x.size() == _p.size()); + x &= _p.box(); + if (x.is_empty()) return; + Polytope tmp (_p); + tmp.meet_with_box(x); + x &= tmp.box(true); +} + diff --git a/src/core/contractors/codac2_CtcPolytope.h b/src/core/contractors/codac2_CtcPolytope.h new file mode 100644 index 000000000..ed2430fe1 --- /dev/null +++ b/src/core/contractors/codac2_CtcPolytope.h @@ -0,0 +1,31 @@ +/** + * \file codac2_CtcPolytope.h + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include "codac2_Polytope.h" +#include "codac2_CtcWrapper.h" + +namespace codac2 +{ + class CtcPolytope : public Ctc + { + public: + + CtcPolytope(const Polytope& p); + + void contract(IntervalVector& x) const; + inline const Polytope& p() const { return _p; }; + + protected: + + const Polytope _p; + }; + +} diff --git a/src/core/domains/polytope/codac2_Facet.cpp b/src/core/domains/polytope/codac2_Facet.cpp new file mode 100644 index 000000000..f85e460f6 --- /dev/null +++ b/src/core/domains/polytope/codac2_Facet.cpp @@ -0,0 +1,542 @@ +/** + * \file codac2_Facet.cpp + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + + +#include /* uses swap */ +#include + +#include "codac2_Index.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_Interval.h" +#include "codac2_IntervalVector.h" +#include "codac2_IntervalRow.h" +#include "codac2_BoolInterval.h" +#include "codac2_Facet.h" + +namespace codac2 { + +namespace Facet_ { + +polytope_inclrel relation_Box(const Facet &f, const IntervalVector &b, bool strict) { + const Row &row = f.first.get_row(); + const double rhs = f.second.get_rhs(); + const bool eqcst = f.second.is_eqcst(); + if (eqcst && strict) return inclrel_notinclude | inclrel_disjoint; + if (b.is_empty()) return inclrel_includes | inclrel_disjoint; + IntervalVector a(b); + /* check the vertex that maximizes row */ + for (Index i=0;i0) a[i]=a[i].ub(); else a[i]=a[i].lb(); + } + Interval maxv = row.dot(a)-rhs; + polytope_inclrel r1=0; + if (maxv.ub()<=0.0 && (!strict || maxv.ub()<0.0)) { + if (!eqcst) return inclrel_includes | inclrel_intersects; + else if (maxv.ub()<0.0) return inclrel_notinclude | inclrel_disjoint; + r1 = inclrel_includes; + } else if (maxv.lb()<=0.0 && (!strict || maxv.lb()<0.0)) { + r1 = inclrel_mayinclude; + } else { + r1 = inclrel_notinclude; + } + /* check the vertex that minimizes row */ + for (Index i=0;i0) a[i]=b[i].lb(); else a[i]=b[i].ub(); + } + Interval minv = row.dot(a)-rhs; + if (minv.lb()>0.0 || (strict && minv.lb()>=0.0)) + return inclrel_notinclude | inclrel_disjoint; + if (!eqcst) { + if (minv.ub()>0.0 || (strict && minv.ub()>=0.0)) + return r1 | inclrel_mayintersect; + return r1 | inclrel_intersects; + } else { /* eqcst==true -> strict == false */ + if (r1[INCLUDES]) { /* maxv.ub == 0 */ + if (minv.lb()>=0.0) return + inclrel_includes | inclrel_intersects; + return (maxv.lb()==0.0 ? inclrel_intersects : inclrel_mayintersect) + | (minv.ub()<0.0 ? inclrel_notinclude : inclrel_mayinclude); + } + if (minv.ub()<0.0) r1 = inclrel_notinclude; + if (maxv.lb()>=0.0 && minv.ub()<=0.0) + return r1 | inclrel_intersects; + else return r1 | inclrel_mayintersect; + } +} + +void contract_Box(const Facet &f, IntervalVector &b) { + const Row &row = f.first.get_row(); + const double rhs = f.second.get_rhs(); + /* use MulOp:bwd */ + IntervalRow x1(row); + MulOp::bwd(f.second.is_eqcst() ? Interval(rhs,rhs) : Interval(-oo,rhs),x1,b); + if (b[0].is_empty()) b.set_empty(); +} + +void contract_out_Box(const Facet &f, IntervalVector &b) { + const Row &row = f.first.get_row(); + const double rhs = f.second.get_rhs(); + if (f.first.is_coord()) { + Index c = f.first.gt_dim(); + Interval val = Interval(rhs)/row[c]; + if (f.second.is_eqcst()) { + if (!val.is_degenerated()) return; /* can't use a value */ + if (val==b[c]) b.set_empty(); + } else { + if (row[c]>0.0) { /* x[c]>val.lb() */ + if (b[c].ub()<=val.lb()) { b.set_empty(); } + else if (b[c].lb()<=val.lb()) + { b[c] = Interval(val.lb(),b[c].ub()); } + } + } + return; + } + Interval a = b.dot(row); + if (a.ub()<=rhs) { b.set_empty(); return; } + if (f.second.is_eqcst()) { + if (a.lb()>=rhs) b.set_empty(); + return; + } + /* use MulOp:bwd */ + IntervalRow x1(row); + MulOp::bwd(Interval(rhs,oo),x1,b); + if (b[0].is_empty()) b.set_empty(); +} + +Interval bound_linear_form(const Facet &f, + const Row &row2, const IntervalVector &b) { + if (f.first.is_null()) return Interval(); + const Row &row = f.first.get_row(); + const double rhs = f.second.get_rhs(); + const bool eqcst = f.second.is_eqcst(); + const Index abdim = f.first.gt_dim(); + bool neg = (row[abdim]<0.0); + if (!eqcst && + ((!neg && row2[abdim]<=0.0) || (neg && row2[abdim]>=0))) + return Interval(); + Interval quotient = Interval(row2[abdim])/Interval(row[abdim]); + Interval res=quotient*rhs; + IntervalRow rem=row2-quotient*row; + rem[abdim]=0.0; + return res + rem.dot(b); +} + +} /* end of namespace Facet_ */ + +std::ostream& operator<<(std::ostream& os, const Facet& f) { + os << f.second.Id << "(" << f.first.bdim << "," << f.first.vdim << ") : " << f.first.row << (f.second.eqcst ? "=" : "<=" ) << f.second.rhs; + return os; +} + +void output_normalised_facet(std::ostream& os, const Facet& f) { + os << f.second.Id << " : " << (f.first.row/std::abs(f.first.row[f.first.gt_dim()])) << (f.second.eqcst ? "=" : "<=" ) << (f.second.rhs/std::abs(f.first.row[f.first.gt_dim()])); +} + +std::ostream& operator<<(std::ostream& os, const CollectFacets& cf) { + os << " Collectfacets : " << cf._map.size() << " facets" << std::endl; + for (const auto &f : cf._map) { + output_normalised_facet(os,f); + os << std::endl; + } + os << " end Collectfacets" << std::endl; + return os; +} + +CollectFacets::CollectFacets(const Matrix &mat, const Vector &rhsvect, + const std::vector &eqSet) : dim(mat.cols()) { + assert(mat.rows()==rhsvect.size()); + for (Index i=0;i a =this->insert_facet(mat.row(i),rhsvect[i],false); + if (!a.second) + throw std::domain_error("CollectFacets construction with duplicate rows"); + } + for (Index i: eqSet) { + _eqFacets.push_back(i); + (_allFacets[i])->second.eqcst=true; + } +} + +std::pair CollectFacets::result_insertion + (std::pair res, + double rhs, bool eqcst, ACTION_DUPLICATE act) { + if (res.second) { + Index id = _allFacets.size()+1; + res.first->second.Id=id; + _allFacets.push_back(res.first); + if (eqcst) _eqFacets.push_back(id-1); + return { id, true }; + } else { + if (act==KEEP_RHS) return { 0, false }; + double old_rhs=res.first->second.rhs; + bool old_eqcst = res.first->second.eqcst; + double new_rhs; + bool new_eqcst; + if (act==CHANGE_RHS) { new_eqcst = eqcst; new_rhs = rhs; } + else if (act==MAX_RHS) { + if (old_eqcst && eqcst) { + if (old_rhs==rhs) return { 0, false }; + /* complex case */ + res.first->second.rhs=std::max(old_rhs,rhs); + res.first->second.eqcst=false; + remove_in_eqFacets(res.first->second.Id-1); + std::pair rinsert = + this->insert_facet(-res.first->first.row, + -std::min(old_rhs,rhs),false,MAX_RHS); + /* can't loop, as we insert an inequality */ + return rinsert; + } else { + if (old_rhs>=rhs && !old_eqcst) return { 0, false }; + new_eqcst=false; + new_rhs=std::max(old_rhs,rhs); + } + } else { /* MIN_RHS */ + if (old_rhs<=rhs) { + if (!eqcst) return {0, false}; + if (old_rhssecond.rhs=new_rhs; + res.first->second.eqcst=new_eqcst; + if (new_eqcst && !old_eqcst) + _eqFacets.push_back(res.first->second.Id-1); + else if (!new_eqcst && old_eqcst) + remove_in_eqFacets(res.first->second.Id-1); + return { res.first->second.Id, false }; + } + return { 0, false }; +} + +void CollectFacets::remove_in_eqFacets(Index id) { + for (Index i=0; i<(Index)_eqFacets.size(); i++) { + if (_eqFacets[i]==id) { + _eqFacets[i]=_eqFacets.back(); + _eqFacets.pop_back(); + return; + } + } +} + +std::pair CollectFacets::insert_facet(const Row &row, + double rhs, bool eqcst, ACTION_DUPLICATE act) { + std::pair res = + _map.emplace(std::move(Facet_::make(row,rhs,eqcst,0))); + return result_insertion(res,rhs,eqcst,act); +} +std::pair CollectFacets::insert_facet(Row &&row, double rhs, bool eqcst, ACTION_DUPLICATE act) { + std::pair res = + _map.emplace(std::move(Facet_::make(row,rhs,eqcst,0))); + return result_insertion(res,rhs,eqcst,act); +} +std::pair CollectFacets::insert_facet(const Facet &fct, ACTION_DUPLICATE act) { + std::pair res = _map.insert(fct); + return result_insertion(res,fct.second.rhs,fct.second.eqcst,act); +} + +std::pair + CollectFacets::change_eqFacet(Index id, Row &&nrow, double nrhs, + ACTION_DUPLICATE act) { + mapIterator fcIt = _allFacets[_eqFacets[id]]; + auto nde = _map.extract(fcIt); + nde.key().change_row(nrow); + nde.mapped().rhs=nrhs; + const auto result_insert = _map.insert(std::move(nde)); + if (!result_insert.inserted) { + double old_rhs=result_insert.position->second.rhs; + if (act==CHANGE_RHS || + (act==MAX_RHS && old_rhsnrhs)) { + bool old_eqcst = result_insert.position->second.eqcst; + result_insert.position->second.rhs=nrhs; + result_insert.position->second.eqcst=true; + result_insert.position->second.Id=_eqFacets[id]+1; + _allFacets[_eqFacets[id]] = result_insert.position; + if (old_eqcst) + remove_in_eqFacets(result_insert.position->second.Id-1); + _allFacets[result_insert.position->second.Id-1] = endFacet(); + nb_removed_facets++; + return { result_insert.position, false }; + } else { + _allFacets[_eqFacets[id]] = endFacet(); + nb_removed_facets++; + _eqFacets[id] = _eqFacets.back(); + _eqFacets.pop_back(); + return { this->_map.end(), false }; + } + } + _allFacets[_eqFacets[id]] = result_insert.position; + return { result_insert.position, true }; +} + +std::pair + CollectFacets::change_ineqFacet(Index id, Row &&nrow, double nrhs, + ACTION_DUPLICATE act) { + mapIterator fcIt = _allFacets[id-1]; + auto nde = _map.extract(fcIt); + nde.key().change_row(nrow); + nde.mapped().rhs=nrhs; + const auto result_insert = _map.insert(std::move(nde)); + if (!result_insert.inserted) { + double old_rhs=result_insert.position->second.rhs; + if (act==CHANGE_RHS || + (act==MAX_RHS && old_rhsnrhs)) { + bool old_eqcst = result_insert.position->second.eqcst; + result_insert.position->second.rhs=nrhs; + result_insert.position->second.eqcst=true; + if (old_eqcst) + remove_in_eqFacets(result_insert.position->second.Id-1); + _allFacets[result_insert.position->second.Id-1] = endFacet(); + nb_removed_facets++; + result_insert.position->second.Id=id; + _allFacets[id-1] = result_insert.position; + return { result_insert.position, false }; + } else { + _allFacets[id-1] = endFacet(); + nb_removed_facets++; + return { this->_map.end(), false }; + } + } + _allFacets[id-1] = result_insert.position; + return { result_insert.position, true }; +} + +CollectFacets::mapIterator CollectFacets::change_ineqFacet_rhs(mapIterator it, + double nrhs) { + it->second.rhs=nrhs; + return it; +} + +CollectFacets::mapIterator + CollectFacets::dissociate_eqFacet(Index id, double nbound, + ACTION_DUPLICATE act) { + Index aId = _eqFacets[id]; + /* start with removing id in _eqFacets */ + _eqFacets[id] = _eqFacets.back(); + _eqFacets.pop_back(); + mapIterator fcIt = _allFacets[aId]; + mapIterator ret = this->_map.end(); + double rhs = fcIt->second.rhs; + if (nbound<=rhs) { + fcIt->second.eqcst=false; + if (nbound>-oo) { + std::pair rinsert = + this->insert_facet(-fcIt->first.row,-nbound,false,act); + if (rinsert.first!=0) ret = _allFacets[rinsert.first-1]; + } else ret = fcIt; + } else { + fcIt->second.eqcst=false; + if (nbound<+oo) { + std::pair rinsert = this->insert_facet(-fcIt->first.row, + -fcIt->second.rhs,false,act); + fcIt->second.rhs = nbound; + if (rinsert.first!=0) ret = _allFacets[rinsert.first-1]; + } else { /* we need to extract and reinsert the facet */ + auto nde = _map.extract(fcIt); + nde.key().negate_row(); + double nrhs = -nde.mapped().rhs; + nde.mapped().rhs= nrhs; + nde.mapped().eqcst= false; + const auto result_insert = _map.insert(std::move(nde)); + if (!result_insert.inserted) { + double old_rhs=result_insert.position->second.rhs; + /* CHECKME : not the best idea ? */ + if (act==CHANGE_RHS || + (act==MAX_RHS && old_rhsnrhs)) { + bool old_eqcst = result_insert.position->second.eqcst; + result_insert.position->second.rhs=nrhs; + result_insert.position->second.eqcst=false; + if (old_eqcst) + remove_in_eqFacets(result_insert.position->second.Id-1); + _allFacets[result_insert.position->second.Id-1] = endFacet(); + nb_removed_facets++; + result_insert.position->second.Id=aId+1; + _allFacets[aId] = result_insert.position; + ret=result_insert.position; + } else { + _allFacets[aId] = endFacet(); + nb_removed_facets++; + ret=endFacet(); + } + } else { + _allFacets[aId] = result_insert.position; + ret=result_insert.position; + } + } + } + return ret; +} + +bool CollectFacets::removeFacetById(Index id) { + if (_allFacets[id-1]==_map.end()) return false; + if (_allFacets[id-1]->second.eqcst) { + this->remove_in_eqFacets(id-1); + } + _map.erase(_allFacets[id-1]); + _allFacets[id-1]=endFacet(); + nb_removed_facets++; + return true; +} + +IntervalVector CollectFacets::extractBox() { + assert_release(this->getDim()!=-1); + /* check emptiness with null row */ + Row row = Row::zero(this->dim); + mapIterator it = _map.find(FacetBase(row)); + if (it!=_map.end()) { + FacetRhs &rhs = it->second; + if (rhs.rhs<0.0 || (rhs.rhs>0.0 && rhs.eqcst)) + return IntervalVector::constant(this->getDim(),Interval::empty()); + _allFacets[it->second.Id-1]= endFacet(); + nb_removed_facets++; + if (rhs.eqcst) remove_in_eqFacets(it->second.Id-1); + _map.erase(it); + } + /* check bounds if IV */ + IntervalVector ret = IntervalVector(this->getDim()); + for (Index i=0;igetDim();i++) { + Interval &bound = ret[i]; + std::pair pIt + = this->range_between_bases(FacetBase::base_range(this->getDim(),i, + false)); + while (pIt.first!=pIt.second && !bound.is_empty()) { + FacetRhs &rhs = pIt.first->second; + Interval val(pIt.first->first.row[i]); + if (rhs.eqcst) bound &= rhs.rhs/val; + else bound = min(bound,rhs.rhs/val); + _allFacets[pIt.first->second.Id-1]=endFacet(); + nb_removed_facets++; + if (rhs.eqcst) remove_in_eqFacets(pIt.first->second.Id-1); + pIt.first = _map.erase(pIt.first); /* pIt points to the next element */ + } + pIt = this->range_between_bases(FacetBase::base_range(this->getDim(),i, + true)); + while (pIt.first!=pIt.second && !bound.is_empty()) { + FacetRhs &rhs = pIt.first->second; + Interval val(pIt.first->first.row[i]); + if (rhs.eqcst) bound &= rhs.rhs/val; + else bound = max(bound,rhs.rhs/val); + _allFacets[pIt.first->second.Id-1]=endFacet(); + nb_removed_facets++; + if (rhs.eqcst) remove_in_eqFacets(pIt.first->second.Id-1); + pIt.first = _map.erase(pIt.first); /* pIt points to the next element */ + } + if (bound.is_empty()) { ret.set_empty(); break; } + } + return ret; +} + +std::vector CollectFacets::renumber() { + /* renumber the constraints, only if there are removed facets */ + if (nb_removed_facets==0) return std::vector(); + std::vector ret(_allFacets.size(),-1); + Index sineq=0, seq=0; + _allFacets.resize(_map.size()); + for (mapIterator it = _map.begin(); it!=_map.end(); ++it) { + if (it->second.eqcst) { + if (seq>=(Index)_eqFacets.size()) + _eqFacets.push_back(sineq); + else + _eqFacets[seq]=sineq; + ++seq; + } + _allFacets[sineq++]=it; + ret[it->second.Id-1]=sineq; + it->second.Id = sineq; + } + if (seq<(Index)_eqFacets.size()) + _eqFacets.resize(seq); + nb_removed_facets=0; + return ret; +} + +void CollectFacets::encompass_vertices(const + std::vector &vertices, + IntervalVector &bbox, bool tight) { + assert_release(this->getDim()!=-1); + if (tight) bbox=IntervalVector::constant(this->getDim(),Interval::empty()); + if (vertices.size()==0) return; + for (mapIterator it = _map.begin(); it!=_map.end(); ++it) { + IntervalRow row = it->first.row; + Interval a = (tight ? Interval::empty() : Interval(it->second.rhs)); + for (const IntervalVector &v : vertices) { + a |= row.dot(v); + } + it->second.rhs=a.ub(); + if (it->second.eqcst && !a.is_degenerated()) { + this->remove_in_eqFacets(it->second.Id-1); + it->second.eqcst=false; + this->insert_facet(-it->first.row,-a.lb(),false); + /* note : possible that the new facet will be visited + later, but that should not be a real problem */ + } + } + for (const IntervalVector &v : vertices) { + bbox |= v; + } +} +void CollectFacets::encompass_vertices(const std::vector &vertices, + IntervalVector &bbox, bool tight) { + assert_release(this->getDim()!=-1); + if (tight) bbox=IntervalVector::constant(this->getDim(),Interval::empty()); + if (vertices.size()==0) return; + for (mapIterator it = _map.begin(); it!=_map.end(); ++it) { + IntervalRow row = it->first.row; + Interval a = (tight ? Interval::empty() : Interval(it->second.rhs)); + for (const Vector &v : vertices) { + a |= row.dot(v); + } + it->second.rhs=a.ub(); + if (it->second.eqcst && !a.is_degenerated()) { + this->remove_in_eqFacets(it->second.Id-1); + it->second.eqcst=false; + this->insert_facet(-it->first.row,-a.lb(),false); + /* note : possible that the new facet will be visited + later, but that should not be a real problem */ + } + } + for (const Vector &v : vertices) { + bbox |= IntervalVector(v); + } +} + +void CollectFacets::encompass_zonotope(const IntervalVector &z, + const IntervalMatrix &A, + const IntervalVector &range, bool tight) { + assert_release(this->getDim()!=-1); + for (mapIterator it = _map.begin(); it!=_map.end(); ++it) { + Interval a = (tight ? Interval::empty() : Interval(it->second.rhs)); + a |= ((z.dot(it->first.row)) + (it->first.row * A).dot(range)); + it->second.rhs=a.ub(); + if (it->second.eqcst && !a.is_degenerated()) { + this->remove_in_eqFacets(it->second.Id-1); + it->second.eqcst=false; + this->insert_facet(-it->first.row,-a.lb(),false); + /* note : possible that the new facet will be visited + later, but that should not be a real problem */ + } + } +} + +int CollectFacets::merge_CollectFacets(const CollectFacets &collect, + ACTION_DUPLICATE act) { + int count=0; + for (const auto &fct : collect._map) { + std::pair r = this->insert_facet(fct,act); + if (r.first==-1) return -1; + if (r.second) count++; + } + return count; +} + +} diff --git a/src/core/domains/polytope/codac2_Facet.h b/src/core/domains/polytope/codac2_Facet.h new file mode 100644 index 000000000..b0f228704 --- /dev/null +++ b/src/core/domains/polytope/codac2_Facet.h @@ -0,0 +1,600 @@ +/** + * \file codac2_Facet.h + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#pragma once + +#include /* uses swap */ +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_Interval.h" +#include "codac2_IntervalVector.h" +#include "codac2_IntervalRow.h" +#include "codac2_Zonotope.h" +#include "codac2_BoolInterval.h" +#include "codac2_operators.h" + +namespace codac2 { + +class FacetBase; +class FacetRhs; +/** \brief a facet is a pair of FacetBase and FacetRhs */ +using Facet = std::pair; + +/** Encapsulation of Row (key of facet maps) */ +class FacetBase { + /* friends */ + friend class CollectFacets; + friend class Polytope; + friend class DDbuildF2V; + friend class DDbuildV2F; + friend struct FacetBaseComp; + + public: + + /** \brief build a new FacetBase from a row + * \param row the Row + */ + FacetBase(const Row &row) : bdim(-1), vdim(0.0), row(row) { + this->compute_key(); + } + + /** \brief build a new FacetBase from a row + * \param row the Row + */ + FacetBase(Row &&row) : bdim(-1), vdim(0.0), row(row) { + this->compute_key(); + } + + FacetBase(const FacetBase &fc) = default; + FacetBase(FacetBase &&fc) = default; + + /** \brief get the row + * \return the row */ + inline const Row &get_row() const { return row; } + + /** \brief test if the Row is zero + * \return true if row = 0 */ + inline bool is_null() const { return (bdim==-1); } + + /** \brief test if the row has only one non-zero coordinate + * \return true if row is zero for all coordinates except one */ + inline bool is_coord() const { + if (row.size()==0) return false; + return (bdim%(2*row.size())==0); + } + + /** \brief return the index of the greatest coordinate (in absolute value) + * suppose the row is not zero + * \return a value between 0 and size-1 */ + inline Index gt_dim () const { + if (row.size()==0) return -1; + return (bdim/(2*row.size()))%row.size(); + } + + inline void swap(FacetBase &fb) { + this->row.swap(fb.row); + std::swap(this->bdim,fb.bdim); + std::swap(this->bdim,fb.bdim); + } + + private: + + Index bdim; /* if A is the index of the greatest abs value and B + the second greatest (=A if 0 outside) + v[A]>0 v[B]>0 A*(2*dim)+A-B (+2*dim if B>A) + v[A]>0 v[B]<0 A*(2*dim)+A-(B+dim)+2*dim + v[A]<0 v[B]>0 (A+dim)*(2*dim)+A+dim-B + v[A]<0 v[B]<0 (A+dim)*(2*dim)+(A+dim)-(B+dim) + (+2*dim if B>A) + noting Adim = A if v[A]>0 and A+dim if v[A]<0 + and Bdim = B if v[B]>0 and B+dim if v[B]<0 + bdim = Adim*(2*dim+1) - Bdim (+2*dim if Adim>Bdim) + if null vector, value = -1 */ + double vdim; /* |v[B]|/|v[A]| (0 if no value for B) */ + Row row; + + + /* create a "min" and "max" FacetBase for constraints of the form + a_i x_i <= b_i (or - a_i x_i <= b_i if neg is true */ + static std::pair + base_range(Index dim, Index i, bool neg) { + Row r = Row::zero(dim); + Index d = (i+(neg?dim:0))*(2*dim); + return std::make_pair + (FacetBase(d,-1.0,r),FacetBase(d,1.0,r)); + } + + inline void negate_row() { + this->row=-row; + Index u = 2*row.size()*row.size(); + if (this->bdim>=u) this->bdim-=u; else this->bdim+=u; + } + + inline void change_row(Row &&row) { + this->row=row; this->bdim=-1; this->vdim=0.0; + this->compute_key(); + } + inline void change_row(const Row &row) { + this->row=row; this->bdim=-1; this->vdim=0.0; + this->compute_key(); + } + + FacetBase(Index bdim, double vdim, const Row &row) : + bdim(bdim), vdim(vdim), row(row) { } + + void compute_key() { + if (row.size()==0) return; + double val=row[0]; + double valabs=std::abs(val); + if (val<0.0) bdim=row.size(); else if (val>0.0) bdim=0; + Index b2dim = -1; + double val2abs=0.0; + for (Index i=1;ivalabs) { + val2abs=valabs; + valabs=val_i; + val=row[i]; + b2dim=bdim; + if (val<0.0) bdim=row.size()+i; else bdim=i; + } else if (val_i>val2abs) { + val2abs=val_i; + if (row[i]<0.0) b2dim=row.size()+i; else b2dim=i; + } + } + if (valabs==0.0) return; /* bdim=-1 */ + if (b2dim==-1) { bdim = 2*row.size()*bdim; vdim=0.0; return; } + if (b2dim>bdim) + bdim = bdim * (2*row.size()+1) + 2*row.size() - b2dim; + else + bdim = bdim * (2*row.size()+1) - b2dim; + vdim = val2abs/valabs; + } + friend std::ostream& operator<<(std::ostream& os, const Facet& f); + friend void output_normalised_facet(std::ostream& os, const Facet& f); + +}; + +/** Second part of the facet: nature of the facet (equality or not), + * right-hand side. An Id number is also kept for internal use inside + * a collection of facets. This number may change when the collection is + * modified */ +class FacetRhs { + /* friends */ + friend class CollectFacets; + friend class Polytope; + friend class DDbuildF2V; + friend class DDbuildV2F; + + public: + /** \brief value of the rhs + * \return the rhs */ + double get_rhs() const { return rhs; } + + /** \brief nature of the facet (= or <=) + * \return true if equality, false if inequality */ + bool is_eqcst() const { return eqcst; } + + /** \brief id number of the facet + * \return id number */ + Index get_Id() const { return Id; } + + FacetRhs(double rhs, bool eqcst, Index Id=0) : rhs(rhs), eqcst(eqcst), Id(Id) + {}; + FacetRhs(const FacetRhs &rhs) = default; + + void swap(FacetRhs &r) { + std::swap(rhs,r.rhs); + std::swap(eqcst,r.eqcst); + } + + private : + + double rhs; + bool eqcst; + Index Id; + friend std::ostream& operator<<(std::ostream& os, const Facet& f); + friend void output_normalised_facet(std::ostream& os, const Facet& f); + +}; + + +/** namespace Facet_ is used to encapsulate Facet functions */ +namespace Facet_ { + + inline Facet make(const Row &row, double rhs, bool eqcst, Index Id) { + return std::make_pair(FacetBase(row),FacetRhs(rhs,eqcst,Id)); + } + inline Facet make(Row &&row, double rhs, bool eqcst, Index Id) { + return std::make_pair(FacetBase(row),FacetRhs(rhs,eqcst,Id)); + } + + /** \brief create a new facet from a row, rhs, eqcst + * \param row the row + * \param rhs the rhs + * \param eqcst true for equalities, false for <= + * \return the new Facet */ + inline Facet make(const Row &row, double rhs, bool eqcst) { + return std::make_pair(FacetBase(row),FacetRhs(rhs,eqcst,0)); + } + + /** \brief create a new facet from a row, rhs, eqcst + * \param row the row + * \param rhs the rhs + * \param eqcst true for equalities, false for <= + * \return the new Facet */ + inline Facet make(Row &&row, double rhs, bool eqcst) { + return std::make_pair(FacetBase(row),FacetRhs(rhs,eqcst,0)); + } + + inline Facet make(const FacetBase &base, double rhs, bool eqcst) { + return std::make_pair(base,FacetRhs(rhs,eqcst,0)); + } + + inline void swap_Facet(Facet &f1,Facet &f2) { + f1.first.swap(f2.first); + f1.second.swap(f2.second); + } + + + /** \brief Potential inclusion relation between sets + * INCLUDES : A is included in B + * MAYINCLUDE : A may be included in B + * NOTINCLUDE : A\\B is not empty + * INTERSECTS : A inter B is non empty + * MAYINTERSECT : A inter B may be non empty + * DISJOINT : A inter B is empty + */ + enum INCLREL { + INCLUDES, + MAYINCLUDE, + NOTINCLUDE, + INTERSECTS, + MAYINTERSECT, + DISJOINT, + SIZEINCLREL + }; + using polytope_inclrel = std::bitset; + const polytope_inclrel inclrel_includes = 1<rhs.bdim) return false; + if (lhs.vdimrhs.vdim) return false; + for (int i=0;irhs.row[i]) return false; + return false; + } +}; + +/** encapsulate a collection of facets with double index : + by the base vector, and by Id number. Renumbering can happen + when facets are removed. */ +class CollectFacets { + public: + using mapType = std::map; + using mapIterator = mapType::iterator; + using mapCIterator = mapType::const_iterator; + + /** generate an empty collection, without dimension */ + CollectFacets() : dim(-1) {}; + /** generate an empty collection, with dimension + * @param dim the dimension + */ + CollectFacets(Index dim) : dim(dim) {}; + /** generate the set of facets from a matrix and a vector, i.e. + * mat X <= rhs (eqSet states which rows of mat are equalities) + * @param mat the matrix + * @param rhsvect the vector of rhs + * @param eqSet the set of equalities */ + CollectFacets(const Matrix &mat, const Vector &rhsvect, + const std::vector &eqSet); + + /** @brief create a copy of the CollectFacets, updating _allFacets + * @param cf the CollectFacets + */ + CollectFacets(const CollectFacets& cf); + + /** move the elements of CollectFacets. Perform a renumbering of Id if + * nb_removed_facets > 0 + * @param cf the CollectFacets + */ + CollectFacets(CollectFacets&& cf); + + ~CollectFacets() = default; + + /** return the dimension of the facets + * @return the dimension */ + Index getDim() const; + /** return the number of facets, as the size of _allfacets */ + Index nbfcts() const; + /** return the number of equality facets, as the size of _eqfacets */ + Index nbeqfcts() const; + /** return the map */ + const mapType &get_map() const; + + enum ACTION_DUPLICATE { /* operation to do when a modification + create a duplicate facet */ + KEEP_RHS, /* change nothing (default behaviour) */ + CHANGE_RHS, /* change the value in all cases */ + MAX_RHS, /* maximize the RHS. if two equalities, + duplicate the constraint. if one equality, + change into inequality */ + MIN_RHS /* minimize the RHS. if needed, return "empty" (-1) */ + }; + + + /** insert a new facet row X <= rhs (or row X = rhs) + * @param row Row + * @param rhs right-hand side + * @param eqcst equality constraint + * @param act action if an rhs with the same row exists + * @return the Id of the constraint, and true if a new facet + * has been inserted, false if old facet */ + std::pair insert_facet(const Row &row, + double rhs, bool eqcst, ACTION_DUPLICATE act=KEEP_RHS); + /** insert a new facet row X <= rhs (or row X = rhs) + * @param row Row + * @param rhs right-hand side + * @param eqcst equality constraint + * @param act action if an rhs with the same row exists + * @return the Id of the new constraint, and true if a new facet + * has been inserted, false if old facet. If old facet and no change, + * returns 0 */ + std::pair insert_facet(Row &&row, double rhs, bool eqcst, + ACTION_DUPLICATE act=KEEP_RHS); + /** insert a new facet row X <= rhs or row X = rhs, + * using the Facet constuction + * @param facet the facet + * @param act action if an rhs with the same row exists + * @return the Id of the constraint, and true if a new facet + * has been inserted, false if old facet. If old facet and no change + * returns 0 */ + std::pair insert_facet(const Facet &facet, + ACTION_DUPLICATE act=KEEP_RHS); + + /** get a facet by Id, starting from 0 (i.e. facet of index Id+1) + * @param id the Id of the facet, minus 1 + * @return an iterator on the facet */ + mapIterator operator[](Index id); + /** get a equality facet by the equality Id, + * starting from 0 (i.e. facet of index eqIndex(Id+1)) + * @param eqId the equality Id of the facet + * @return an iterator on the facet */ + mapIterator get_eqFacet(Index eqId); + /** get the end iterator of the facet + * @return an iterator on the facet */ + mapIterator endFacet(); + /** get the end iterator of the facet + * @return an iterator on the facet */ + mapCIterator endFacet() const; + /** get the first iterator + * @return an iterator on the facet */ + mapIterator beginFacet(); + /** get the first iterator + * @return an iterator on the facet */ + mapCIterator beginFacet() const; + + + /** change a eqfacet, keeping its Id (unless the new row creates + * a duplicate) + * @param eqId (the eqId of the facet) + * @param nrow (the new row of the facet) + * @param nrhs (the new right-hand side) + * @param act action if duplicate + * @return an iterator on the new facet, true if the modified facet + * has been inserted, false if there was a duplicate. In this case + * if the existing facet is not changed, return endFacet() */ + std::pair change_eqFacet(Index eqId, Row&& nrow, + double nrhs, ACTION_DUPLICATE act=KEEP_RHS); + /** change a facet, keeping its Id (unless the new row creates + * a duplicate) + * @param id the id of the facet starting from 1 + * @param nrow the new row of the facet + * @param nrhs the new right-hand side + * @param act action if duplicate + * @return an iterator on the new facet, true if the modified facet + * has been inserted, false if there was a duplicate. In this + * cas, without change, return endFacet() */ + std::pair change_ineqFacet(Index id, Row&& nrow, + double nrhs, ACTION_DUPLICATE act=KEEP_RHS); + + /** change the rhs of a facet + * @param it the iterator of the facet + * @param nrhs the new right-hand side + * @return the iterator of the facet (=it) */ + mapIterator change_ineqFacet_rhs(mapIterator it, double nrhs); + + /** transform an equality facet to an unequality facet, and add + a new facet with a new bound nbound + if nboundrhs, put row x <= nbound and add -row x <= -rhs + if nbound = -oo or +oo, do not modify/add the new facet + @param eqId the equality Id of the equality facet + @param nbound the new bound + @param act action if duplicate + @return depends on the construction + nbound=-oo => current facet iterator + -oo new facet iterator (for nbound), or endFacet if + insertion "failed" + rhs new facet iterator (for rhs), or endFacet if + insertion "failed" + nbound=+oo => new facet iterator (for rhs), or endFacet if + insertion "failed" */ + mapIterator dissociate_eqFacet(Index eqId , double nbound, + ACTION_DUPLICATE act=MAX_RHS); + + /** insert a complete collection of facets + @param collect the collection + @param act handling of the duplicate +  @return the number of added facets, -1 if empty + */ + int merge_CollectFacets(const CollectFacets &collect, + ACTION_DUPLICATE act=MIN_RHS); + + /** remove a facet by its Id + * @param id Id of the facet (starting from 1) + * @return true if the facet existed and has been removed */ + bool removeFacetById(Index id); + + /** extract the constraints corresponding to the bounding box, + * or a possible "emptiness constraint". + * remove some contraints, may put endFacet() for some Id */ + IntervalVector extractBox(); + + /** \brief renumber the set, if _allFacets have removed facets. + * do not modify the iterators, but modify the indices of the + * facets and equality facets. + * \return the matching old Id => new Id */ + std::vector renumber(); + + /** \brief modify the rhs such that the polyhedron encompasses a list + * of intervalVector. Also adapt the bbox. + * if tight=true, may tighten the bound, and not only increase them. + * may add facets if needed to dissociate equality facets. + * \param vertices the list of vertices + * \param bbox the bbox (to be updated) + * \param tight if false, the rhs are only increased */ + void encompass_vertices(const std::vector &vertices, + IntervalVector &bbox, bool tight); + void encompass_vertices(const std::vector &vertices, + IntervalVector &bbox, bool tight); + + /** \brief modify the rhs such that the polyhedron encompasses an + * "interval zonotope" z+A v. + * if tight=true, may tighten the bound, and not only increase them. + * may add facets if needed to dissociate equality facets. + * \brief z center of the interval zonotope + * \brief A matrix of the interval zonotope + * \brief v interval vector for the interval zonotope + * \brief tight if false, the rhs are only increased */ + void encompass_zonotope(const IntervalVector &z, const IntervalMatrix &A, + const IntervalVector &range, bool tight); + + /** \brief dump the set of facets, for debugging purposes */ + friend std::ostream& operator<<(std::ostream& os, + const CollectFacets& cf); + + + private: + mutable Index dim; + mapType _map; + std::vector _allFacets; + std::vector _eqFacets; /* index in _allFacets */ + int nb_removed_facets=0; /* number of removed facets in _allFacets */ + + /* look for and remove a value in eqFacet */ + void remove_in_eqFacets(Index id); + /* finish an insertion process */ + std::pair result_insertion(std::pair res, + double rhs, bool eqcst, ACTION_DUPLICATE act); + /* return a range of iterators between two bases */ + std::pair range_between_bases + (const std::pair &bases); + +}; + +inline CollectFacets::CollectFacets(const CollectFacets& cf) : + dim(cf.getDim()), _map(cf._map), + _allFacets(cf._allFacets.size(),_map.end()), + _eqFacets(cf._eqFacets), nb_removed_facets(cf.nb_removed_facets) +{ + /* redirect the facets */ + for (mapIterator it = _map.begin(); it!=_map.end(); ++it) { + _allFacets[it->second.Id-1]=it; + } +} + +inline CollectFacets::CollectFacets(CollectFacets&& cf) : + dim(cf.getDim()), _map(std::move(cf._map)), + _allFacets(std::move(cf._allFacets)), + _eqFacets(cf._eqFacets), nb_removed_facets(cf.nb_removed_facets) +{ + if (nb_removed_facets>0) this->renumber(); +} + +inline Index CollectFacets::getDim() const { + if (this->dim!=-1) return this->dim; + if (_map.size()==0) return -1; + this->dim=(_map.begin()->first.row.size()); + return this->dim; +} + +inline const CollectFacets::mapType &CollectFacets::get_map() const + { return _map; } + +inline CollectFacets::mapIterator CollectFacets::operator[](Index id) + { return this->_allFacets[id]; } +inline CollectFacets::mapIterator CollectFacets::get_eqFacet(Index id) + { return this->_allFacets[this->_eqFacets[id]]; } + +inline Index CollectFacets::nbfcts() const { + return this->_map.size(); +} + +inline Index CollectFacets::nbeqfcts() const { + return this->_eqFacets.size(); +} + +inline CollectFacets::mapIterator CollectFacets::endFacet() { + return this->_map.end(); +} +inline CollectFacets::mapCIterator CollectFacets::endFacet() const { + return this->_map.cend(); +} +inline CollectFacets::mapIterator CollectFacets::beginFacet() { + return this->_map.begin(); +} +inline CollectFacets::mapCIterator CollectFacets::beginFacet() const { + return this->_map.cbegin(); +} + + +inline std::pair + CollectFacets::range_between_bases + (const std::pair &bases) { + return std::pair + (_map.lower_bound(bases.first),_map.lower_bound(bases.second)); +} + +} diff --git a/src/core/domains/polytope/codac2_Polytope.cpp b/src/core/domains/polytope/codac2_Polytope.cpp new file mode 100644 index 000000000..d19af892d --- /dev/null +++ b/src/core/domains/polytope/codac2_Polytope.cpp @@ -0,0 +1,1249 @@ +/** + * \file codac2_Polytope.cpp + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Matrix.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_IntervalRow.h" +#include "codac2_IntervalVector.h" +#include "codac2_IntervalMatrix.h" +#include "codac2_inversion.h" +#include "codac2_Parallelepiped.h" +#include "codac2_Zonotope.h" +#include "codac2_Polytope.h" +#include "codac2_Facet.h" + +using namespace codac2; + + +namespace codac2 { + +Polytope::Polytope() : _dim(-1), state(pol_state_empty | pol_state(1<(dim)), + state(empty ? pol_state_empty : pol_state_init) +{ +} + +Polytope::Polytope(const IntervalVector &box) : _dim(box.size()), _box(box), + _facets(std::make_shared(_dim)), + state(box.is_empty() ? pol_state_empty : pol_state_init) + +{ +} + +Polytope::Polytope(const Polytope &P) : _dim(P._dim), _box(P._box), + _facets(P.state[EMPTY] ? std::make_shared(_dim) : + std::make_shared(*(P._facets))), + state(P.state & (pol_state_empty | + pol_state(1<_dim = P._dim; + this->_box = P._box; + this->_facets = P.state[EMPTY] ? std::make_shared(_dim) : + std::make_shared(*(P._facets)); + this->_DDbuildF2V=nullptr; + this->_DDbuildV2F=nullptr; + this->state = P.state & (pol_state_empty | + pol_state(1< &vertices) : + Polytope() +{ + if (vertices.empty()) return; + _dim=vertices[0].size(); + /* build the V2F form */ + _DDbuildV2F = std::make_unique(1,vertices[0]); + for (Index i=1;i<(Index)vertices.size();i++) { + _DDbuildV2F->add_vertex(i+1,vertices[i]); + } + /* get the CollectFacets */ + _facets=_DDbuildV2F->getFacets(); + _box = _facets->extractBox(); + _facets->encompass_vertices(vertices,_box,true); + state = pol_state_init; + _facets->renumber(); + _DDbuildV2F.reset(); /* encompass may have added constraints */ +} + +Polytope::Polytope(const std::vector &vertices) : + Polytope() +{ + _dim=vertices[0].size(); + /* build the V2F form */ + _DDbuildV2F = std::make_unique(1,vertices[0].mid()); + for (Index i=1;i<(Index)vertices.size();i++) { + _DDbuildV2F->add_vertex(i+1,vertices[i].mid()); + } + /* get the CollectFacets */ + _facets=_DDbuildV2F->getFacets(); + _box = _facets->extractBox(); + _facets->encompass_vertices(vertices,_box,true); + _facets->renumber(); + state = pol_state_init; + _DDbuildV2F.reset(); /* encompass may have added constraints */ +} + +Polytope::Polytope(const std::vector &vertices, + const CollectFacets &facetsform) : Polytope() +{ + _dim = facetsform.getDim(); + if (_dim==-1) return; + _facets=std::make_shared(facetsform); + _box = IntervalVector::constant(_dim,Interval::empty()); + if (vertices.empty()) { + state = pol_state_empty; + return; + } + for (Index i=0;i<(Index)vertices.size();i++) { + _box |= vertices[i]; + } + _facets->encompass_vertices(vertices,_box,true); + state = pol_state_init; +} + + +Polytope::Polytope(const IntervalVector &box, + const std::vector> &facets, + bool minimize) : Polytope(box) { + for (const std::pair &cst : facets) { + this->add_constraint(cst); + } + state = 0; + if (minimize) this->minimize_constraints(); +} + + + +Polytope::Polytope(const Parallelepiped &par) : + Polytope(par.box()) { + /* we can use the vertices, or inverse the shape, + we consider the inversion */ + IntervalMatrix u = inverse_enclosure(par.A); + for (Index i=0;iadd_constraint_band(IntervalRow(u.row(i)), + u.row(i).dot(par.z)+Interval(-1.0,1.0),0.0); + } + IntervalVector range = IntervalVector::constant(par.A.cols(), + Interval(-1.0,1.0)); + _box = par.z+par.A*range; /* FIXME : use box */ +} + +Polytope::Polytope(const Zonotope &zon) : + Polytope() { + _dim=zon.z.size(); + /* we use the vertices (note : can do better) */ + IntervalVector v = IntervalVector::constant(zon.A.cols(),-1.0); + this->_DDbuildV2F = std::make_unique(1,(zon.z+zon.A*v).mid()); + Index idV = 2; + Index a = 0; + while (a_DDbuildV2F->add_vertex(idV++,IntervalVector(zon.z+zon.A*v)); + } + } + /* get the CollectFacets */ + _facets=_DDbuildV2F->getFacets(); + IntervalVector range = IntervalVector::constant(zon.A.cols(), + Interval(-1.0,1.0)); + _box = zon.z+zon.A*range; /* FIXME : use box */ + _facets->encompass_zonotope(IntervalVector(zon.z), + IntervalMatrix(zon.A), range, true); + _DDbuildV2F.reset(); /* encompass may have added constraints */ + state = pol_state_init; +} + +Polytope::Polytope(CollectFacets &&facets) : _dim(facets.getDim()), state(0) { + _box = facets.extractBox(); + _facets = std::make_shared(std::move(facets)); + /* include renumbering if needed */ + state = 0; +} + +Polytope::Polytope(IntervalVector &&box, + CollectFacets &&facets) : _dim(facets.getDim()), + _box(std::move(box)), state(0) { + _box &= facets.extractBox(); + _facets = std::make_shared(std::move(facets)); +} + + +Polytope Polytope::from_ineFile(const char *filename) { + std::shared_ptr fcts + = std::make_shared(read_ineFile(filename)); + if (fcts->getDim()==-1) return Polytope(); + Polytope p(fcts->getDim(),false); + p._facets = fcts; + p.state = 0; + return p; +} + +Polytope Polytope::from_extFile(const char *filename) { + std::vector vts = read_extFile(filename); + return Polytope(vts); +} + + +Interval Polytope::fast_bound(const FacetBase &base) const { + if (base.is_null()) return Interval::zero(); + if (state[EMPTY]) return Interval(); + if (is_box()) return base.row.dot(_box).ub(); + Index gdim = base.gt_dim(); + Interval res; + if (base.row[gdim]>0.0) { + res= Interval(_box[gdim].ub())*base.row[gdim]; + } else { + res= Interval(_box[gdim].lb())*base.row[gdim]; + } + if (!base.is_coord()) { + Row bcopy = base.row; + bcopy[gdim]=0.0; + res += bcopy.dot(_box); + if (_facets->get_map().empty()) return res; + std::pair + pIt = _facets->get_map().equal_range(base); + if (pIt.first!=pIt.second) { /* we found an exact key */ + return Interval(pIt.first->second.rhs); + } else { + if (pIt.first!=_facets->get_map().begin()) --pIt.first; + Interval a1 = Facet_::bound_linear_form(*(pIt.first), base.row, _box); + if (pIt.second!=_facets->endFacet()) { + Interval a2 = Facet_::bound_linear_form(*(pIt.second), base.row, _box); + if (a2.ub()nbeqfcts();i++) { + CollectFacets::mapIterator it = _facets->get_eqFacet(i); + Interval a = Facet_::bound_linear_form(*it, base.row, _box); + if (a.ub() &vertices=_DDbuildF2V->get_vertices(); + for (const DDvertex &vt : vertices) { + IntervalVector vtB = _DDbuildF2V->compute_vertex(vt.vertex); + Interval a = base.row.dot(vtB); + if (a.ub()build_DDbuildF2V(); + double a = -oo; + if (state[EMPTY]) return a; + const std::forward_list &vertices=_DDbuildF2V->get_vertices(); + for (const DDvertex &vt : vertices) { + IntervalVector vtB = _DDbuildF2V->compute_vertex(vt.vertex); + Interval b = r.dot(vtB); + a = std::max(a,b.ub()); + } + return a; +} + +void Polytope::build_DDbuildF2V() const { + if (state[F2VFORM]) return; + if (state[EMPTY]) { + _DDbuildF2V=nullptr; + return; + } + /* TODO : insert box progressively ? */ + _DDbuildF2V = std::make_unique(_dim,_box,_facets,true); + for (CollectFacets::mapCIterator it = _facets->get_map().begin(); + it!=_facets->get_map().end(); ++it) { + if ((*it).second.eqcst) continue; + _DDbuildF2V->add_facet(it); + if (_DDbuildF2V->is_empty()) { + this->set_empty_private(); + return; + } + state[NOTEMPTY]=true; + } + state[F2VFORM]=true; +} + +void Polytope::minimize_constraints() const { + if (state[MINIMIZED]) return; + if (state[EMPTY]) return; + this->build_DDbuildF2V(); + if (state[EMPTY]) return; + std::set redundant = _DDbuildF2V->redundantFacets(); + _box &= _DDbuildF2V->build_bbox(); + state[BOXUPDATED]=true; + bool changed=false; + for (const Index &a : redundant) { + if (a<=0) continue; + _facets->removeFacetById(a); + changed=true; + } + if (changed) { + std::vector corresp = _facets->renumber(); + _DDbuildF2V->update_renumber(corresp); + state[V2FFORM]=false; + } + state[MINIMIZED]=true; +} + +void Polytope::update_box() const { + if (state[BOXUPDATED]) return; + if (state[EMPTY]) return; + if (is_box()) { state[BOXUPDATED]=true; return; } + assert_release(!state[BOXUPDATED]); + this->build_DDbuildF2V(); + if (state[EMPTY]) return; + _box &= _DDbuildF2V->build_bbox(); + state[BOXUPDATED]=true; +} + +bool Polytope::check_empty() const { + if (state[NOTEMPTY]) return false; + if (state[EMPTY]) return true; + if (is_box()) { + if (_box.is_empty()) { + state[EMPTY]=true; + return true; + } else { + state[NOTEMPTY]=true; + return false; + } + } + this->build_DDbuildF2V(); + return state[EMPTY]; +} + +BoolInterval Polytope::contains(const IntervalVector& p) const { + assert_release(p.size()==_dim); + if (p.is_empty()) return BoolInterval::TRUE; + if (state[EMPTY]) return BoolInterval::FALSE; + if (!p.is_subset(_box)) return BoolInterval::FALSE; + BoolInterval r = BoolInterval::TRUE; + for (const auto &fct : _facets->get_map()) { + Facet_::polytope_inclrel pincl = + Facet_::relation_Box(fct,p); + if (pincl[Facet_::NOTINCLUDE]) return BoolInterval::FALSE; + if (pincl[Facet_::MAYINCLUDE]) r = BoolInterval::UNKNOWN; + } + return r; +} + +BoolInterval Polytope::interior_contains(const IntervalVector& p) const { + assert_release(p.size()==_dim); + if (p.is_empty()) return BoolInterval::TRUE; + if (state[EMPTY]) return BoolInterval::FALSE; + if (!p.is_interior_subset(_box)) return BoolInterval::FALSE; + BoolInterval r = BoolInterval::TRUE; + for (const auto &fct : _facets->get_map()) { + Facet_::polytope_inclrel pincl = + Facet_::relation_Box(fct,p,true); + if (pincl[Facet_::NOTINCLUDE]) return BoolInterval::FALSE; + if (pincl[Facet_::MAYINCLUDE]) r = BoolInterval::UNKNOWN; + } + return r; +} + +bool Polytope::box_is_subset(const IntervalVector& x) const { + return this->box().is_subset(x); +} + +bool Polytope::is_subset(const Polytope& P) + const { + const IntervalVector &b2 = P.box(true); + const IntervalVector &b1 = this->box(true); + if (!b1.is_subset(b2)) return false; + for (const auto &fctP : P._facets->get_map()) { + if (this->fast_bound(fctP.first).ub() <= fctP.second.rhs) continue; + double l1 = this->bound_row(fctP.first.row); + if (l1<=fctP.second.rhs) continue; + return false; + } + return true; +} + +bool Polytope::is_superset(const Polytope& P) + const { + return P.is_subset((*this)); +} + +bool Polytope::is_interior_subset(const Polytope& P) + const { + const IntervalVector &b2 = P.box(true); + const IntervalVector &b1 = this->box(true); + if (!b1.is_interior_subset(b2)) return false; + for (const auto &fctP : P._facets->get_map()) { + if (this->fast_bound(fctP.first).ub() < fctP.second.rhs) continue; + double l1 = this->bound_row(fctP.first.row); + if (l1box(true); + if (!b1.intersects(p)) return false; + Polytope q(*this); /* TODO : copy DDbuildF2V */ + q.meet_with_box(p); + return !(q.check_empty()); +} + +bool Polytope::intersects(const Polytope& p) const { + assert(!status[INVALID]); + const IntervalVector &b1 = this->box(); + if (!b1.intersects(p.box())) return false; + Polytope q(*this); /* TODO : copy DDbuildF2V */ + q.meet_with_polytope(p); + return !(q.check_empty()); +} + +bool Polytope::is_disjoint(const Polytope& p) const { + return (!this->intersects(p)); +} + +void Polytope::contract_out_Box(IntervalVector &b) const { + assert(!status[INVALID]); + const IntervalVector &b1 = this->_box; + /* we can compute the diff, but as the goal is to return one interval vector, + we do it "by hand" to avoid the generation of too many elements. + also, we quit as soon a the vector is "full" */ + IntervalVector result = IntervalVector::constant(_dim,Interval::empty()); + for (Index i=0;i<_dim;i++) { + std::vector vt = b[i].diff(b1[i]); + if (vt.empty()) continue; + IntervalVector r1(b); + r1[i] = vt[0]; + if (vt.size()>1) r1[i] |= vt[1]; + if (r1[i]==b[i]) return; + result |= r1; + if (b.is_subset(result)) return; + } + for (auto &facet : _facets->get_map()) { + IntervalVector r1(b); + Facet_::contract_out_Box(facet,r1); + result |= r1; + if (b.is_subset(result)) return; + } + b = result; +} + +#if 0 + +polytope_inclrel Polytope::relation_Box(const IntervalVector& p) const { + assert_release(_dim>=0); + if (p.is_empty()) return (1<& facets = _clpForm->getFacets(); + for (const Facet &facet : facets) { + polytope_inclrel res2 = facet.relation_Box(p); + if (res2[DISJOINT]) return (1<relation_Box(p); + return (r[INCLUDES] ? BoolInterval::TRUE : + (r[NOTINCLUDE] ? BoolInterval::FALSE : BoolInterval::UNKNOWN)); +} + +BoolInterval Polytope::intersects(const IntervalVector& p) const { + polytope_inclrel r = this->relation_Box(p); + return (r[INTERSECTS] ? BoolInterval::TRUE : + (r[DISJOINT] ? BoolInterval::FALSE : BoolInterval::UNKNOWN)); +} + +const Interval& Polytope::operator[](Index i) const { + assert_release(i<=_dim); + return _box[i]; +} + +double Polytope::;istance_cst(const Facet &fc) const { + assert_release(_dim==fc.row.size()); + if (_empty) return -oo; + if (!_clpForm) return (fc.row.dot(_box)-fc.rhs).ub(); + _clpForm->setObjective(fc.row); + LPclp::lp_result_stat stat = _clpForm->solve(false,0); + if (stat[LPclp::EMPTY]) return -oo; + return (_clpForm->getValobj()-fc.rhs).ub(); +} + +#endif + +Vector Polytope::mid() const { + assert_release(!state[INVALID]); + this->update_box(); + return (_box.mid()); +} + +Vector Polytope::mid_in() const { + assert_release(!state[INVALID]); + if (this->is_box()) return (_box.mid()); + std::vector vt = this->vertices(); + if (vt.empty()) return _box.mid(); + Vector ret = vt[0].mid(); + for (Index i=1;i<(Index)vt.size();i++) + ret += vt[i].mid(); + return ret/vt.size(); +} + +void Polytope::clear() { + assert_release(!state[INVALID]); + state=pol_state_init; + _box = IntervalVector::Zero(_dim); + _facets=nullptr; + _DDbuildF2V=nullptr; + _DDbuildV2F=nullptr; +} + + +bool Polytope::add_constraint(const std::pair& facet, + double tolerance) { + if (state[EMPTY]) return false; + FacetBase base = FacetBase(facet.first); + if (base.is_null()) { + if (facet.second>=tolerance) return false; + this->set_empty(); return true; + } + Interval act = this->fast_bound(base); + if (act.ub()<=facet.second+tolerance) return false; + if (base.is_coord()) { + Index gdim = base.gt_dim(); + Interval val = Interval(facet.second)/facet.first[gdim]; + if (facet.first[gdim]>0.0) { + _box[gdim] &= Interval(-oo,val.ub()); + if (_box[gdim].is_empty()) { + this->set_empty_private(); + return true; + } +#if 0 + if (state[CLPFORM]) { + this->_clpForm->set_bbox(_box); + } +#endif + if (state[F2VFORM]) { + int ret = _DDbuildF2V->add_bound_var(gdim,true,val.ub()); + if (ret==-1) { + this->set_empty_private(); + return true; + } + if (ret==0) return false; + } + } else { + _box[gdim] &= Interval(-val.ub(),+oo); + if (_box[gdim].is_empty()) { + this->set_empty_private(); + return true; + } +#if 0 + if (state[CLPFORM]) { + this->_clpForm->set_bbox(_box); + } +#endif + if (state[F2VFORM]) { + int ret = _DDbuildF2V->add_bound_var(gdim,false,-val.ub()); + if (ret==-1) { + this->set_empty_private(); + return true; + } + if (ret==0) return false; + } + } + state[V2FFORM]=false; + state[MINIMIZED]=false; + state[NOTEMPTY]=false; + return true; + } + auto res= _facets->insert_facet(facet.first, facet.second, + false, CollectFacets::MIN_RHS); + if (res.first==-1) { this->set_empty_private(); return true; } + if (res.first==0) return false; +#if 0 + if (state[CLPFORM]) { + this->_clpForm->updateConstraint(res.first); + } +#endif + if (state[F2VFORM]) { + int ret = this->_DDbuildF2V->add_facet((*_facets)[res.first-1]); + if (ret==-1) { + this->set_empty_private(); + return true; + } + if (ret==0) return false; + } + state[V2FFORM]=state[MINIMIZED]=state[BOXUPDATED]=false; + state[NOTEMPTY]=false; + return true; +} + +bool Polytope::add_constraint(const IntervalRow &cst, double rhs, + double tolerance) { + if (state[EMPTY]) return false; + Row cstmid = cst.mid(); + IntervalRow rem = cst-cstmid; + Interval d = rhs+rem.dot(_box); + if (d.ub()==+oo) return false; + return this->add_constraint(std::pair(cstmid, d.ub()), tolerance); +} + +std::pair Polytope::add_constraint_band(const IntervalRow &cst, + const Interval &rhs, double tolerance) { + if (state[EMPTY]) return { false, false }; + Row cstmid = cst.mid(); + IntervalRow rem = cst-cstmid; + Interval d = rhs+rem.dot(_box); + return this->add_constraint_band(cstmid,d,tolerance); +} + +std::pair Polytope::add_constraint_band(const Row &cst, + const Interval &rhs, double tolerance) { + if (state[EMPTY]) return { false, false }; + bool rub, rlb; + if (rhs.ub()==+oo) rub=false; + else rub = this->add_constraint(std::pair(cst, rhs.ub()), tolerance); + if (rhs.lb()==-oo) rlb=false; + else rlb = this->add_constraint(std::pair(-cst, -rhs.lb()), tolerance); + return { rlb, rub }; +} + +int Polytope::add_equality(const std::pair& facet) { + if (state[EMPTY]) return -1; + FacetBase base(facet.first); + if (base.is_null()) { + if (facet.second==0.0) return 0; + this->set_empty(); return -1; + } + if (base.is_coord()) { + Index gdim = base.gt_dim(); + Interval val = Interval(facet.second)/facet.first[gdim]; + _box[gdim] &= val; + if (_box[gdim].is_empty()) { + this->set_empty_private(); + return -1; + } +#if 0 + if (state[CLPFORM]) { + this->_clpForm->set_bbox(_box); + } +#endif + state[F2VFORM]=false; + state[V2FFORM]=false; + state[MINIMIZED]=false; + state[NOTEMPTY]=false; + return 1; + } + FacetBase negBase(base); + negBase.negate_row(); + Interval actub = this->fast_bound(base); + Interval actlb = this->fast_bound(negBase); + if (actub.ub()set_empty(); return -1; } + if (actlb.ub()<-facet.second) { this->set_empty(); return -1; } + auto res= _facets->insert_facet(facet.first, facet.second, + true, CollectFacets::MIN_RHS); + if (res.first==-1) { this->set_empty_private(); return -1; } + if (res.first==0) return 0; +#if 0 + if (state[CLPFORM]) { + this->_clpForm->updateConstraint(res.first); + } +#endif + state[F2VFORM]=false; + state[V2FFORM]=state[MINIMIZED]=state[BOXUPDATED]=false; + state[NOTEMPTY]=false; + return 1; +} + +int Polytope::meet_with_box(const IntervalVector &b) { + assert(!state[INVALID]); + if (_box.is_subset(b)) return 0; + if (_box.is_disjoint(b)) { this->set_empty(); return -1; } + _box &= b; +#if 0 + if (state[CLPFORM]) { + this->_clpForm->set_bbox(_box); + } +#endif + if (state[F2VFORM]) { + if (!b.is_degenerated()) { + int ret = this->_DDbuildF2V->add_constraint_box(b); + if (ret==-1) { + this->set_empty(); + return -1; + } + if (ret==0) return 0; + } else { /* we do not keep the representation */ + state[F2VFORM]=false; + _DDbuildV2F=nullptr; + } + } + state[V2FFORM]=state[MINIMIZED]=state[BOXUPDATED]=false; + state[NOTEMPTY]=false; + return 1; +} + +Polytope &Polytope::operator&=(const IntervalVector &b) { + this->meet_with_box(b); + return (*this); +} + +Polytope Polytope::meet_with_hyperplane(Index dm, double x) const { + assert(!state[INVALID]); + if (!_box[dm].contains(x)) { return Polytope(this->_dim,true); } + IntervalVector nbox(_box); + nbox[dm]=x; + CollectFacets cf(_dim); + for (auto &facet : _facets->get_map()) { + Row r(facet.first.row); + Interval rhs = facet.second.rhs - nbox[dm]*r[dm]; + r[dm]=0.0; + if (!facet.second.eqcst || rhs.is_degenerated()) { + std::pair ret = + cf.insert_facet(std::move(r),rhs.ub(),facet.second.eqcst, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(this->_dim, true); + } else { + std::pair ret = + cf.insert_facet(r,rhs.ub(),false, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(this->_dim, true); + ret = cf.insert_facet(-r,-rhs.lb(),false, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(this->_dim, true); + } + } + return Polytope(std::move(nbox),std::move(cf)); +} + +int Polytope::meet_with_polytope(const Polytope &P) { + assert_release(!P.state[INVALID]); + assert_release(!state[INVALID]); + if (state[EMPTY]) return -1; + if (P.state[EMPTY]) { this->set_empty(); return -1; } + if (P._facets->nbfcts()==0) { return this->meet_with_box(P._box); } + int ret=0; + if (!_box.is_subset(P._box)) { + _box &= P._box; + ret=1; + if (_box.is_empty()) { this->set_empty(); return -1; } + } + int ret2= this->_facets->merge_CollectFacets(*(P._facets)); + if (ret2==-1) { this->set_empty(); return -1; } + if (ret2==0 && ret==0) return 0; + state[V2FFORM]=state[MINIMIZED]=state[BOXUPDATED]=false; + state[NOTEMPTY]=false; + state[F2VFORM]=false; +#if 0 + state[CLPFORM]=false; +#endif + return 1; +} + +Polytope &Polytope::homothety(const IntervalVector &c, double delta) { + assert_release(!state[INVALID]); + assert(delta>0); + Interval deltaI (delta); /* for interval computation */ + if (state[EMPTY]) return (*this); + _box = (1.0-deltaI) * c + delta * _box; +#if 0 + if (state[CLPFORM]) { + this->_clpForm->set_bbox(_box); + } +#endif + for (CollectFacets::mapIterator it = _facets->beginFacet() ; + it != _facets->endFacet(); ++it) { + if (it->second.eqcst) continue; + Interval nrhs = deltaI*it->second.rhs + + (1.0-deltaI)*it->first.row.dot(c); + _facets->change_ineqFacet_rhs(it,nrhs.ub()); + } + for (Index i=_facets->nbeqfcts()-1;i>=0;i--) { + CollectFacets::mapIterator it = _facets->get_eqFacet(i); + Interval nrhs = deltaI*it->second.rhs + + (1.0-deltaI)*it->first.row.dot(c); + _facets->change_ineqFacet_rhs(it,nrhs.ub()); + if (!nrhs.is_degenerated()) { + _facets->dissociate_eqFacet(i,nrhs.lb()); + } + } + state[V2FFORM]=false; + state[MINIMIZED]=false; + state[F2VFORM]=false; +#if 0 + state[CLPFORM]=false; /* TODO : possible if c is punctual */ +#endif + return (*this); +} + +Polytope Polytope::reverse_affine_transform(const IntervalMatrix &M, + const IntervalVector &P, const IntervalVector &bbox) const { + assert_release(!state[INVALID]); + assert_release(M.rows()==_dim); + assert_release(bbox.size()==M.cols()); + assert_release(P.size()==_dim); + if (state[EMPTY]) return Polytope(M.cols(),true); + CollectFacets cf(M.cols()); + /* first the box */ + for (Index i=0;i<_dim;i++) { + IntervalRow rI = M.row(i); + Interval rhs = _box[i] - P[i]; + Row r = rI.mid(); + rI = rI - r; + rhs -= rI.dot(bbox); + if (rhs.is_degenerated()) { + if (rhs.ub()==+oo) continue; + std::pair ret = + cf.insert_facet(std::move(r),rhs.ub(),true, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(M.cols(), true); + } else { + if (rhs.ub()!=+oo) { + std::pair ret = + cf.insert_facet(r,rhs.ub(),false, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(M.cols(), true); + } + if (rhs.lb()!=-oo) { + std::pair ret = + cf.insert_facet(-r,-rhs.lb(),false, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(M.cols(), true); + } + } + } + for (auto &facet : _facets->get_map()) { + IntervalRow rI = facet.first.row * M; + Interval rhs = facet.second.rhs - facet.first.row.dot(P); + Row r = rI.mid(); + rI = rI - r; + rhs -= rI.dot(bbox); + if (!facet.second.eqcst || rhs.is_degenerated()) { + if (rhs.ub()==+oo) continue; + std::pair ret = + cf.insert_facet(std::move(r),rhs.ub(),facet.second.eqcst, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(M.cols(), true); + } else { + if (rhs.ub()!=+oo) { + std::pair ret = + cf.insert_facet(r,rhs.ub(),false, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(M.cols(), true); + } + if (rhs.lb()!=-oo) { + std::pair ret = + cf.insert_facet(-r,-rhs.lb(),false, + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(M.cols(), true); + } + } + } + return Polytope(IntervalVector(bbox),std::move(cf)); +} + +Polytope Polytope::affine_transform(const IntervalMatrix &M, + const IntervalVector &P, const IntervalVector &B, + bool use_direct) const { + assert_release(!state[INVALID]); + if (state[EMPTY]) return Polytope(_dim,true); + assert_release(M.cols()==_dim); + assert_release(M.rows()==P.size()); + assert_release(B.size()==P.size()); + /* x' \in [M] x + [P] + becomes x' \in Mmid x + ([P] + ([M]-Mmid) B) */ + Matrix Mmid = M.mid(); + IntervalVector P2 = P + (M - Mmid)*_box; + if (Mmid.rows()==_dim) { + IntervalMatrix Minv = inverse_enclosure(Mmid); + if (!Minv.is_unbounded()) { + IntervalVector M2 = B & (M*_box+P); + if (!use_direct) + return this->reverse_affine_transform(Minv,-Minv*P2,M2); + else { + Polytope R = this->reverse_affine_transform(Minv,-Minv*P2,M2); + return (R &= this->direct_affine_transform(M,P)); + } + } + } else if (M.rows()>_dim) { + IntvFullPivLU LUdec(Matrix(Mmid.transpose())); + /* we use the tranpose for solve as we did not define + leftSolve */ + if (LUdec.is_surjective()==BoolInterval::TRUE) { + IntervalVector M2 = B & (M*_box+P); + IntervalMatrix inv = + LUdec.solve(IntervalMatrix::Identity(_dim,_dim)); + IntervalMatrix leftInv = inv.transpose(); + Polytope R = this->reverse_affine_transform(leftInv,-leftInv*P2,M2); + IntervalMatrix leftKer = LUdec.kernel().transpose(); + for (Index k = 0; k(r.mid(),u.mid())); + else R.add_constraint_band(r,u); + } + if (!use_direct) + return R; + else + return (R &= this->direct_affine_transform(M,P)); + } + } + Polytope R = this->direct_affine_transform(M,P); + R.meet_with_box(B); + return R; +} + +Polytope Polytope::direct_affine_transform(const IntervalMatrix &M, + const IntervalVector &P) const { + assert_release(!state[INVALID]); + if (state[EMPTY]) return Polytope(_dim,true); + this->build_DDbuildF2V(); + IntervalVector nBox = M*this->box()+P; + std::vector resultat; + if (state[EMPTY]) return Polytope(*this); + const std::forward_list &vertices=_DDbuildF2V->get_vertices(); + for (const DDvertex &vt : vertices) { + IntervalVector a = M*_DDbuildF2V->compute_vertex(vt.vertex)+P; + resultat.push_back(a); + } + Polytope res(resultat); + return (res &= nBox); +} + +Polytope Polytope::time_elapse(const IntervalVector &P, + const Interval &T) const { + this->build_DDbuildF2V(); + IntervalVector nBox = this->box()+T*P; + std::vector resultat; + if (state[EMPTY]) return Polytope(*this); + const std::forward_list &vertices=_DDbuildF2V->get_vertices(); + for (const DDvertex &vt : vertices) { + IntervalVector vert = _DDbuildF2V->compute_vertex(vt.vertex); + resultat.push_back(vert+T.lb()*P); + resultat.push_back(vert+T.ub()*P); + } + Polytope res(resultat); + return (res &= nBox); +} + +Polytope &Polytope::operator&=(const Polytope &P) { + this->meet_with_polytope(P); + return (*this); +} + +Polytope Polytope::union_of_polytopes(std::initializer_list lst) { + if (lst.size()==0) return Polytope(); + if (lst.size()==1) return Polytope(*(lst.begin())); + Index dim = lst.begin()->dim(); + IntervalVector box = IntervalVector::constant(dim,Interval::empty()); + std::vector lvert; + for (auto &P : lst) { + std::vector tmp = + P.vertices(); + lvert.insert(lvert.end(),tmp.begin(),tmp.end()); + box |= P.box(); + } + if (lvert.size()==0) { + return Polytope(dim,true); + } + Polytope ret= Polytope(lvert); + return (ret &= box); +} + +Polytope Polytope::union_of_polytopes(const std::vector &lst) { + if (lst.size()==0) return Polytope(); + if (lst.size()==1) return Polytope(*(lst.begin())); + Index dim = lst.begin()->dim(); + IntervalVector box = IntervalVector::constant(dim,Interval::empty()); + std::vector lvert; + for (auto &P : lst) { + std::vector tmp = + P.vertices(); + lvert.insert(lvert.end(),tmp.begin(),tmp.end()); + box |= P.box(); + } + if (lvert.size()==0) { + return Polytope(dim,true); + } + Polytope ret= Polytope(lvert); + return (ret &= box); +} + +int Polytope::join_with_polytope(const Polytope &P) { + assert(!state[INVALID]); + assert(!P.state[INVALID]); + if (P.state[EMPTY]) return 0; + if (state[EMPTY]) { *this=P; return 1; } + *this=union_of_polytopes( { *this, P } ); + return 1; +} + +Polytope &Polytope::operator|=(const Polytope &P) { + assert(!state[INVALID]); + assert(!P.state[INVALID]); + if (P.state[EMPTY]) return *this; + if (state[EMPTY]) { *this=P; return *this; } + *this=union_of_polytopes( { *this, P } ); + return *this; +} + +Polytope &Polytope::inflate(double rad) { + assert(!state[INVALID]); + return (this->inflate(IntervalVector::constant(_dim,Interval(-rad,rad)))); +} + +Polytope &Polytope::inflate_ball(double rad) { + assert(!state[INVALID]); + if (state[EMPTY]) return *this; + if (rad<=0.0) return *this; + _box.inflate(rad); + /* first, we consider the inequality facets */ + for (CollectFacets::mapIterator it = _facets->beginFacet() ; + it != _facets->endFacet(); ++it) { + if (it->second.eqcst) continue; + double nrm = (it->first.row).norm()*rad; + double nrhs = it->second.rhs+nrm; + _facets->change_ineqFacet_rhs(it,nrhs); + } + /* now, we consider the equality facets */ + /* decreasing order because we "destroy" the eqFacet array */ + for (Index i=_facets->nbeqfcts()-1;i>=0;i--) { + CollectFacets::mapIterator it = _facets->get_eqFacet(i); + double nrm = (it->first.row).norm()*rad; + Interval rhs2 = it->second.rhs+Interval(-nrm,nrm); + _facets->change_ineqFacet_rhs(it,rhs2.ub()); + _facets->dissociate_eqFacet(i,rhs2.lb()); + } + /* TODO : update */ + state[V2FFORM]=false; + state[F2VFORM]=false; + _DDbuildF2V=nullptr; + _DDbuildV2F=nullptr; +#if 0 + state[CLPFORM]=false; + _clpForm=nullptr; +#endif + return *this; +} + +Polytope &Polytope::inflate(const IntervalVector& box) { + assert(!state[INVALID]); + if (state[EMPTY]) return *this; + if (box.is_empty()) { this->set_empty(); return (*this); } + _box += box; + /* first, we consider the inequality facets */ + for (CollectFacets::mapIterator it = _facets->beginFacet() ; + it != _facets->endFacet(); ++it) { + if (it->second.eqcst) continue; + double nrhs = (it->second.rhs+box.dot(it->first.row)).ub(); + _facets->change_ineqFacet_rhs(it,nrhs); + } + /* now, we consider the equality facets */ + /* decreasing order because we "destroy" the eqFacet array + however, some facets are kept as equalities if there is no + expansion */ + for (Index i=_facets->nbeqfcts()-1;i>=0;i--) { + CollectFacets::mapIterator it = _facets->get_eqFacet(i); + Interval rhs2 = it->second.rhs+box.dot(it->first.row); + _facets->change_ineqFacet_rhs(it,rhs2.ub()); + if (!rhs2.is_degenerated()) + _facets->dissociate_eqFacet(i,rhs2.lb()); + } + /* TODO : update */ + state[V2FFORM]=false; + state[F2VFORM]=false; + _DDbuildF2V=nullptr; + _DDbuildV2F=nullptr; +#if 0 + _clpForm=nullptr; + state[CLPFORM]=false; +#endif + return (*this); +} + + +Polytope &Polytope::unflat(Index dm, double rad) { + assert(!state[INVALID]); + if (state[EMPTY]) return *this; + if (rad<=0.0) { return *this; } + _box[dm].inflate(rad); + /* first, we consider the inequality facets */ + Interval radI(-rad,rad); /* for interval computations */ + for (CollectFacets::mapIterator it = _facets->beginFacet() ; + it != _facets->endFacet(); ++it) { + if (it->second.eqcst) continue; + if (it->first.row[dm]==0.0) continue; + double nrhs = (it->second.rhs+radI*it->first.row[dm]).ub(); + _facets->change_ineqFacet_rhs(it,nrhs); + } + /* now, we consider the equality facets */ + /* decreasing order because we "destroy" the eqFacet array + however, some facets are kept as equalities if there is no + expansion */ + for (Index i=_facets->nbeqfcts()-1;i>=0;i--) { + CollectFacets::mapIterator it = _facets->get_eqFacet(i); + if (it->first.row[dm]==0.0) continue; + Interval rhs2 = it->second.rhs+radI*it->first.row[dm]; + _facets->change_ineqFacet_rhs(it,rhs2.ub()); + if (!rhs2.is_degenerated()) + _facets->dissociate_eqFacet(i,rhs2.lb()); + } + /* TODO : update */ + state[V2FFORM]=false; + state[F2VFORM]=false; + _DDbuildF2V=nullptr; + _DDbuildV2F=nullptr; +#if 0 + state[CLPFORM]=false; + _clpForm=nullptr; +#endif + return (*this); +} + +Polytope operator+ (const Polytope &p1, const Polytope &p2) { + assert(!p1.state[INVALID]); + assert(!p2.state[INVALID]); + if (p1.is_box() && p2.is_box()) return Polytope(p1.box()+p2.box()); + std::vector vt1 = p1.vertices(); + std::vector vt2 = p2.vertices(); + IntervalVector bres = p1.box() + p2.box(); + std::vector result; + for (auto &v1 : vt1) + for (auto &v2 : vt2) + result.push_back(v1+v2); + Polytope ret = Polytope(result); + return (ret &= bres); +} + +Polytope Polytope::operator-() const { + assert_release(!state[INVALID]); + if (state[EMPTY]) return Polytope(_dim,true); + CollectFacets cf(_dim); + IntervalVector bbox(-_box); + for (auto &facet : _facets->get_map()) { + FacetBase nf = facet.first; + nf.negate_row(); + std::pair ret = + cf.insert_facet(Facet_::make(nf,-facet.second.rhs, + facet.second.eqcst), + CollectFacets::MIN_RHS); + if (ret.first==-1) return Polytope(_dim, true); /* should not happen */ + } + Polytope ret = Polytope(IntervalVector(bbox),std::move(cf)); + ret.state=this->state & pol_state_init; + return ret; +} + +std::ostream& operator<<(std::ostream& os, + const Polytope &P) { + if (P.state[Polytope::EMPTY]) { + os << "Polytope(" /* << P.state */ << " empty dim " << P._dim << ")" << std::endl; + return os; + } + os << "Polytope(" /* << P.state */ << " bbox " << P._box << ") : " << std::endl; + os << *P._facets; + os << "EndPolytope" << std::endl; + return os; +} + +std::vector Polytope::vertices() const { + assert(!state[INVALID]); + if (this->is_box()) { + if (_box.is_empty()) return std::vector(); + std::vector ret; + Vector a = _box.lb(); + ret.push_back(IntervalVector(a)); + Index i=0; + while (i<_box.size()) { + if (_box[i].is_degenerated()) { i++; continue; } + if (a[i]==_box[i].lb()) { + a[i]=_box[i].ub(); + ret.push_back(IntervalVector(a)); + i=0; continue; + } + a[i]=_box[i].lb(); i++; + } + return ret; + } + this->build_DDbuildF2V(); + if (state[EMPTY]) return std::vector(); + std::vector ret; + for (const DDvertex &vtx : _DDbuildF2V->get_vertices()) { + ret.push_back(_DDbuildF2V->compute_vertex(vtx.vertex)); + } + return ret; +} + +std::vector> Polytope::vertices_3Dfacets() const { + this->build_DDbuildF2V(); + if (state[EMPTY]) return std::vector>(); + return build_3Dfacets(*_DDbuildF2V); +} + +std::vector Polytope::vertices_2Dfacet() const { + this->build_DDbuildF2V(); + if (state[EMPTY]) return std::vector(); + return build_2Dfacet(*_DDbuildF2V); +} + +} + diff --git a/src/core/domains/polytope/codac2_Polytope.h b/src/core/domains/polytope/codac2_Polytope.h new file mode 100644 index 000000000..cc3d8cc44 --- /dev/null +++ b/src/core/domains/polytope/codac2_Polytope.h @@ -0,0 +1,684 @@ +/** + * \file codac2_Polytope.h + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Matrix.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_IntervalRow.h" +#include "codac2_IntervalVector.h" +#include "codac2_IntervalMatrix.h" +#include "codac2_BoolInterval.h" +#include "codac2_Polytope_util.h" +#include "codac2_Facet.h" +#include "codac2_Polytope_dd.h" +#include "codac2_Parallelepiped.h" +#include "codac2_Zonotope.h" + +namespace codac2 { + /** + * \class Polytope + * Represents a bounded convex polytope as a set of constraints + * and a bounding box (the bounding box is part of the constraints) */ + class Polytope + { + public: + /** + * basic constuctor + */ + Polytope(); + + /** + * basic constuctor with a dimension, full + */ + explicit Polytope(Index dim); + + /** + * basic constuctor with a dimension, empty or full + */ + explicit Polytope(Index dim, bool empty); + + /** + * Constructs a bounded polytope representing a box + * from a bounded interval vector + */ + explicit Polytope(const IntervalVector &box); + + /** + * Copy (no minimization) + */ + Polytope(const Polytope &P); + +// Polytope(Polytope &&P) = default; + + /** + * Constructs a polytope from a set of vertices + * \param vertices the vertices + */ + explicit Polytope(const std::vector &vertices); + + /** + * Constructs a polytope from a set of vertices as intervalboxes + * \param vertices the vertices + */ + explicit Polytope(const std::vector &vertices); + + /** + * Constructs a polytope from a set of vertices (IntervalVector) + * and a set of linear forms, described as a CollectFacets + * \param vertices the vertices + * \param facetsform the facets description (the CollectFacets) + */ + Polytope(const std::vector &vertices, + const CollectFacets &facetsform); + + /** + * from a parallelepiped + * \param par the parallelepiped + */ + explicit Polytope(const Parallelepiped &par); + + /** + * from a zonotope + * \param zon the zonotope + */ + explicit Polytope(const Zonotope &zon); + + /** + * Constructs a polytope from a box and a CollectFacets which + * will be moved in the polytope. Perform a "box extraction" and a + * renumbering on the collection. + * \param box the box + * \param facets the collection of facets + */ + Polytope(IntervalVector &&box, CollectFacets &&facets); + + /** + * Constructs a polytope from a CollectFacets which + * will be moved in the polytope. Perform a "box extraction" and a + * renumbering on the collection. + * \param facets the collection of facets + */ + Polytope(CollectFacets &&facets); + + /** + * Constructs a polytope from a bounding box and + * a set of constraints (inequalities) + * + */ + Polytope(const IntervalVector &box, + const std::vector> &facets, + bool minimize=false); + + /** + * build a polytope from an ine File (cddlib format) + * \param filename the file */ + static Polytope from_ineFile(const char *filename); + /** + * build a polytope from an ext File (cddlib format) + * \param filename the file */ + static Polytope from_extFile(const char *filename); + + /** convex union of polytopes + * use vertices to join the polytopes, as such the result may be + * too large + * \param lst a non-empty initializer_list of valid polytopes of + * same dimension + * \return the polytope */ + static Polytope union_of_polytopes(std::initializer_list lst); + + /** convex union of polytopes + * use vertices to join the polytopes, as such the result may be + * too large + * \param lst a non-empty initializer_list of valid polytopes of + * same dimension + * \return the polytope */ + static Polytope union_of_polytopes(const std::vector &lst); + + /** + * Destructor + */ + virtual ~Polytope() = default; + + /** + * copy assignment operator + * \param P copy + */ + Polytope &operator=(const Polytope &P); + + Polytope &operator=(Polytope &&P) = default; + + + /***************** ACCESS **************************/ + /** + * Get the dimension of the space + * + * \return the dimension of the space + */ + Index dim() const; + Index size() const; + + /** + * Get the dimension of the polytope (ie the dimension + * of the space minus the equalities). -1 means that the polytope + * is empty. + * + * \return the dimension of the polytope. + */ +// Index dim_polytope() const; + + /** + * Returns true if the polytope is a box (i.e. there are no + * other facets). + * + * \return true if the polytope is a box, false otherwise. + */ + bool is_box() const; + + /** + * Get the number of constraints + * + * \return the number of facets (including equalities) + */ + Index nb_facets() const; + + /** + * Get the number of equality constraints + * + * \return the number of equalities + */ + Index nb_eq_facets() const; + + /** + * is empty + */ + bool is_empty(bool check=true) const; + + /** + * has a flat dimension + */ + bool is_flat() const; + + /** ``component'' : interval of the bounding box + * (only const !!!). Update the bbox if needed. + * \return enclosing of the component + */ + const Interval& operator[](Index i) const; + + /** + * ``middle'' of the bbox + * \return the middle of the box + */ + Vector mid() const; + + /** + * point inside the polytope using the vertices + * \return a barycenter of the middle of the vertices + */ + Vector mid_in() const; + + /** + * bound a constraint (fast), either by a "close" one (+box) + * or by the vertices if they are computed. + * no update or LP are done + * return an interval I : row X is guaranteed to be less than I.ub() + * and I.lb gives an indication of the precision of the computation + * (i.e. if diam(I) is low, a constraint close to row was used + * \param base the constraint, expressed as a FacetBase + * \return the interval I */ + Interval fast_bound(const FacetBase &base) const; + + /** + * bound a constraint, using DD + * \param r the row + * \return the max + */ + virtual double bound_row(const Row &r) const; + +#if 0 + /** ``distance'' from the satisfaction of a constraint + * ie the inflation of the rhs needed to satisfy the constraint + * \param fc the constraint + * \return an interval + */ + double distance_cst (const Facet &fc) const; + + /** + * relationships with a box (fast check) + * \param p the box + * \return polytope_inclrel checking inclusion and intersection */ + Facet_::polytope_inclrel relation_Box(const IntervalVector& p) const; +#endif + + /** + * Checks whether the polytope contains a given point, or + * includes a box + * + * \param p The point or box to check, enclosed in an IntervalVector. + * \return BoolInterval indicating possible containment. + */ + BoolInterval contains(const IntervalVector& p) const; + + /** + * Checks whether the interior of the polytope contains a given + * point, or includes a box + * + * \param p The point or box to check, enclosed in an IntervalVector. + * \return BoolInterval indicating possible containment. + */ + BoolInterval interior_contains(const IntervalVector& p) const; + + /** + * Checks inclusion in another polytope + * \param P the polytope + * \return true if inclusion is guaranteed + */ + virtual bool is_subset(const Polytope &P) const; + + /** + * Checks enclosure of another polytope + * \param P the polytope + * \return true if (*this) encloses P + */ + bool is_superset(const Polytope &P) const; + + /** + * Checks inclusion in the interior of another polytope + * \param P the polytope + * \return true if inclusion is guaranteed + */ + virtual bool is_interior_subset(const Polytope &P) const; + + /** intersects a box + * \param x the box (IntervalVector) + * \return false if the intersection is guaranteed to be empty + */ + bool intersects(const IntervalVector& x) const; + + /** intersects a polytope + * \param p the polytope + * \return false if the intersection is guaranteed to be empty + */ + bool intersects(const Polytope &p) const; + + /** test if it is disjoint with a polytope + * \param p the polytope + * \return true if the intersection is guaranteed to be empty + */ + bool is_disjoint(const Polytope &p) const; + + /** compute the box hull of B diff (*this) + * \param b a box, contracted. + */ + void contract_out_Box(IntervalVector &b) const; + + /** minimize the constraints, removing (possibly) redundant + * constraints. + */ + virtual void minimize_constraints() const; + + + /************* Box access *************************/ + + /** + * get the current bounding box of the polytope + * (note: it may not be tightest bounding box, it is the + * bounding box used for approximations results) + * + * \return The current bounding box + */ + const IntervalVector &box(bool tight=true) const; + + /** + * test if bounding box is included + * \return true if the polytope is included in x + */ + bool box_is_subset(const IntervalVector& x) const; + + /** + * tests if the polytope is bisectable + * (i.e. its bounding box is bisectable). + * do not update the box + * \return true if the bbox is bisectable + */ + bool is_bisectable() const; + + + /************* Modification **********************/ + /* keeping the current polytope */ + + /** set to empty */ + void set_empty(); + + /** set to (singleton) 0^dim */ + virtual void clear(); + + /** add a inequality (pair row X <= rhs) + * do basic checks, but do not minimize the system + * \param facet the constraint + * \param tolerance tolerance for redundancy checking (CLP only) + * \return false if the (basic) checks showed the constraint to + * be redundant (or inside tolerance), true if added */ + bool add_constraint(const std::pair& facet, + double tolerance=0.0); + + /** add a inequality with intervalVector (cst X <= rhs) + * using the bounding box + * do basic checks, but do not minimize the system + * \param cst the row constraint + *  \param rhs the rhs + * \param tolerance tolerance for redundancy checking (CLP only) + * \return false if the (basic) checks showed the constraint to + * be redundant (or inside tolerance), true if added */ + bool add_constraint(const IntervalRow &cst, double rhs, + double tolerance=0.0); + + /** two inequalities with intervalVector (cst X in rhs) + * using the bounding box + * do basic checks, but do not minimize the system. Also, + * should not be used for equalities (cst and rhs punctual) + * \param cst the row constraint + *  \param rhs the rhs + * \param tolerance tolerance for redundancy checking (CLP only) + * \return pair of booleans (one for each constraints possibly added */ + std::pair add_constraint_band(const IntervalRow &cst, + const Interval &rhs, double tolerance=0.0); + std::pair add_constraint_band(const Row &cst, + const Interval &rhs, double tolerance=0.0); + + /** add a equality (pair row X = rhs) + * do basic checks, but do not minimize the system + * \param facet the constraint + * \return -1 if empty, 0 no change (not probable), 1 change */ + int add_equality(const std::pair& facet); + + /** intersect with a box. + * \param b the box + * \return 0 if nothing done, 1 changed made, -1 results is empty */ + int meet_with_box(const IntervalVector &b); + + /** intersect with a box. + * \param b the box + * \return *this */ + Polytope &operator&=(const IntervalVector &b); + + /** intersect with a polytope. + * \param P the polytope + * \return 0 if nothing done, 1 changed made, -1 results is empty */ + int meet_with_polytope(const Polytope &P); + + /** intersect with a polytope. + * \param P the polytope + * \return *this */ + Polytope &operator&=(const Polytope &P); + + /** union with a polytope + * \param P the polytope + * \return 0 if nothing done, 1 changed made, -1 results is empty */ + int join_with_polytope(const Polytope &P); + + /** union with a polytope + * \param P the polytope + * \return *this */ + Polytope &operator|=(const Polytope &P); + + /** inflation by a cube, keeping the shape (not optimal) + * this <- this + [-rad,rad]^d + * \param rad radius of the box + * \return *this (keep minimized, except rounding) + */ + Polytope &inflate(double rad); + + /** inflation with a box, keeping the shape (not optimal) + * this <- this + box + * \param box the box + * \return *this (keep minimized, except rounding) + */ + Polytope &inflate(const IntervalVector& box); + + /** inflation by a ball, keeping the shape (not optimal) + * this <- this + B_d(rad) (note: vector norm is computed on double) + * \param rad radius of the ball + * \return *this (keep minimized, except rounding) + */ + Polytope &inflate_ball(double rad); + + /** expansion of a dimension (not optimal) + * equivalent to inflation by a box 0,...,[-rad,rad],...0 + * \param dm index of the dimension + * \param rad radius + * \return *this (keep minimized, except rounding) */ + Polytope &unflat(Index dm, double rad); + + /** centered homothety (optimal if c is punctual) + * x <- [c] + delta*(x-[c]) ( or (1-delta)[c] + delta*x ) + * note : delta must be > 0 + * \param c ``center'' + * \param delta expansion + * \return *this + */ + Polytope &homothety(const IntervalVector &c, double delta); + + /* creating a new polytope */ + + /** return a polytope intersected with an hyperplan + * in a specific coordinate (rebuild the constraints without + * any dependence on this coordinate in a new polytope) + * \param dm index of the coordinate + * \param x the value of the coordinate + * \return the polytope */ + Polytope meet_with_hyperplane(Index dm, double x) const; + + /** reverse affine transformation + * compute { x' in Box s.t. ([M] x' + [P]) in the initial polytope } + * The bounding box is needed to handle interval in [M]. + * \param M non-empty matrix + * \param P non-empty vector + * \param bbox bouding-box of the new polytope + * \return a new polytope */ + Polytope reverse_affine_transform(const IntervalMatrix &M, + const IntervalVector &P, const IntervalVector &bbox) const; + + /** affine transformation + * compute { [B] cap ([M] x + [P]) s.t. x in the initial polytope } + * If [M] is bijective, use inverse to compute the inverse of M + * If [M] is full-rank and injective, use the + * LU decomposition of transpose(M) to compute + * the left inverse of M and the kernel of M. + * Otherwise, call direct-affine-transform. + * if use_direct is true, use both methods if possible + * \param M matrix + * \param P non-empty vector + * \param B box limiting the result + * \param use_direct use direct_affine_transform in addition of + * the other method + * \return a new polytope */ + Polytope affine_transform(const IntervalMatrix &M, + const IntervalVector &P, const IntervalVector &B, + bool use_direct=false) const; + + /** affine transformation + * compute { [M] x + [P] s.t x in the initial polytope } + * Generate IntervalVector from the vertices, but not the + * the vertices of the IntervalVector themselves (hence not + * optimal if [M] and [P] are not punctual). As it uses vertices, + * it may be expensive. Use affine_transform to avoid use of + * vertices (for injective transformations). + * \param M non-empty matrix + * \param P non-empty vector + * \return a new polytope */ + Polytope direct_affine_transform(const IntervalMatrix &M, + const IntervalVector &P) const; + + /** time-elapse transformation + * compute { x + t [P] s.t x in the initial polytope and t in T } + * Generate IntervalVector from the vertices, but not the + * the vertices of the IntervalVector themselves (hence not + * optimal if [P] is not punctual). As it uses vertices, + * it may be expensive. An inexpensive but less optimal approach is + * to use inflate(T*P). + * \param P non-empty vector + * \param T non-empty interval + * \return a new polytope */ + Polytope time_elapse(const IntervalVector &P, const Interval &T) const; + + /** sum of two polytopes, based on the sum of vertices + * (expensive computations) + * \param p1 first polytope + * \param p2 second polytope + * \return a new polytope */ + friend Polytope operator+ (const Polytope &p1, const Polytope &p2); + + /** unary minus + * \return the negation of the current polytope */ + Polytope operator- () const; + + + /*********** Printing and other ``global access'' *********/ + /** + * output the bbox and the set of facets + * \param os the output stream + * \param P the polytope + * \return the stream + */ + friend std::ostream& operator<<(std::ostream& os, + const Polytope &P); + + /** + * Return the set of facets, as a map of pairs + * FacetBase/FacetRhs. Note: calling any method + * of the polytope (even const methods) can change the map. The + * bounding box is not included in the map (use box()). For + * empty polytopes, the set is empty. + * + * \return the facets + */ + const CollectFacets::mapType &facets() const; + + /** + * Computes a set of vertices enclosing the polytope + * + * \return set of vertices, as IntervalVectors + */ + std::vector vertices() const; + + /** + * Computes the set of 3D facets + * + * \return set of set of vertices, as Vectors + */ + std::vector> vertices_3Dfacets() const; + + /** + * Computes the vertices as a 2D facet + * + * \return set of set of vertices, as Vectors + * (mid of ``real'' vertices'') + */ + std::vector vertices_2Dfacet() const; + + + + + protected: + Index _dim; /* dimension */ + mutable IntervalVector _box; /* bounding box */ + mutable std::shared_ptr _facets; + /* "native" formulation , may be shared by other formulations */ + mutable std::unique_ptr _DDbuildF2V; + /* DDbuildF2V formulation, if used */ + mutable std::unique_ptr _DDbuildV2F; + /* DDbuildV2F formulation, generally not used */ + + /** state of the polytope */ + enum POLSTATE { + EMPTY, /* is empty */ + NOTEMPTY, /* is NOT empty */ + MINIMIZED, /* minimized : redundant constraints removed */ + BOXUPDATED, /* box is minimal */ + F2VFORM, /* has an up-to-date F2V form */ + V2FFORM, /* has an up-to-date V2F form */ + INVALID, /* not initialised */ + SIZE_POLSTATE + }; + using pol_state = std::bitset; + static constexpr pol_state pol_state_empty = + 1<update_box(); + return _box; +} +inline Index Polytope::dim() const { return _dim; } +inline Index Polytope::size() const { return _dim; } +inline bool Polytope::is_empty(bool check) const { + if (!check) return state[EMPTY]; + return this->check_empty(); } +inline bool Polytope::is_flat() const { return state[EMPTY] || + _box.is_degenerated() || (_facets && _facets->nbeqfcts()>0); } +inline bool Polytope::is_box() const { + return (!_facets || _facets->nbfcts()==0); +} +inline Index Polytope::nb_facets() const { + if (state[EMPTY]) return -1; + if (!_facets) return (2*_dim); + return (2*_dim + _facets->nbfcts()); + } + +inline Index Polytope::nb_eq_facets() const { + if (state[EMPTY]) return -1; + Index cnt=0; + for (Index i=0;i<_dim;i++) + if (_box[i].is_degenerated()) cnt++; + return (cnt + _facets->nbeqfcts()); + } + +inline void Polytope::set_empty_private() const { + state = pol_state_empty; + _box.set_empty(); + _DDbuildF2V=nullptr; + _DDbuildV2F=nullptr; + _facets=nullptr; +} + + +inline void Polytope::set_empty() { +} + +inline const CollectFacets::mapType &Polytope::facets() const +{ return this->_facets->get_map(); } + +std::ostream& operator<<(std::ostream& os, const Polytope &P); +Polytope operator+ (const Polytope &p1, const Polytope &p2); + +} diff --git a/src/core/domains/polytope/codac2_Polytope_dd.cpp b/src/core/domains/polytope/codac2_Polytope_dd.cpp new file mode 100644 index 000000000..35c73b095 --- /dev/null +++ b/src/core/domains/polytope/codac2_Polytope_dd.cpp @@ -0,0 +1,984 @@ +/** + * \file codac2_Polytope_dd.cpp implementation of the dd algorithm + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_Facet.h" +#include "codac2_Polytope_dd.h" +#include "codac2_IntvFullPivLU.h" +#include "codac2_PosetMax.h" + +#undef DEBUG_CODAC2_DD + +namespace codac2 { + + +/********************/ +/*** DDbuildF2V *****/ +/********************/ + +DDbuildF2V::DDbuildF2V(Index dim, const IntervalVector &box, + std::shared_ptr facetsptr, bool include_box) : dim(dim), + fdim(dim), bbox(box), facetsptr(facetsptr), + empty(false), flat(false), + nbIn(0) { + int nb_eq_dim= (facetsptr==nullptr ? 0 : facetsptr->nbeqfcts()); + if (include_box) { + for (int i=0;i0) { + /* managing the equalities */ + /* we first check the existence of a solution, and get a finite + number of idependant equalities */ + Matrix MatEQ = Matrix::Zero(nb_eq_dim,dim); + IntervalVector RhsEQ = + IntervalVector::Zero(nb_eq_dim); + int base=0; + if (include_box) { + for (int i=0;inbeqfcts();i++) { + MatEQ.row(base+i) = facetsptr->get_eqFacet(i)->first.row; + RhsEQ[base+i] = facetsptr->get_eqFacet(i)->second.rhs; + } + } + IntvFullPivLU LUdec(MatEQ); +#if 0 + if (LUdec.is_surjective)!=BoolInterval::TRUE) { + /* TODO : check emptiness */ + /* IntervalMatrix eqsolve = LUdec.solve(RhsEQ); + if (eqsolve.is_empty()) { empty=true; return; } */ + /* FIXME : NOT CORRECT, as solve DOES NOT + guarantee the emptiness of solutions. We need + to use e.g. a bounding box... wait for + pull request */ + } else +#endif + { + this->fdim = this->dim - nb_eq_dim; + this->M_EQ = IntervalMatrix::Zero(this->dim,this->fdim+1); + this->M_EQ.col(0) = LUdec.solve(RhsEQ); + this->M_EQ.rightCols(this->fdim) = LUdec.kernel(); + this->flat=true; + } +#ifdef DEBUG_CODAC2_DD + std::cout << "M_EQ : " << M_EQ << std::endl; +#endif + } + /* initial vertice : + (1 0 0 ... ), + rays */ + IntervalVector v = IntervalVector::Zero(this->fdim+1); + v[0]=1.0; + this->addDDvertex(v); + for (int i=0;ifdim+1); + line[i+1]=1.0; + lines.push_back(line); + } + if (include_box) this->add_constraint_box(bbox); +} + +int DDbuildF2V::add_bound_var(Index c, bool mx, double rhs) { + if (empty) return -1; + if (mx && rhs==+oo) return 0; + if (!mx && rhs==-oo) return 0; + IntervalRow facet = IntervalRow::Zero(this->fdim+1); + facet[0] = -rhs; + if (flat) { + facet += this->M_EQ.row(c); + } else { + facet[c+1] = 1.0; + } + if (!mx) facet=-facet; + Index idFacet = -c-1-(mx ? 0 : this->dim); + return this->add_facet(idFacet,facet); +} + +int DDbuildF2V::add_constraint_box(const IntervalVector &box) { + int ret=0; + for (Index i=0;iadd_bound_var(i,true,box[i].ub()); + if (ret1==-1) return -1; + if (ret1==1) ret=1; + } + for (Index i=0;iadd_bound_var(i,false,box[i].lb()); + if (ret1==-1) return -1; + if (ret1==1) ret=1; + } + return ret; +} + +std::forward_list::iterator DDbuildF2V::addDDvertex(const IntervalVector &v) { + nbIn++; + vertices.emplace_front(DDvertex(v,nbIn)); + return vertices.begin(); +} +std::forward_list::iterator DDbuildF2V::addDDvertex(IntervalVector &&v) { + nbIn++; + vertices.emplace_front(DDvertex(v,nbIn)); + return vertices.begin(); +} + +inline std::vector fctSon(Index fct,const std::vector fcts1, + const std::vector fcts2) { + std::vector nfcts; + Index a1=0; Index a2=0; + bool idFacetAdded=(fct<0); + while (a1<(Index)fcts1.size() && a2<(Index)fcts2.size()) { + if (fcts1[a1]fct && !idFacetAdded) { + nfcts.push_back(fct); + idFacetAdded=true; + } + nfcts.push_back(fcts1[a1]); + a1++; a2++; + } + } + if (!idFacetAdded) nfcts.push_back(fct); + + return nfcts; +} + +#if 0 +/* p1 has lambda>rhs, p2 has lambda nfcts; + Index a1=0; Index a2=0; + bool idFacetAdded=false; + while (a1<(Index)p1.fcts.size() && a2<(Index)p2.fcts.size()) { + if (p1.fcts[a1]idFacet && !idFacetAdded) { + nfcts.push_back(idFacet); + idFacetAdded=true; + } + nfcts.push_back(p1.fcts[a1]); + a1++; a2++; + } + } + if (!idFacetAdded) nfcts.push_back(idFacet); + std::set lnks; + lnks.insert(p2.Id); + Index newId=this->addDDvertex(newV,lnks,nfcts); + p2.replaceLnks(p1.Id,newId); + return newId; +} +#endif + +IntervalVector DDbuildF2V::compute_vertex(const IntervalVector &vect) const { + if (!flat) { + if (vect[0]==0.0) return vect.tail(this->fdim); + return vect.tail(this->fdim)/vect[0]; + } + if (vect[0]==0.0) return this->M_EQ*vect; + return this->M_EQ*(vect/vect[0]); +} + + +struct CmpDD { + bool operator()(const std::forward_list::iterator &a, + const std::forward_list::iterator &b) { + return (b->vertex[0].lb()*a->lambda.lb()<=a->vertex[0].lb()*b->lambda.lb()); + } +}; + +bool DDbuildF2V::addDDlink(const std::forward_list::iterator& V1, + const std::forward_list::iterator& V2) { + if (V1->addLnk(V2,true)) return V2->addLnk(V1,false); + return false; +// return (V1->addLnk(V2,true) && V2->addLnk(V1,false)); +} + +int DDbuildF2V::add_facet(Index idFacet, const IntervalRow &facet) { + /* consider lines */ + Index bestline=-1; + double maxmig=0.0; + for (Index i=0;i<(Index)lines.size();i++) { + Interval u = facet.dot(lines[i]); + if (u.contains(0.0)) continue; + double mg=u.mig(); + if (mg>maxmig) { + bestline=i; + maxmig=mg; + } + } + if (bestline!=-1) { + Index codeFacet = reffacets.size(); + reffacets.push_back(idFacet); + if (bestline!=(Index)lines.size()-1) { + lines[bestline].swap(lines.back()); + } + IntervalVector& bstline = lines.back(); + Interval ubest = facet.dot(bstline); + /* traitement des lines */ + for (Index i=0;i<(Index)lines.size()-1;i++) { + Interval u = facet.dot(lines[i]); + lines[i] = ubest*lines[i] - u*bstline; + reduce_vector_without_rounding(lines[i]); + } + /* traitement des sommets */ + for (DDvertex &v : vertices) { + double lower = facet.dot(v.vertex).lb(); + if (ubest.lb()>0.0) { + if (lower>0.0) + v.vertex = ubest.ub()*v.vertex - lower*bstline; + else + v.vertex = ubest.lb()*v.vertex - lower*bstline; + } else { + if (lower>0.0) + v.vertex = -ubest.lb()*v.vertex + lower*bstline; + else + v.vertex = -ubest.ub()*v.vertex + lower*bstline; + } + v.reduce_vector(); + v.addFct(codeFacet); + } + /* ajout du nouveau sommet et des liens */ + std::forward_list::iterator itVertex = vertices.begin(); + std::forward_list::iterator newVertex = + this->addDDvertex(ubest.lb()<0.0 ? bstline : (-bstline)); + for (;itVertex!=vertices.end();++itVertex) { + this->addDDlink(itVertex,newVertex); + } + for (Index a=0;aaddFct(a); + } + /* ce nouveau sommet a une face spéciale, "négation" de la contrainte, + de code 0 */ + newVertex->addFct(reffacets.size()); + reffacets.push_back(0); + /* suppression de la ligne */ + this->lines.pop_back(); + return 1; + } + /* end of line treatments , the following + is when no line is available */ + int nbGT=0; + int nbGE=0; + int nbTot=0; + bool notempty=false; + for (DDvertex &vt : vertices) { + vt.lambda = facet.dot(vt.vertex); + notempty |= (vt.lambda.lb()<=0.0); + nbTot++; + if (vt.lambda.ub()<0.0) { + vt.status = DDvertex::VERTEXLT; + continue; + } + if (vt.lambda.lb()<0.0) { + vt.status = DDvertex::VERTEXON; + continue; + } + nbGE++; + if (vt.lambda.lb()tolerance) { + vt.status = DDvertex::VERTEXGE; + continue; + } + nbGT++; + vt.status=DDvertex::VERTEXGT; + } + if (!notempty) { this->empty=true; vertices.clear(); return -1; } + if (nbGT==0) return 0; + + std::forward_list::iterator itVertex = vertices.begin(); + std::priority_queue::iterator, + std::vector::iterator>, + CmpDD> stack; /* note : absolutely not clear that we need + a priority queue */ + Index refFacet=-1; + while (itVertex!=vertices.end()) { +// bool inBoundary; + std::forward_list::iterator actVertex; + if (stack.size()>0) { + actVertex = stack.top(); + stack.pop(); + if (actVertex->status!=DDvertex::VERTEXSTACK) continue; +// inBoundary=true; + } else { + actVertex=itVertex; + ++itVertex; + if (actVertex->status!=DDvertex::VERTEXGT) continue; + refFacet = reffacets.size(); /* new facet identification */ + reffacets.push_back(idFacet); +// inBoundary=false; + } +#ifdef DEBUG_CODAC2_DD + std::cout << " elimination sommet : " << actVertex->Id << ", " << actVertex->links.size() << "liens \n"; +#endif + actVertex->status = DDvertex::VERTEXREM; + std::vector::iterator, + std::vector>> + adjacentVertices; + for (std::forward_list::iterator lnk : actVertex->links) { + DDvertex &destlink = *lnk; +#ifdef DEBUG_CODAC2_DD + std::cout << " lien : " << destlink.Id << "\n"; +#endif + if (destlink.status==DDvertex::VERTEXGT) { + destlink.status=DDvertex::VERTEXSTACK; + stack.push(lnk); + } + if (destlink.status==DDvertex::VERTEXGE || + destlink.status==DDvertex::VERTEXSTACK) { +#ifdef DEBUG_CODAC2_DD + std::cout << " GE ou STACK " << destlink.lambda << "\n"; +#endif + /* add the facet to the vertex, if it does not already exists */ + adjacentVertices.push_back(std::pair(lnk,fctSon(-1,lnk->fcts,actVertex->fcts))); + destlink.addFct(refFacet); + destlink.removeLnk(actVertex); + continue; + } + assert(destlink.status!=DDvertex::VERTEXREM); + /* should not happen */ + IntervalVector partV=destlink.vertex; + if (destlink.status==DDvertex::VERTEXON) { + IntervalRow x1(facet); + MulOp::bwd(Interval(-oo,0.0),x1,partV); + if (partV.is_empty()) { +#ifdef DEBUG_CODAC2_DD + std::cout << " Sommet conserve " << facet << " " << destlink.vertex << " " << destlink.lambda << "\n"; +#endif + adjacentVertices.push_back(std::pair(lnk,fctSon(-1,lnk->fcts,actVertex->fcts))); + destlink.addFct(refFacet); + destlink.removeLnk(actVertex); + continue; + } + } + if (-destlink.lambda.lb()/actVertex->lambda.lb()lambda.lb() << "\n"; +#endif + /* extend the vertex, and add fct */ + adjacentVertices.push_back(std::pair(lnk,fctSon(-1,lnk->fcts,actVertex->fcts))); + destlink.vertex |= + (partV - + (destlink.lambda.lb()/actVertex->lambda)*actVertex->vertex); + destlink.addFct(refFacet); + destlink.removeLnk(actVertex); + continue; + } + /* creation of a new vertex */ + IntervalVector newV = -destlink.lambda.lb()*actVertex->vertex + + actVertex->lambda.lb()*partV; + assert(!newV.is_empty()); + std::forward_list::iterator newVx = + this->addDDvertex(newV); +#ifdef DEBUG_CODAC2_DD + std::cout << " creation new V " << newV << "id : " << newVx->Id << " coef " << (-destlink.lambda.lb()/actVertex->lambda.lb()) << "\n"; +#endif + newVx->lambda=0.0; + destlink.changeLnk(actVertex,newVx); + newVx->addLnk(lnk,false); + newVx->fcts = fctSon(-1,destlink.fcts,actVertex->fcts); + adjacentVertices.push_back(std::pair(newVx,newVx->fcts)); + newVx->addFct(refFacet); + newVx->status=DDvertex::VERTEXGE; + } +#ifdef DEBUG_CODAC2_DD + std::cout << " fin traitement liens : " << adjacentVertices.size() << "\n"; +#endif + PosetMaxSet::iterator, + std::forward_list::iterator>> maxset; + for (Index i=0;i<=(Index)adjacentVertices.size()-1;i++) { +#ifdef DEBUG_CODAC2_DD + std::cout << " sommet : " << adjacentVertices[i].first->Id << " : "; + for (Index u : adjacentVertices[i].second) { + std::cout << u << ","; + } + std::cout << std::endl; +#endif + for (Index j=i+1;j<(Index)adjacentVertices.size();j++) { + PosetElem::iterator, + std::forward_list::iterator>> + actuel(std::pair(adjacentVertices[i].first, + adjacentVertices[j].first), + adjacentVertices[i].second, + adjacentVertices[j].second,refFacet); + if (actuel.elems.size()Id << "\n"; +#endif + for (const PosetElem::iterator, + std::forward_list::iterator>> + &pm : maxset) { +#ifdef DEBUG_CODAC2_DD + if (pm.elems.size()Id << " " << pm.a.second->Id << " " << pm.elems.size() << std::endl; + } + std::cout << pm.a.first->Id << " " << pm.a.second->Id << "\n"; +#endif + [[maybe_unused]] bool res = this->addDDlink(pm.a.first,pm.a.second); +#ifdef DEBUG_CODAC2_DD + std::cout << (res ? " added " : " already here ") << "\n"; +#endif + + } + } + + vertices.remove_if([](const DDvertex &v) + { return v.status==DDvertex::VERTEXREM; }); + return 1; +} + + +int DDbuildF2V::add_facet(Index id) { + return this->add_facet((*facetsptr)[id]); +} + + +int DDbuildF2V::add_facet(CollectFacets::mapCIterator itFacet) { + Index idFacet = itFacet->second.Id; + if (empty) return 0; + if (itFacet->second.eqcst) + throw std::domain_error("add_facet with an equality facet."); + /* the first step is to take the equalities constraints into account */ + IntervalRow facet = IntervalRow::Zero(this->fdim+1); + facet[0] = -itFacet->second.rhs; + if (flat) { + facet += itFacet->first.row * this->M_EQ; + } else { + facet.tail(this->fdim) = itFacet->first.row; + } + return this->add_facet(idFacet,facet); +} + +/* detection of "irredundant facets" */ +std::set DDbuildF2V::redundantFacets() { + std::map> mapping; + std::set resultat; /* set of existing facets */ + for (Index i = -2*dim;i<0;i++) { + resultat.insert(i); + } + for (const auto &u : facetsptr->get_map()) { + resultat.insert(u.second.Id); + } + for (const DDvertex &vt : vertices) { + /* vt.Id is in decreasing order , therefore we consider + the negation of vt.Id to be in increasing order */ + Index id = vt.Id; + for (Index ft : vt.fcts) { + std::vector &elem = mapping[ft]; + elem.push_back(-id); /* in decreasing order */ + } + } + PosetMaxSet facetElems; + for (std::pair> &p : mapping) { +// std::cout << "insert " << p.first << ":" << reffacets[p.first] << " " << p.second.size() << " : "; +// for (Index q : p.second) std::cout << q << "\n"; + facetElems.addElement(p.first,p.second); /* swap p.second */ + } + /* to get the result of facetElems, we store booleans for each facets */ + std::vector redundant(reffacets.size(),true); + for (const PosetElem &elm : facetElems) { + redundant[elm.a]=false; + resultat.erase(reffacets[elm.a]); + /* note : this approach (fill resultat then remove + irredundant facets) is ok when one facet may "appear" + several times */ + } + /* renumbering reffacets */ + std::vector new_reffacets; + Index nb_irredundant=0; + for (unsigned int i=0;i &ret) { + for (Index &a : reffacets) { + if (a<=0) continue; + if (ret[a-1]==-1) { + a=0; + } else { + a=ret[a-1]; + } + } +} + +/* build the bbox of all vertices */ +IntervalVector DDbuildF2V::build_bbox() const { + IntervalVector b = IntervalVector::empty(dim); + for (const DDvertex &a : vertices) { + b |= this->compute_vertex(a.vertex); + } + return b; +} + + +/* normalisation for printing of vertice... good idea ? */ +IntervalVector normalise(const IntervalVector &v) { + IntervalVector res(v); + if (res[0].lb()<=0.0) return res; + return res/res[0].lb(); +} + +std::ostream& operator<<(std::ostream& os, const DDbuildF2V& build) { + os << "DDbuildF2V (" << /*build.vertices.size() <<*/ " vtx) : " << std::endl; + os << " " << build.lines.size() << " lines : " << std::endl; + for (const auto &l : build.lines) { + os << " " << l << std::endl; + } + Index nbvx=0, nbrays=0; + for (const auto &vt : build.vertices) { + if (vt.vertex[0]!=0) nbvx++; else nbrays++; + os << " " << vt.Id << " : " << normalise(vt.vertex) << std::endl; + os << " fct : ("; + for (const auto &a : vt.fcts) { + os << a << ":" << build.reffacets[a] << ","; + } + os << ")" << std::endl; + os << " lnks: ("; + for (const auto &a : vt.links) { + os << a->Id << ","; + } + os << ")" << std::endl; + } + os << "endDDbuildF2V ( "<< nbrays << " rays, " << nbvx << " vertices, " << (nbrays+nbvx) << " total)" << std::endl; + return os; +} + +/*****************************/ +/*** DDbuildF2V::EqFacet *****/ +/*****************************/ +/* management of the equality constraints for DDbuildF2V */ +#if 0 +void DDbuildF2V::EqFacet::adapt_eqconstraint(const IntervalVector &bbox, + Row &cst, Interval &rhs) const { + Interval alpha(cst[dim]); cst[dim]=0.0; + IntervalRow newcst = cst + alpha*this->eqcst; + cst = newcst.mid(); + rhs += (newcst-cst).dot(bbox)-alpha*valcst; +} +void DDbuildF2V::EqFacet::adapt_ineqconstraint(IntervalRow &cst, double &rhs) + const { + Interval alpha(cst[dim]); cst[dim]=0.0; + cst = cst + alpha*this->eqcst; + rhs = (rhs-alpha*valcst).ub(); +} +DDbuildF2V::EqFacet + DDbuildF2V::EqFacet::build_EqFacet(const IntervalVector &bbox, + const Row &cst, Interval &rhs) { + /* identification of the dimension */ + Index dim = -1; + double val=0; + for (Index i=0;ival) { + dim=i; + val = v2; + } + } + if (dim==-1) throw + std::domain_error("zero vector in build_EqFacet"); + IntervalRow newcst(cst); + val=cst[dim]; + newcst[dim]=0.0; + newcst = newcst/val; + rhs = rhs/val - (newcst-newcst.mid()).dot(bbox); + return EqFacet(dim,-newcst.mid(),rhs); +} + +void DDbuildF2V::EqFacet::adapt_box(IntervalVector &bbox) const { + Interval eval = this->eqcst.dot(bbox) + valcst; + if ((bbox[this->dim] & eval).is_empty()) { + bbox.set_empty(); + } + bbox[this->dim]=0; +} +#endif + +/*****************************/ +/*** DDbuildV2F::EqFacet *****/ +/*****************************/ + +DDbuildV2F::DDbuildV2F(Index idVertex, const Vector &vertex) : + facetsptr(std::make_shared()) { + Row r = Row::Zero(vertex.size()); + for (Index i=0;iinsert_facet(r,vertex[i],true); + r[i]=0.0; + } + idVertices.push_back(idVertex); +} + +std::forward_list::iterator DDbuildV2F::addFacetSon( + const std::forward_list::iterator &It1, + std::forward_list::iterator &It2, Index idVertex) { + const CollectFacets::mapIterator &p1 = It1->facetIt; + const CollectFacets::mapIterator &p2 = It2->facetIt; +#ifdef DEBUG_CODAC2_DD + std::cout << " p1.facet.row " << p1->first.row << std::endl; + std::cout << " p2.facet.row " << p2->first.row << std::endl; +#endif + Row newV = - It2->lambda*p1->first.row + It1->lambda*p2->first.row; + double newRhs = - It2->lambda*p1->second.rhs + It1->lambda*p2->second.rhs; +#ifdef DEBUG_CODAC2_DD + std::cout << " newV " << newV << std::endl; + std::cout << " newRhs " << newRhs << std::endl; +#endif + /* computation of vertices */ + std::vector nvtx; + Index a1=0; Index a2=0; + bool idVertexAdded=false; + while (a1<(Index)It1->vtx.size() && a2<(Index)It2->vtx.size()) { + if (It1->vtx[a1]vtx[a2]) { + a1++; if (a1==(Index)It1->vtx.size()) break; + } else if (It2->vtx[a2]vtx[a1]) { + a2++; if (a2==(Index)It2->vtx.size()) break; + } else { + if (It1->vtx[a1]>idVertex && !idVertexAdded) { + nvtx.push_back(idVertex); + idVertexAdded=true; + } + nvtx.push_back(It1->vtx[a1]); + a1++; a2++; + } + } + if (!idVertexAdded) nvtx.push_back(idVertex); + std::vector::iterator> lnks; + lnks.push_back(It2); + std::forward_list::iterator + newF=this->addDDfacet(reduce_facet(newV,newRhs),lnks,nvtx); + if (newF!=facets.end()) It2->changeLnk(It1,newF); + else It2->removeLnk(It1); + return newF; +} + +int DDbuildV2F::add_vertex(Index idVertex, const Vector& vertex) { + return this->add_vertex(idVertex,IntervalVector(vertex)); +} + +int DDbuildV2F::add_vertex(Index idVertex, const IntervalVector& vertex) { + int cnt=0; + /* first check equalities */ + if (facetsptr->nbeqfcts()>0) { + Index eqViolated=-1; + CollectFacets::mapIterator eqViolatedIt= facetsptr->endFacet(); + Interval delta; + for (Index i=facetsptr->nbeqfcts()-1;i>=0;i--) { + CollectFacets::mapIterator fcIt = facetsptr->get_eqFacet(i); + Interval calc = fcIt->first.row.dot(vertex)-fcIt->second.rhs; + if (calc.mig()tolerance) continue; + if (eqViolated==-1) { + eqViolated=i; + delta=calc; + eqViolatedIt = fcIt; + continue; + } else { + double correction = -(calc/delta).mid(); + bool success; + std::tie(fcIt,success) = facetsptr->change_eqFacet(i, + fcIt->first.row + correction*eqViolatedIt->first.row, + fcIt->second.rhs + correction*eqViolatedIt->second.rhs); + /* if success is false, the eqFacet index may have changed */ + if (!success && + eqViolated==facetsptr->nbeqfcts()) eqViolated=i; + } + } + if (eqViolated!=-1) { + if (facets.empty()) { + this->nbIn=0; + std::vector::iterator> lnk; + std::vector vtx = { idVertex }; + Interval nrhsI = eqViolatedIt->first.row.dot(vertex); + double nrhs = delta.lb()>0.0 ? nrhsI.ub() : nrhsI.lb(); + CollectFacets::mapIterator newIt = + facetsptr->dissociate_eqFacet(eqViolated,nrhs); + assert_release(newIt!=facetsptr->endFacet()); + if (delta.lb()>0.0) { + auto a = this->addDDfacet(eqViolatedIt,lnk,vtx); + vtx = idVertices; + lnk = { a }; + auto b = this->addDDfacet(newIt,lnk,vtx); + a->addLnk(b,false); + } else { + auto a = this->addDDfacet(newIt,lnk,vtx); + vtx = idVertices; + lnk = { a }; + auto b = this->addDDfacet(eqViolatedIt,lnk,vtx); + a->addLnk(b,false); + } + cnt=2; + } else { + /* add the remaining inequality, linked with all existing facets */ + std::vector::iterator> lnk; + for (std::forward_list::iterator it = + facets.begin(); it!=facets.end(); ++it) { + it->status=DDfacet::FACETOUT; /* to remove FACETNEW */ + lnk.push_back(it); + } + Row saveEqrow = eqViolatedIt->first.row; + double saveRhs = eqViolatedIt->second.rhs; + CollectFacets::mapIterator nIt = + facetsptr->dissociate_eqFacet(eqViolated, + (delta.lb()>0.0 ? +oo : -oo)); +// std::cout << "dissociated : " << (delta.lb()>0.0 ? +oo : -oo) << +// " : " << (nIt == facetsptr->endFacet()) << "\n"; + eqViolatedIt=nIt; /* in the following, the case nIt==facets.end() + is not treated, as it should not happen */ + std::vector vtx = idVertices; + + std::forward_list::iterator u = + (eqViolatedIt != facetsptr->endFacet() ? + this->addDDfacet(eqViolatedIt,lnk,vtx) : + facets.end()); + /* change all existing facets to take the new vertex into account */ + /* no facet should be removed, theoretically, but we + * nevertheless consider this possibility */ + std::forward_list::iterator it_f=facets.before_begin(); + std::forward_list::iterator nxt = facets.begin(); + while (nxt!=facets.end()) { + DDfacet &f = (*nxt); + if (f.status!=DDfacet::FACETOUT) { it_f=nxt; ++nxt; continue; } + Interval q = f.facetIt->second.rhs- + f.facetIt->first.row.dot(vertex); + double correction = (q/delta).mid(); + if (u!=facets.end()) f.addLnk(u); + f.addVtx(idVertex); + bool success; + std::tie(f.facetIt, success) = + facetsptr->change_ineqFacet(f.facetIt->second.Id, + f.facetIt->first.row + correction*saveEqrow, + f.facetIt->second.rhs + correction*saveRhs); + if (!success) { this->removeFacet(nxt,true); + facets.erase_after(it_f); + nxt=it_f; ++nxt; continue; } + it_f=nxt; ++nxt; + } + cnt=(u!=facets.end() ? 1 : 0); + } + if (facetsptr->nbeqfcts()>0) idVertices.push_back(idVertex); + return cnt; + } else idVertices.push_back(idVertex); /* and follows */ + } + /* no equalities removed */ + for (DDfacet &d : facets) { + if (d.status==DDfacet::FACETREM) continue; + Interval calc = d.facetIt->first.row.dot(vertex)-d.facetIt->second.rhs; + d.lambda = calc.lb(); + if (calc.mig()tolerance) { d.status=DDfacet::FACETON; } + else if (calc.lb()>0.0) { d.status=DDfacet::FACETOUT; } + else { d.status=DDfacet::FACETIN; } + } + for (std::forward_list::iterator it = facets.begin(); + it!=facets.end();++it) { + DDfacet &d = (*it); + if (d.status==DDfacet::FACETON) { + d.addVtx(idVertex); + continue; + } + if (d.status==DDfacet::FACETIN || d.status==DDfacet::FACETNEW) continue; + if (d.links.empty()) { /* not supposed to happen */ + /* very specific case of a 1-dimension polyhedra */ + d.vtx.clear(); + d.vtx.push_back(idVertex); + facetsptr->change_ineqFacet_rhs(d.facetIt, + d.facetIt->first.row.dot(vertex).ub()); + continue; + } + for (std::forward_list::iterator idL : d.links) { + if (idL==facets.end()) continue; /* can it happen ? */ + DDfacet &d2 = (*idL); + if (d2.status==DDfacet::FACETON) { + d2.removeLnk(it); + continue; + } + if (d2.status!=DDfacet::FACETIN) { + continue; /* do nothing for out */ + } +#ifdef DEBUG_CODAC2_DD + std::cout << " addFacet " << d.facetIt->second.Id << " " << d2.facetIt->second.Id << " " << d.status << " " << d2.status << std::endl; + std::cout << " (" << d.lambda << " " << d2.lambda << ") " << std::endl; +#endif + [[maybe_unused]] std::forward_list::iterator + newIt = this->addFacetSon(it,idL,idVertex); +#ifdef DEBUG_CODAC2_DD + if (newIt==facets.end()) + std::cout << " ... FAILED! " << std::endl; + else + std::cout << " ... done " << newIt->facetIt->second.Id << std::endl; +#endif + cnt++; + } + facetsptr->removeFacetById(it->facetIt->second.Id); + this->removeFacet(it,false); + } + /*construction of new links for the eq facets */ + PosetMaxSet::iterator, + std::forward_list::iterator>> maxset; + for (std::forward_list::iterator it = facets.begin(); + it!=facets.end();++it) { + DDfacet &d = (*it); + if (d.status==DDfacet::FACETON || d.status==DDfacet::FACETNEW) { + for (std::forward_list::iterator it2 = next(it); + it2!=facets.end(); ++it2) { + DDfacet &d2 = (*it2); + if (d2.status==DDfacet::FACETON || d2.status==DDfacet::FACETNEW) { + PosetElem::iterator, + std::forward_list::iterator>> + actuel(std::pair(it,it2),d.vtx,d2.vtx); + maxset.addElement(std::move(actuel)); + } + } + } + } + for (const PosetElem::iterator, + std::forward_list::iterator>> &pm : maxset) { +// std::cout << " maxset : " << pm.a << " " << pm.b << std::endl; + pm.a.first->addLnk(pm.a.second); + pm.a.second->addLnk(pm.a.first); + } + facets.remove_if([](const DDfacet &v) + { return v.status==DDfacet::FACETREM; }); + facetsptr->renumber(); + return cnt; +} + +std::forward_list::iterator + DDbuildV2F::addDDfacet(const CollectFacets::mapIterator &it, + std::vector::iterator> &links, + std::vector &vtx) { + facets.emplace_front(it,links,vtx); + return facets.begin(); +} + +std::forward_list::iterator + DDbuildV2F::addDDfacet(const std::pair &row_rhs, + std::vector::iterator> &links, + std::vector &vtx) { + std::pair result = + facetsptr->insert_facet(row_rhs.first,row_rhs.second,false); + if (!result.second) { + return facets.end(); + } + facets.emplace_front((*facetsptr)[result.first-1],links,vtx); + return facets.begin(); +} + +std::forward_list::iterator + DDbuildV2F::addDDfacet(std::pair &&row_rhs, + std::vector::iterator> &links, + std::vector &vtx) { + std::pair result = + facetsptr->insert_facet(row_rhs.first,row_rhs.second,false); + if (!result.second) { + return facets.end(); + } + facets.emplace_front((*facetsptr)[result.first-1],links,vtx); + return facets.begin(); +} + +void DDbuildV2F::removeFacet(std::forward_list::iterator It, + bool removeLnks) { + DDfacet &d = (*It); + if (removeLnks) { + for (std::forward_list::iterator it : d.links) { + it->removeLnk(It); + } + } + d.facetIt=facetsptr->endFacet(); + d.status=DDfacet::FACETREM; +} + +std::ostream& operator<<(std::ostream& os, const DDbuildV2F& build) { + os << "DDbuildV2F (" << build.facetsptr->nbeqfcts() << " eq, " << + build.facetsptr->nbfcts() << " total fcts) : " << std::endl; + for (Index i=0;inbeqfcts();++i) { + const Facet &eqf = (*build.facetsptr->get_eqFacet(i)); + os << eqf.first.get_row() << " = " << eqf.second.get_rhs() << std::endl; + } + for (const auto &ifacet : build.facets) { + auto &ieqf = *(ifacet.facetIt); + os << ieqf.second.get_Id() << " : " << ieqf.first.get_row() << " <= " << ieqf.second.get_rhs() << std::endl; + os << " vtx: ("; + for (Index a : ifacet.vtx) { + os << a << ","; + } + os << ")" << std::endl; + os << " lnks: ("; + for (std::forward_list::iterator a : ifacet.links) { + os << a->facetIt->second.get_Id() << ","; + } + os << ")" << std::endl; + } + os << "endDDbuildV2F" << std::endl; + return os; +} + +} + diff --git a/src/core/domains/polytope/codac2_Polytope_dd.h b/src/core/domains/polytope/codac2_Polytope_dd.h new file mode 100644 index 000000000..29636dbf7 --- /dev/null +++ b/src/core/domains/polytope/codac2_Polytope_dd.h @@ -0,0 +1,373 @@ +/** + * \file codac2_Polytope_dd.h implementation of the dd algorithms + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_Facet.h" + + +namespace codac2 { + + +/** \brief multiply a vector by a power of 2 (without rounding) + * such that 2^(-10) <= norm(vec) <= 2^10 + * \param vec the vector (modified) + * \return the multiplying factor */ +inline double reduce_vector_without_rounding(IntervalVector &vec) { + int fx; + if (frexp(vec[0].mag(),&fx)==0.0) fx=INT_MIN; + for (Index i=1;ifx) fx=fx2; + } + assert_release (fx!=INT_MIN); + if (std::abs(fx)<10) return 1.0; + double mult=exp2(-fx); + for (Index i=0;ifx) fx=fx2; + } + assert_release (fx!=INT_MIN); + if (std::abs(fx)<10) return 1.0; + double mult=exp2(-fx); + for (Index i=0;i fcts; + /* note : double redirection here, a ``real'' facet may have several + Ids */ + + enum ddvertexstatus { + VERTEXNEW, /* created */ + VERTEXLT, /* < (GT,LT) => new vertex */ + VERTEXON, /* <,=,> (GT,LE) => for now, equivalent to LT */ + VERTEXGE, /* = and possible small > + keep the vertex. no change with (GT,GE) */ + VERTEXGT, /* vertex GT, untreated and not in stack */ + VERTEXSTACK, /* vertex in stack for later treatment */ + VERTEXREM /* vertex GT treated, to be removed */ + }; + ddvertexstatus status; + + std::vector::iterator> links; + + DDvertex(const IntervalVector &v, Index Id) : Id(Id), + vertex(v), status(VERTEXNEW) { + reduce_vector(); + } + + DDvertex(IntervalVector &&v, Index Id, + std::vector::iterator> &nlinks) : + Id(Id), + vertex(v), status(VERTEXNEW) { + this->links.swap(nlinks); + reduce_vector(); + } + + inline void reduce_vector() { + reduce_vector_without_rounding(vertex); + } + + bool addLnk(const std::forward_list::iterator &a, + bool check=false) { + if (check) { + for (auto b : links) { + if (a==b) return false; + } + } + this->links.push_back(a); + return true; + } + +// void addLnkVertex(const std::forward_list::iterator &V1); + + void removeLnk(const std::forward_list::iterator &a) { + for (Index i=0;i<(Index)links.size();i++) { + if (links[i]==a) { + links[i]=links[links.size()-1]; + links.pop_back(); + return; + } + } + } + + void changeLnk(const std::forward_list::iterator &a, + const std::forward_list::iterator &b) { + for (Index i=0;i<(Index)links.size();i++) { + if (links[i]==a) { links[i]=b; return; } + } + /* should not happen */ + links.push_back(b); + } + + bool addFct(Index a) { + auto rev_it = this->fcts.rbegin(); + while (rev_it != this->fcts.rend() && (*rev_it)>a) rev_it++; + if (rev_it!=this->fcts.rend() && (*rev_it)==a) return false; + this->fcts.insert(rev_it.base(),a); + return true; + } + +}; + +/** \brief Structure to build an enclosing set of IntervalVector + * from a set of facets */ +class DDbuildF2V { + public: + /* create the build with a collection of facet, but + use only the equality facets */ + DDbuildF2V(Index dim, const IntervalVector &box, + std::shared_ptr facetsptr, + bool include_box=true); + + /* add an inequality facet in the build */ + /* return 1 : changed, 0 : not changed, -1 : empty */ + int add_facet(Index id); + int add_facet(CollectFacets::mapCIterator itFacet); + + int add_constraint_box(const IntervalVector &box); + int add_bound_var(Index c, bool mx, double rhs); + + IntervalVector compute_vertex(const IntervalVector &vect) const; + Index get_dim() const; + Index get_fdim() const; + + bool is_empty() const; + const IntervalMatrix &get_M_EQ() const; + + /* return a set of the Id of (detected as) redundant facets , + and remove these facets of the representation */ + std::set redundantFacets(); + + /* update reffacets following a renumbering */ + void update_renumber(const std::vector &ref); + + /* return the bounding box of all vertices */ + IntervalVector build_bbox() const; + + const std::vector &get_lines() const; + const std::forward_list &get_vertices() const; + const std::vector &get_reffacets() const; + + friend std::ostream& operator<<(std::ostream& os, const DDbuildF2V& build); + + private: + int add_facet(Index idFacet, const IntervalRow &facet); + bool addDDlink(const std::forward_list::iterator& V1, + const std::forward_list::iterator& V2); + + std::forward_list::iterator addDDvertex(const IntervalVector &v); + std::forward_list::iterator addDDvertex(IntervalVector &&v); +/* std::map::iterator addVertexSon(const DDvertex &p1, + DDvertex &p2, double rhs, Index idFacet); */ +/* void removeVertex(Index Id, bool removeLnks); */ + + Index dim; + Index fdim; /* dimension - nb independant equalities constraints */ + IntervalVector bbox; + + std::shared_ptr facetsptr; + + /* if fdim lines; /* lines (for unbounded polyhedron) */ + /* rays can be handled as classical vertices */ + + std::forward_list vertices; +#if 0 + std::forward_list total_links; +#endif + std::vector reffacets; /* indirection for facets */ + + bool empty, flat; + Index nbIn; + double tolerance = 1e-9; +}; + +inline Index DDbuildF2V::get_dim() const { return this->dim; } +inline Index DDbuildF2V::get_fdim() const { return this->fdim; } +inline bool DDbuildF2V::is_empty() const { return this->empty; } +inline const IntervalMatrix &DDbuildF2V::get_M_EQ() const { return this->M_EQ; } +inline const std::vector &DDbuildF2V::get_lines() const + { return this->lines; } +inline const std::forward_list &DDbuildF2V::get_vertices() const + { return this->vertices; } +inline const std::vector &DDbuildF2V::get_reffacets() const + { return this->reffacets; } + + +/** \brief structure for facets in BuildF2V */ +struct DDfacet { + CollectFacets::mapIterator facetIt; + double lambda; + + enum ddfacetstatus { + FACETNEW, + FACETIN, + FACETON, + FACETOUT, + FACETREM + }; + ddfacetstatus status; + + std::vector::iterator> links; + std::vector vtx; + + DDfacet(CollectFacets::mapIterator facetIt, + std::vector::iterator> &nlinks, + std::vector &nvtx) : facetIt(facetIt), status(FACETNEW) { + this->links.swap(nlinks); + this->vtx.swap(nvtx); + } + + bool addLnk(std::forward_list::iterator a, bool check=false) { + if (check) { + for (auto b : links) { + if (a==b) return false; + } + } + this->links.push_back(a); + return true; + } + + void removeLnk(const std::forward_list::iterator &a) { + for (Index i=0;i<(Index)links.size();i++) { + if (links[i]==a) { + links[i]=links[links.size()-1]; + links.pop_back(); + return; + } + } + } + + void changeLnk(const std::forward_list::iterator &a, + const std::forward_list::iterator &b) { + for (Index i=0;i<(Index)links.size();i++) { + if (links[i]==a) { links[i]=b; return; } + } + /* should not happen */ + links.push_back(b); + } + + void addVtx(Index a) { + auto it = this->vtx.begin(); + while (it != this->vtx.end() && (*it)vtx.end() && (*it)==a) return; + this->vtx.insert(it,a); + } + +}; + +/** \brief construction of an (approximate) set of facets from a + * a set of vertices */ +class DDbuildV2F { + public: + DDbuildV2F(Index idVertex, const Vector &vertex); + + int add_vertex(Index idVertex, const Vector& vertex); + int add_vertex(Index idVertex, const IntervalVector& vertex); + + std::shared_ptr getFacets(); + + friend std::ostream& operator<<(std::ostream& os, const DDbuildV2F& build); + + + private: + /* can NOT return facets.end() */ + std::forward_list::iterator addDDfacet( + const CollectFacets::mapIterator &it, + std::vector::iterator> &links, + std::vector &vtx); + /* can return facets.end() if insertion failed */ + std::forward_list::iterator addDDfacet( + std::pair &&row_rhs, + std::vector::iterator> &links, + std::vector &vtx); + /* can return facets.end() if insertion failed */ + std::forward_list::iterator addDDfacet( + const std::pair &row_rhs, + std::vector::iterator> &links, + std::vector &vtx); + /* can return facets.end() if insertion failed */ + std::forward_list::iterator + addFacetSon(const std::forward_list::iterator &It1, + std::forward_list::iterator &It2, Index idVertex); + void removeFacet(std::forward_list::iterator It, + bool removeLnks); + + inline std::pair reduce_facet(const Row &row, double rhs) { + std::pair ret(row,rhs); + double mult = reduce_row_without_rounding(ret.first); + ret.second *= mult; + return ret; + } + + std::shared_ptr facetsptr; + std::forward_list facets; + std::vector idVertices; /* list of vertices entered, useful as + long as there are eqfacets */ + Index nbIn=0; + double tolerance =1e-8; + +}; + +inline std::shared_ptr DDbuildV2F::getFacets() { return facetsptr; } + +} + diff --git a/src/core/domains/polytope/codac2_Polytope_util.cpp b/src/core/domains/polytope/codac2_Polytope_util.cpp new file mode 100644 index 000000000..c7baf4ce6 --- /dev/null +++ b/src/core/domains/polytope/codac2_Polytope_util.cpp @@ -0,0 +1,194 @@ +/** + * \file codac2_Polytope_util.cpp utilities for polytope + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_Facet.h" +#include "codac2_Polytope_dd.h" +#include "codac2_Polytope_util.h" + + +namespace codac2 { + +/** read a .ine file and return a list of facets + * this is a very crude function, which assumes the format of the + * file is correct, and no real check is done. Linearities are not treated. */ +CollectFacets read_ineFile(const char *filename) { + std::ifstream fichier(filename); + assert(!fichier.fail()); + /* skip everything until begin */ + std::string line; + while (getline(fichier,line)) { + std::cout << line << std::endl; + if (line.compare("begin")==0) break; + } + assert(fichier.good()); + Index nbfacets, dim; + char dummy[50]; + fichier >> nbfacets >> dim >> dummy; + std::cout << nbfacets << " facettes " << dim << " dimensions" << std::endl; + dim = dim-1; /* non-normalised dimension */ + CollectFacets result; + for (Index i=0;i> rhs; + for (Index j=0;j> row[j]; + } + row = -row; + result.insert_facet(row,rhs,false); + } + fichier.close(); + return result; +} + +/** read a .ext file and return a list of vertices + * this is a very crude function, which assumes the format of the + * file is correct, and no real check is done. Rays (or lines) + * are not treated for now. */ +std::vector read_extFile(const char *filename) { + std::ifstream fichier(filename); + assert(!fichier.fail()); + /* skip everything until begin */ + std::string line; + while (getline(fichier,line)) { + std::cout << line << std::endl; + if (line.compare("begin")==0) break; + } + assert(fichier.good()); + Index nbvertices, dim; + char dummy[50]; + fichier >> nbvertices >> dim >> dummy; + std::cout << nbvertices << " sommets " << (dim-1) << " dimensions" << std::endl; + dim = dim-1; /* non-normalised dimension */ + std::vector result; + for (Index i=0;i> rhs; /* 'rhs' is supposed to be 1 */ + for (Index j=0;j> vrt[j]; + } + result.push_back(vrt); + } + fichier.close(); + return result; +} + +/* create a 2D facet in a 3D polyhedron, return the number of vertices built. + fct = -1 : any facet */ +Index fill_3Dfacet(const DDbuildF2V &build, + std::vector &tofill, Index fct, double bound) { + const std::forward_list &vertices=build.get_vertices(); + std::forward_list::const_iterator it=vertices.begin(); + while (it != vertices.end()) { + if (it->vertex[0]!=0.0 && + (fct==-1 || std::find(it->fcts.begin(), + it->fcts.end(),fct)!=it->fcts.end())) + break; + ++it; + } + if (it==vertices.end()) return 0; + + std::vector::const_iterator> used; + used.push_back(it); + Vector vt = build.compute_vertex(it->vertex).mid(); + tofill.push_back(vt); + while (it != vertices.end()) { + std::forward_list::const_iterator it2; + Index i=0; + while (i<(Index)it->links.size()) { + it2 = it->links[i]; + if ((fct==-1 || + std::find(it2->fcts.begin(), + it2->fcts.end(),fct)!=it2->fcts.end()) + && (std::find(used.begin(),used.end(),it2)==used.end())) break; + i++; + } + if (i==(Index)it->links.size()) break; + it=it2; + if (it->vertex[0]==0.0) /* ray */ { + vt += bound*build.compute_vertex(it->vertex).mid(); + } else { + vt = build.compute_vertex(it->vertex).mid(); + } + tofill.push_back(vt); + used.push_back(it); +// std::cout << "add vert : " << it->Id << "\n"; + } + return (Index)tofill.size(); +} + +std::vector> build_3Dfacets(const DDbuildF2V &build, double bound) { + assert(build.get_dim()==3); + if (build.is_empty()) return std::vector>(); + if (build.get_fdim()==0) { /* the only point has no coordinate */ + std::vector> result(1); + Vector v = Vector::constant(1,1.0); + result[0].push_back((build.get_M_EQ()).mid()*v); + return result; + } else + if (build.get_fdim()<3) { /* line or face are managed as a single face */ + std::vector> result(1); + std::vector &face = result[0]; + fill_3Dfacet(build,face,-1,bound); + /* look for a initial vertex */ + return result; + } + std::vector> result; + for (Index i=0;i<(Index)build.get_reffacets().size();i++) { + if (build.get_reffacets()[i]==0) continue; + std::vector nfac; + Index nbV = fill_3Dfacet(build, nfac, i, bound); +// std::cout << "facet build, " << nbV << " vertices.\n"; + if (nbV<3) continue; + result.push_back(nfac); + } + return result; +} + +std::vector build_2Dfacet(const DDbuildF2V &build, double bound) { + assert(build.get_dim()==2); + if (build.is_empty()) return std::vector(); + if (build.get_fdim()==0) { /* the only point has no coordinate */ + std::vector result(1); + Vector v = Vector::constant(1,1.0); + result[0]=(build.get_M_EQ()).mid()*v; + return result; + } else + { /* equivalent to finding one facet of a 3Dfacet */ + std::vector face; + fill_3Dfacet(build,face,-1,bound); + /* look for a initial vertex */ + return face; + } +} + +} + + diff --git a/src/core/domains/polytope/codac2_Polytope_util.h b/src/core/domains/polytope/codac2_Polytope_util.h new file mode 100644 index 000000000..399fb7633 --- /dev/null +++ b/src/core/domains/polytope/codac2_Polytope_util.h @@ -0,0 +1,69 @@ +/** + * \file codac2_Polytope_util.h utilities for polytope + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "codac2_Index.h" +#include "codac2_Vector.h" +#include "codac2_Row.h" +#include "codac2_Facet.h" +#include "codac2_Polytope_dd.h" + + +namespace codac2 { + +/** read a .ine file and return a list of facets + * this is a very crude function, which assumes the format of the + * file is correct. Linearities are not treated for now. + * @param filename name of the file + * @return the CollectFacets structure containing the facets */ +CollectFacets read_ineFile(const char *filename); + +/** read a .ext file and return a list of vertices + * this is a very crude function, which assumes the format of the + * file is correct. Rays are not treated for now. + * @param filename name of the file + * @return a set of vertices */ +std::vector read_extFile(const char *filename); + +/** construct the list of facets of a 3D-polyhedron from the buildF2V + * structure + * use the mid of each vertex. For rays/line use an arbitrary bound + * @param build the buildF2V structure + * @param bound the bound for rays/line (not correctly handled) + * @return the set of set of vertices + */ +std::vector> build_3Dfacets(const DDbuildF2V &build, double bound=50.0); + +/** construct the facet of a 2D-polyhedron from the buildF2V + * structure + * use the mid of each vertex. For rays/line use an arbitrary bound + * @param build the buildF2V structure + * @param bound the bound for rays/line (not correctly handled) + * @return the set of set of vertices + */ +std::vector build_2Dfacet(const DDbuildF2V &build, double bound=50.0); + +} + diff --git a/src/core/domains/polytope/codac2_PosetMax.h b/src/core/domains/polytope/codac2_PosetMax.h new file mode 100644 index 000000000..c2f443777 --- /dev/null +++ b/src/core/domains/polytope/codac2_PosetMax.h @@ -0,0 +1,178 @@ +/** + * \file codac2_PosetMax.h Managing a set of maximal elements in a poset + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#pragma once + +#include +#include +#include +#include + +#include "codac2_Index.h" + +namespace codac2 { + /* the poset is an (ordered) vector of Index */ + /* used in .h because of its templated version */ + +/*** elements ****************************************************/ +/*** we may store more that set of elems to speed up computation */ + +/** \brief Pair of an element of type T and an (ordered) vector of + * Index, along with the partial (subset) order on the vector */ +template +struct PosetElem { + T a; + std::vector elems; + + /** \brief empty element */ + PosetElem() {} + + /** \brief Element with empty list + * \param a key of type T */ + PosetElem(const T &a) : a(a) { + } + + /** \brief Construct an element (a,e1) + * \param a key + * \param e1 set of indices (supposed to be ordered) */ + PosetElem(const T &a, const std::vector &e1) : a(a), elems(e1) { + } + + /** \brief build an element (a, e1 inter e2) + * \param a key + * \param e1 ordered set + * \param e2 ordered set */ + PosetElem(const T &a, const std::vector &e1, const std::vector &e2) : a(a) { + Index i=0,j=0; + while(i<(Index)e1.size() && j<(Index)e2.size()) { + if (e1[i]==e2[j]) { + elems.push_back(e1[i]); + i++; j++; + } else if (e1[i] &e1, const std::vector &e2, const Index rejectF) : a(a) { + Index i=0,j=0; + while(i<(Index)e1.size() && j<(Index)e2.size()) { + if (e1[i]==e2[j]) { + if (e1[i]!=rejectF) + elems.push_back(e1[i]); + i++; j++; + } else if (e1[i]a,p2.a); + this->elems.swap(p2.elems); + } + + /** \brief comparison between the set of indices + * \param pm the element to be compared + * \return -1 if (*this).elems is included in pm.elems + * 0 if both sets are equal + * 1 if (*this).elems is a superset of pm.elems + * 2 if both sets are not comparable + */ + int comp(const PosetElem &pm) const { + /* -1 if elems included in pm.elems; 0 if ==, + +1 elems includes pm.elems, +2 if not comparable */ + const std::vector &e1=elems; + const std::vector &e2=pm.elems; + std::vector::size_type i=0,j=0; + bool leq=(e1.size()<=e2.size()), geq=(e1.size()>=e2.size()); + bool same_length = (e1.size()==e2.size()); + /* leq=F,geq=T => +1 + leq=F,geq=F => 2 + leq=T,geq=F => -1 + leq=T,geq=T => 0 */ + while(ie2[j]) { + if (!leq || same_length) return 2; + geq=false; j++; + } else { i++; j++; } + } + if (i +struct PosetMaxSet : std::vector> { + + /** \brief build and add a new element + * \param a key + * \param elems the set of indices (ordered) + * \return true if the element is added, false otherwise + */ + bool addElement(const T &a, std::vector &elems) { + PosetElem e(a); + std::swap(e.elems,elems); + return this->addElement(std::move(e)); + } + + /** \brief move in a new element (if maximal) + * \param elem the element + * \return true if the element is added, false otherwise + */ + bool addElement(PosetElem &&elem) { + bool ismax=true; + unsigned long int k=0,l=this->size(); + while (ksize()) { + this->push_back(std::move(elem)); + } else { + (*this)[l] = std::move(elem); + this->resize(l+1); + } + } else if (l!=this->size()) { /* should not happen */ + this->resize(l); + } + return ismax; + } +}; + +} diff --git a/src/core/functions/analytic/codac2_AnalyticFunction.h b/src/core/functions/analytic/codac2_AnalyticFunction.h index 2a527fbf1..3c4a020c1 100644 --- a/src/core/functions/analytic/codac2_AnalyticFunction.h +++ b/src/core/functions/analytic/codac2_AnalyticFunction.h @@ -21,6 +21,7 @@ #include "codac2_vec.h" #include "codac2_Wrapper.h" #include "codac2_Parallelepiped.h" +#include "codac2_Polytope.h" #include "codac2_peibos_tools.h" namespace codac2 @@ -221,6 +222,48 @@ namespace codac2 return Parallelepiped(z, A_inf); } + template + requires std::is_same_v + Polytope polytope_eval(const Args&... x) const + { + this->check_valid_inputs(x...); + + /* use centered evaluation */ + auto x_ = eval_(x...); + if(x_.da.size() == 0 || !x_.def_domain) { + return Polytope(eval(EvalMode::NATURAL, x...)); + } + auto flatten_x = IntervalVector(cart_prod(x...)); + assert(x_.da.rows() == x_.a.size() && + x_.da.cols() == flatten_x.size()); + + IntervalVector centered_flatten = flatten_x-flatten_x.mid(); + Polytope res = Polytope(centered_flatten).affine_transform + (x_.da, x_.m, x_.a); + return res; + } + + template + requires std::is_same_v + Polytope polytope_eval(const Polytope& input) const + { + IntervalVector box = input.box(); + this->check_valid_inputs(box); + + /* use centered evaluation */ + auto x_ = eval_(box); + if(x_.da.size() == 0 || !x_.def_domain) + return Polytope(eval(EvalMode::NATURAL, box)); + assert(x_.da.rows() == x_.a.size() && + x_.da.cols() == box.dim()); + + Polytope input_copy(input); + input_copy.inflate(IntervalVector(-box.mid())); + Polytope res = input_copy.direct_affine_transform(x_.da, x_.m); + res.meet_with_box(x_.a); + return res; + } + Index output_size() const { if constexpr(std::is_same_v) @@ -340,4 +383,4 @@ namespace codac2 AnalyticFunction(const FunctionArgsList&, const T&) -> AnalyticFunction::Type>; -} \ No newline at end of file +} diff --git a/src/core/matrices/codac2_IntvFullPivLU.h b/src/core/matrices/codac2_IntvFullPivLU.h index 108cfaa39..c65591e57 100644 --- a/src/core/matrices/codac2_IntvFullPivLU.h +++ b/src/core/matrices/codac2_IntvFullPivLU.h @@ -48,7 +48,7 @@ namespace codac2 /** * \brief Check if the matrix is injective, - * *i.e.* its rank is equal to its number of rows. + * *i.e.* its rank is equal to its number of cols. * * \return a ``BoolInterval`` */ @@ -64,7 +64,7 @@ namespace codac2 /** * \brief Check if the matrix is surjective - * *i.e.* its rank is equal to its number of cols. + * *i.e.* its rank is equal to its number of rows. * * \return a ``BoolInterval`` */ @@ -339,4 +339,4 @@ inline Derived IntvFullPivLU::coimage return ret; } -} \ No newline at end of file +} diff --git a/src/core/matrices/codac2_inversion.h b/src/core/matrices/codac2_inversion.h index 9f014c7b9..cbdd6e655 100644 --- a/src/core/matrices/codac2_inversion.h +++ b/src/core/matrices/codac2_inversion.h @@ -34,14 +34,13 @@ namespace codac2 template inline IntervalMatrix inverse_correction(const Eigen::MatrixBase& A, const Eigen::MatrixBase& B) { - assert_release(A.is_squared()); - assert_release(B.is_squared()); + assert_release(B.cols()==A.rows()); + assert_release(A.cols()==B.rows()); auto A_ = A.template cast(); auto B_ = B.template cast(); - Index N = A_.rows(); - assert_release(N==B_.rows()); + Index N = (O == LEFT_INV ? A_.cols() : A_.rows()); auto Id = IntervalMatrix::Identity(N,N); auto erMat = [&]() { if constexpr(O == LEFT_INV) return -B_*A_+Id; else return -A_*B_+Id; }(); @@ -60,12 +59,10 @@ namespace codac2 for (Index c=0;c + inline IntervalMatrix left_inverse_enclosure + (const Eigen::MatrixBase& A) + { + Index N = A.rows(); + Index M = A.cols(); + assert_release(N>M); + /* we use the transpose to solve, then we correct with LEFT_INV */ + if constexpr(std::is_same_v) { + return inverse_correction(A, + (A.mid().transpose()).fullPivLu() + .solve(Matrix::Identity(M,M)).transpose()); + } else { + return inverse_correction(A, + (A.transpose()).fullPivLu() + .solve(Matrix::Identity(M,M)).transpose()); + } + } + + /** + * \brief Enclosure of a right inverse of a (full-rank) matrix + * expression, possibly an interval matrix. + * + * \pre \f$\mathbf{A}\f$ has more columns than rows. + * + * \param A A non-square matrix expression + * \return The enclosure of a right-inverse. Can have + * \f$(-\infty,\infty)\f$ coefficients if the inversion "failed", + * either because the matrix is (almost) not full-rank, or the full pivot + * LU decomposition did not use the ''correct'' pivots + * (for Interval matrix) */ + template + inline IntervalMatrix right_inverse_enclosure + (const Eigen::MatrixBase& A) + { + Index N = A.rows(); + Index M = A.cols(); + assert_release(M>N); + /* we solve directly, and transpose for the inverse_correction (?) */ + if constexpr(std::is_same_v) { + return inverse_correction (A.transpose(), + (A.mid()).fullPivLu().solve(Matrix::Identity(N,N)).transpose()) + .transpose(); + } else { + return inverse_correction (A.transpose(), + A.fullPivLu().solve(Matrix::Identity(N,N)).transpose()) + .transpose(); + } + } + /** * \brief Compute an upper bound of \f$\left([\mathbf{A}]+[\mathbf{A}]^2+[\mathbf{A}]^3+\dots\right)\f$, * with \f$[\mathbf{A}]\f$ a matrix of intervals as an "error term" (uses only bounds on coefficients). diff --git a/src/core/separators/codac2_SepPolytope.cpp b/src/core/separators/codac2_SepPolytope.cpp new file mode 100644 index 000000000..f0e65a109 --- /dev/null +++ b/src/core/separators/codac2_SepPolytope.cpp @@ -0,0 +1,30 @@ +/** + * SepPolygon.cpp + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Simon Rohou, Benoit Desrochers + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include "codac2_SepPolytope.h" + +using namespace std; + +namespace codac2 +{ + SepPolytope::SepPolytope(const Polytope& p) + : Sep(p.size()), _ctc(p) { } + + BoxPair SepPolytope::separate(const IntervalVector& x) const + { + auto x_in = x; + auto x_out = x; + _ctc.contract(x_out); + _ctc.p().contract_out_Box(x_in); + + assert((x_in | x_out) == x); + return { x_in, x_out }; + } +} diff --git a/src/core/separators/codac2_SepPolytope.h b/src/core/separators/codac2_SepPolytope.h new file mode 100644 index 000000000..049724d00 --- /dev/null +++ b/src/core/separators/codac2_SepPolytope.h @@ -0,0 +1,29 @@ +/** + * \file codac2_SepPolygon.h + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Simon Rohou, Benoit Desrochers + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include "codac2_Polytope.h" +#include "codac2_CtcPolytope.h" +#include "codac2_SepWrapper.h" + +namespace codac2 +{ + class SepPolytope : public Sep + { + public: + + SepPolytope(const Polytope& p); + BoxPair separate(const IntervalVector& x) const; + + protected: + const CtcPolytope _ctc; + + }; +} diff --git a/src/graphics/figures/codac2_Figure2D.cpp b/src/graphics/figures/codac2_Figure2D.cpp index ed2351f5c..b904fcfd4 100644 --- a/src/graphics/figures/codac2_Figure2D.cpp +++ b/src/graphics/figures/codac2_Figure2D.cpp @@ -270,6 +270,14 @@ void Figure2D::draw_zonotope(const Zonotope& z, const StyleProperties& style) output_fig->draw_polygon(vertices,style); } +void Figure2D::draw_polytope(const Polytope& P, const StyleProperties& style) +{ + assert_release(P.size() == 2); + vector w = P.vertices_2Dfacet(); + for(const auto& output_fig : _output_figures) + output_fig->draw_polygon(w,style); +} + void Figure2D::draw_parallelepiped(const Parallelepiped& p, const StyleProperties& style) { diff --git a/src/graphics/figures/codac2_Figure2D.h b/src/graphics/figures/codac2_Figure2D.h index b4d19f18e..dc328d39c 100644 --- a/src/graphics/figures/codac2_Figure2D.h +++ b/src/graphics/figures/codac2_Figure2D.h @@ -21,6 +21,7 @@ #include "codac2_PavingStyle.h" #include "codac2_Ellipsoid.h" #include "codac2_Polygon.h" +#include "codac2_Polytope.h" #include "codac2_SlicedTube.h" #include "codac2_Ctc.h" #include "codac2_Sep.h" @@ -304,13 +305,21 @@ namespace codac2 void draw_parallelepiped(const Parallelepiped& p, const StyleProperties& style = StyleProperties()); /** - * \brief Draws a zonotope z+sum_i [-1,1] A_i on the figure + * \brief draws a zonotope z+sum_i [-1,1] a_i on the figure * - * \param z Zonotope to draw (center and shape matrix) - * \param style Style of the zonotope (edge color and fill color) + * \param z zonotope to draw (center and shape matrix) + * \param style style of the zonotope (edge color and fill color) */ void draw_zonotope(const Zonotope& z, const StyleProperties& style = StyleProperties()); + /** + * \brief draws a 2D convex polytope on the figure + * + * \param P polytope to draw + * \param style style of the polytope (edge color and fill color) + */ + void draw_polytope(const Polytope& P, const StyleProperties& style = StyleProperties()); + /** * \brief Draws a pie on the figure * @@ -882,6 +891,18 @@ namespace codac2 selected_fig()->draw_ellipsoid(e,style); } + /** + * \brief draws a 2D convex polytope on the figure + * + * \param P polytope to draw + * \param style style of the polytope (edge color and fill color) + */ + static void draw_polytope(const Polytope& P, const StyleProperties& style = StyleProperties()) + { + auto_init(); + selected_fig()->draw_polytope(P,style); + } + /** * \brief Draws a trajectory on the figure * @@ -1260,4 +1281,4 @@ namespace codac2 }; } -#include "codac2_Figure2D_pave.h" \ No newline at end of file +#include "codac2_Figure2D_pave.h" diff --git a/src/graphics/figures/codac2_Figure3D.cpp b/src/graphics/figures/codac2_Figure3D.cpp index 86d7fcdc0..7ff230039 100644 --- a/src/graphics/figures/codac2_Figure3D.cpp +++ b/src/graphics/figures/codac2_Figure3D.cpp @@ -137,7 +137,7 @@ void Figure3D::draw_box(const IntervalVector& x, const StyleProperties& style) } void Figure3D::draw_zonotope(const Zonotope& z, const StyleProperties& style) { - assert_release(z.z.size() == 3); + assert_release(z.z.size() == 3); Matrix id = Matrix::Identity(3,3); this->set_style_internal(style); lock_style=true; @@ -193,6 +193,19 @@ void Figure3D::draw_zonotope(const Zonotope& z, const StyleProperties& style) { lock_style=false; } +void Figure3D::draw_polytope(const Polytope& P, const StyleProperties& style) { + assert_release(P.size() == 3); + this->set_style_internal(style); + lock_style=true; + const Vector center = Vector::zero(3); + const Matrix transfo = Matrix::Identity(3,3); + std::vector> facets3D=P.vertices_3Dfacets(); + for (const std::vector &vec : facets3D) { + this->draw_polygon(center,transfo,vec,style); + } + lock_style=false; +} + void Figure3D::draw_arrow(const Vector& c, const Matrix &A, const StyleProperties& style) { diff --git a/src/graphics/figures/codac2_Figure3D.h b/src/graphics/figures/codac2_Figure3D.h index 3887bdee1..f2450d4db 100644 --- a/src/graphics/figures/codac2_Figure3D.h +++ b/src/graphics/figures/codac2_Figure3D.h @@ -20,6 +20,7 @@ #include "codac2_IntervalMatrix.h" #include "codac2_Ellipsoid.h" #include "codac2_Parallelepiped.h" +#include "codac2_Polytope.h" namespace codac2 { @@ -123,6 +124,14 @@ namespace codac2 */ void draw_zonotope(const Zonotope& z, const StyleProperties& style = { Color::dark_gray(0.5) }); + /** + * \brief Draws a polytope on the figure + * + * \param P Polytope to draw + * \param style Style of the polytope (edge color) + */ + void draw_polytope(const Polytope& P, const StyleProperties& style = { Color::dark_gray(0.5) }); + /** * \brief Draws a box on the figure * diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5185c4be8..3e3c2b77d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -39,6 +39,7 @@ list(APPEND SRC_TESTS # listing files without extension core/contractors/codac2_tests_CtcInverseNotIn core/contractors/codac2_tests_CtcLazy core/contractors/codac2_tests_CtcPolygon + core/contractors/codac2_tests_CtcPolytope core/contractors/codac2_tests_CtcSegment core/contractors/codac2_tests_linear_ctc ../doc/manual/manual/contractors/geometric/src @@ -51,6 +52,7 @@ list(APPEND SRC_TESTS # listing files without extension core/domains/interval/codac2_tests_IntervalVector core/domains/zonotope/codac2_tests_Parallelepiped core/domains/zonotope/codac2_tests_Parallelepiped_eval + core/domains/polytope/codac2_tests_Polytope core/domains/tube/codac2_tests_TDomain core/domains/tube/codac2_tests_Slice core/domains/tube/codac2_tests_Slice_polygon @@ -78,6 +80,7 @@ list(APPEND SRC_TESTS # listing files without extension core/matrices/codac2_tests_inversion core/matrices/codac2_tests_IntFullPivLU core/matrices/codac2_tests_GaussJordan + ../doc/manual/manual/linear/src core/operators/codac2_tests_operators @@ -88,6 +91,7 @@ list(APPEND SRC_TESTS # listing files without extension core/separators/codac2_tests_SepCtcBoundary core/separators/codac2_tests_SepInverse core/separators/codac2_tests_SepPolygon + core/separators/codac2_tests_SepPolytope core/separators/codac2_tests_SepProj core/separators/codac2_tests_SepTransform diff --git a/tests/core/contractors/codac2_tests_CtcPolytope.cpp b/tests/core/contractors/codac2_tests_CtcPolytope.cpp new file mode 100644 index 000000000..ddc71cb82 --- /dev/null +++ b/tests/core/contractors/codac2_tests_CtcPolytope.cpp @@ -0,0 +1,44 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Simon Rohou + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("CtcPolytope") +{ + std::vector> facets + { { {1,1,1}, 3.0 }, + { {-1,1,1}, 3.0 }, + { {1,-1,1}, 3.0 }, + { {1,1,-1}, 3.0 }, + { {-1,-1,1}, 3.0 }, + { {-1,1,-1}, 3.0 }, + { {1,-1,-1}, 3.0 }, + { {-1,-1,-1}, 3.0 } }; + Polytope p({{-2,2.0},{-2.0,2.0},{-2.0,2.0}}, facets); + CtcPolytope c(p); + + IntervalVector x(3); + c.contract(x); + CHECK(x == IntervalVector({{-2,2},{-2,2},{-2,2}})); + + x = IntervalVector({{1.02,2.0},{1.02,2.0},{1.02,2.0}}); // possible bug + c.contract(x); + CHECK(x.is_empty()); + +// auto p_ctc = pave({{-3,3},{-3,3},{-3,3}}, c, 0.2); +// Figure3D fig_ctc("Paving contractor"); +// fig_ctc.draw_paving(p_ctc); +} + diff --git a/tests/core/contractors/codac2_tests_CtcPolytope.py b/tests/core/contractors/codac2_tests_CtcPolytope.py new file mode 100644 index 000000000..aa538bf22 --- /dev/null +++ b/tests/core/contractors/codac2_tests_CtcPolytope.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2025 +# \author Simon Rohou +# \copyright Copyright 2024 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +from codac import * + +class TestSepPolytope(unittest.TestCase): + + def test_SepPolytope(self): + + # Polytope, defined as a set of facets + facets =[ [ Row([1,1,1]), 3.0 ], + [ Row([-1,1,1]), 3.0 ], + [ Row([1,-1,1]), 3.0 ], + [ Row([1,1,-1]), 3.0 ], + [ Row([-1,-1,1]), 3.0 ], + [ Row([-1,1,-1]), 3.0 ], + [ Row([1,-1,-1]), 3.0 ], + [ Row([-1,-1,-1]), 3.0 ] ] + + p = Polytope(IntervalVector([[-2,2],[-2,2],[-2,2]]),facets,True) + c = CtcPolytope(p) + + x = IntervalVector(3) + c.contract(x) + self.assertTrue(x == IntervalVector([[-2,2],[-2,2],[-2,2]])) + + x = IntervalVector([[1.02,2.0],[1.02,2.0],[1.0,2.0]]) + c.contract(x) + self.assertTrue(x.is_empty()) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/core/domains/polytope/codac2_tests_Polytope.cpp b/tests/core/domains/polytope/codac2_tests_Polytope.cpp new file mode 100644 index 000000000..0dfef94f2 --- /dev/null +++ b/tests/core/domains/polytope/codac2_tests_Polytope.cpp @@ -0,0 +1,99 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2025 + * \author Damien Massé + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("Polytope") +{ + /* polytope with vertices + (0,0,0) , (2,0,1) , (1,0.1,0) , (0,-0.1,1), (1,2,0.5) + inequalities + -x + 20y + 82z <= 80 x - 20y - 2z <= 0 + -20.5x + 10y +z <= 0 -x + 10y - 38z <= 0 + 3.9x + y - 3.8z <= 4 x -10y - 2z <= 0 + + "redundant" inequality + -2.1 x + 2y + 4.2z <= 4 */ + std::vector> facets + { { {-1,20,82}, 80 }, + { {-20.5,10,1} ,0 }, + { {3.9,1,-3.8} ,4 }, + { {1,-20,-2} ,0 }, + { {-1,10,-38} ,0 }, + { {1,-10,-2} ,0 } }; + + SECTION("Definition, box, contains") + { + Polytope p(IntervalVector({{-4,4},{-4,4},{-4,4}}), facets, true); + + CHECK(!p.is_empty()); + CHECK(!p.is_flat()); + CHECK(Approx(p.box(),1e-8) == IntervalVector({{0.,2.},{-0.1,2.},{0.,1.}})); + CHECK(IntervalVector({{0.,2.},{-0.1,2.},{0.,1.}}).is_subset(p.box())); + CHECK(Approx(p.bound_row(Row({-2.1,2,4.2})),1e-8)==4.0); + CHECK(Approx(p.bound_row(Row({1.0,1.0,1.0})),1e-8)==3.5); + CHECK(p.contains(IntervalVector({1.0,1.0,0.5}))==BoolInterval::TRUE); + CHECK(p.contains(IntervalVector({1.0,1.0,61.0/82.0}))==BoolInterval::UNKNOWN); + CHECK(p.contains(IntervalVector({1.0,1.0,1.0}))==BoolInterval::FALSE); + } + + SECTION("Subset") + { + Polytope p(IntervalVector({{-4,4},{-4,4},{-4,4}}), facets, true); + + Polytope q(IntervalVector({{0.2,0.4},{0.0,0.3},{0.3,0.5}})); + CHECK(q.is_subset(p)); + Polytope r(3,true); /* empty polyhedron */ + CHECK(r.is_empty()); + CHECK(r.is_subset(q)); + } + + SECTION("Transformation") + { + Polytope p(IntervalVector({{-4,4},{-4,4},{-4,4}}), facets, true); + Polytope q = p; + q.inflate(0.5); + CHECK(p.is_subset(q)); + Polytope r = p; + r.inflate_ball(0.5); + CHECK(r.is_subset(q)); + Polytope s = p; + s.unflat(2,0.5); + CHECK(s.is_subset(r)); + Polytope t = p; + t.homothety(IntervalVector({0,0,-0.5}),2); + CHECK(Approx(t.bound_row(Row({1.0,1.0,1.0})),1e-8)==7.5); + Polytope u=t; + u.meet_with_polytope(p); + CHECK(u.is_subset(t)); + CHECK(u.is_subset(p)); + Polytope v = Polytope::union_of_polytopes({ p,t }); + CHECK(p.is_subset(v)); + CHECK(t.is_subset(v)); + t.inflate(0.1); + CHECK(!t.is_subset(v)); + /* transformation : rotation of pi/3 degrees (reverse), center (1,0,0) */ + IntervalMatrix M { { cos(PI/3) , sin(PI/3) , 0 }, + { -sin(PI/3), cos(PI/3) , 0 }, + { 0 , 0 , 1 } }; + IntervalVector P { 1-cos(PI/3) , sin(PI/3) , 0 }; + Polytope w = p.reverse_affine_transform(M,P, + IntervalVector({{-4,4},{-4,4},{-4,4}})); + CHECK(Approx(w.box(),1e-5) == + IntervalVector({{1.0-2*sin(PI/3),1.5}, + {-sin(PI/3)-0.1*cos(PI/3),1.}, + {0.,1.}})); + } + +} + diff --git a/tests/core/domains/polytope/codac2_tests_Polytope.py b/tests/core/domains/polytope/codac2_tests_Polytope.py new file mode 100644 index 000000000..728744d18 --- /dev/null +++ b/tests/core/domains/polytope/codac2_tests_Polytope.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2025 +# \author Damien Massé +# \copyright Copyright 2024 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +from codac import * +import sys +import math + + +class TestPolytope(unittest.TestCase): + facets = [[Row([-1,20,82]),80], + [Row([-20.5,10,1]),0], + [Row([3.9,1,-3.8]),4], + [Row([1,-20,-2]),0], + [Row([-1,10,-38]),0], + [Row([1,-10,-2]),0]] + +## polytope with vertices +## (0,0,0) , (2,0,1) , (1,0.1,0) , (0,-0.1,1), (1,2,0.5) +## inequalities +## -x + 20y + 82z <= 80 x - 20y - 2z <= 0 +## -20.5x + 10y +z <= 0 -x + 10y - 38z <= 0 +## 3.9x + y - 3.8z <= 4 x -10y - 2z <= 0 +## + "redundant" inequality +## -2.1 x + 2y + 4.2z <= 4 + + def test_polytope(self): + p = Polytope(IntervalVector([[-4,4],[-4,4],[-4,4]]),self.facets,True) + self.assertFalse(p.is_empty()) + self.assertFalse(p.is_flat()) + self.assertTrue(Approx(p.box(),1e-8)==IntervalVector([[0.,2.],[-0.1,2.],[0.,1.]])) + self.assertTrue(Approx(p.bound_row(Row([-2.1,2.0,4.2])),1e-8)==4.0) + self.assertTrue(Approx(p.bound_row(Row([1.0,1.0,1.0])),1e-8)==3.5) + self.assertTrue(p.contains(IntervalVector([1.0,1.0,0.5]))==BoolInterval.TRUE) +# self.assertTrue(p.contains(IntervalVector([1.0,1.0,61.0/82.0]))==BoolInterval.UNKNOWN) + self.assertTrue(p.contains(IntervalVector([1.0,1.0,1.0]))==BoolInterval.FALSE) + + def test_polytope_subset(self): + p = Polytope(IntervalVector([[-4,4],[-4,4],[-4,4]]),self.facets,True) + q = Polytope(IntervalVector([[0.2,0.4],[0.0,0.3],[0.3,0.5]])) + self.assertTrue(q.is_subset(p)) + r = Polytope(3,True) # empty polyhedron + self.assertTrue(r.is_empty()) + self.assertTrue(r.is_subset(q)) + + def test_polytope_transformation(self): + p = Polytope(IntervalVector([[-4,4],[-4,4],[-4,4]]),self.facets,True) + q = Polytope(p) # or q = Polytope() q.assign(p) + q.inflate(0.5) + self.assertTrue(p.is_subset(q)) + r = Polytope(p) + r.inflate_ball(0.5) + self.assertTrue(r.is_subset(q)) + self.assertFalse(q.is_subset(r)) + s = Polytope(p) + s.unflat(2,0.5) + self.assertTrue(s.is_subset(r)) + t = Polytope(p) + t.homothety(IntervalVector([0,0,-0.5]),2) + self.assertTrue(Approx(t.bound_row(Row([1.0,1.0,1.0])),1e-8)==7.5) + u = Polytope(t) + u.meet_with_polytope(p) + self.assertTrue(u.is_subset(p)) + self.assertTrue(u.is_subset(t)) + v = Polytope.union_of_polytopes([ p,t ]) + self.assertTrue(p.is_subset(v)) + self.assertTrue(t.is_subset(v)) + t.inflate(0.1) + self.assertFalse(t.is_subset(v)) +## transformation : rotation of pi/3 degrees (reverse), center (1,0,0) + M = IntervalMatrix([[cos(PI/3),sin(PI/3),0],[-sin(PI/3),cos(PI/3),0], + [0,0,1]]) + P = IntervalVector([1-cos(PI/3),sin(PI/3),0]) + w = p.reverse_affine_transform(M,P,IntervalVector([[-4,4],[-4,4],[-4,4]])) + self.assertTrue(Approx(w.box(),1e-8)== + IntervalVector([[1.0-2*math.sin(PI/3),1.5], + [-math.sin(PI/3)-0.1*math.cos(PI/3),1], + [0,1]])) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/core/matrices/codac2_tests_inversion.cpp b/tests/core/matrices/codac2_tests_inversion.cpp index d5b553d49..5db990bc9 100644 --- a/tests/core/matrices/codac2_tests_inversion.cpp +++ b/tests/core/matrices/codac2_tests_inversion.cpp @@ -14,10 +14,12 @@ #include #include +#include + using namespace std; using namespace codac2; -TEST_CASE("Matrix") +TEST_CASE("Square matrices") { IntervalMatrix x({ { {0.0,0.0}, {-0.1,-0.1}, {0.2,0.2} }, @@ -52,3 +54,36 @@ TEST_CASE("Matrix") y = inverse_enclosure(w); CHECK((w.template cast()*y).contains(Matrix::Identity(3,3))); } + +TEST_CASE("Non-square matrices") +{ + IntervalMatrix x({ + { {0.0,0.0}, {-0.1,-0.1}, {0.2,0.2}, {0.1,0.1} }, + { {0.0,0.0}, {-0.2,-0.2}, {0.1,0.1}, {0.3,0.3} }, + { {0.1,0.1}, {-0.1,-0.1}, {0.1,0.1}, {-0.1,-0.1} } + }); + + IntervalMatrix z = right_inverse_enclosure(x); + CHECK(Approx(x*z,1e-8)==IntervalMatrix::Identity(3,3)); + + Matrix u( { + { 1,2,3 }, + { 1,3,5 }, + { 2,6,10 }, + { 3,4,5 }, + { 2,1,3 }, + { 6,2,4 } + }); + + IntervalMatrix v=left_inverse_enclosure(u); + CHECK((v*u.template cast()).contains(Matrix::Identity(3,3))); + + Matrix w({ + { 1, 2, 0, 3, 1, 5 }, + { 0, 4, 1, 8, 2, 4 }, + { 0, 1, 0, 4, 5, 1 } + }); + + IntervalMatrix y = right_inverse_enclosure(w); + CHECK((w.template cast()*y).contains(Matrix::Identity(3,3))); +} diff --git a/tests/core/separators/codac2_tests_SepPolytope.cpp b/tests/core/separators/codac2_tests_SepPolytope.cpp new file mode 100644 index 000000000..b9f08b859 --- /dev/null +++ b/tests/core/separators/codac2_tests_SepPolytope.cpp @@ -0,0 +1,44 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Simon Rohou + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include + +using namespace std; +using namespace codac2; + +TEST_CASE("SepPolytope") +{ + std::vector> facets + { { {1,1,1}, 3.0 }, + { {-1,1,1}, 3.0 }, + { {1,-1,1}, 3.0 }, + { {1,1,-1}, 3.0 }, + { {-1,-1,1}, 3.0 }, + { {-1,1,-1}, 3.0 }, + { {1,-1,-1}, 3.0 }, + { {-1,-1,-1}, 3.0 } }; + Polytope p({{-2,2.0},{-2.0,2.0},{-2.0,2.0}}, facets); + SepPolytope s(p); + + IntervalVector x(3); + auto xs = s.separate(x); + CHECK(xs.outer == IntervalVector({{-2,2},{-2,2},{-2,2}})); + CHECK(xs.inner == IntervalVector(3)); + + x = IntervalVector({{1.02,2.0},{1.02,2.0},{1.02,2.0}}); + xs = s.separate(x); + CHECK(xs.inner == x); + CHECK(xs.outer.is_empty()); + +} + diff --git a/tests/core/separators/codac2_tests_SepPolytope.py b/tests/core/separators/codac2_tests_SepPolytope.py new file mode 100644 index 000000000..4b423d9f5 --- /dev/null +++ b/tests/core/separators/codac2_tests_SepPolytope.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2025 +# \author Simon Rohou +# \copyright Copyright 2024 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +from codac import * + +class TestSepPolytope(unittest.TestCase): + + def test_SepPolytope(self): + + # Polytope, defined as a set of facets + facets =[ [ Row([1,1,1]), 3.0 ], + [ Row([-1,1,1]), 3.0 ], + [ Row([1,-1,1]), 3.0 ], + [ Row([1,1,-1]), 3.0 ], + [ Row([-1,-1,1]), 3.0 ], + [ Row([-1,1,-1]), 3.0 ], + [ Row([1,-1,-1]), 3.0 ], + [ Row([-1,-1,-1]), 3.0 ] ] + + p = Polytope(IntervalVector([[-2,2],[-2,2],[-2,2]]),facets,True) + s = SepPolytope(p) + + x = IntervalVector(3) + inner,outer = s.separate(x) + self.assertTrue(inner == IntervalVector(3)) + self.assertTrue(outer == IntervalVector([[-2,2],[-2,2],[-2,2]])) + + x = IntervalVector([[1.02,2.0],[1.02,2.0],[1.0,2.0]]) + inner,outer = s.separate(x) + self.assertTrue(inner == x) + self.assertTrue(outer.is_empty()) + +if __name__ == '__main__': + unittest.main()