Merge branch 'MLouedec-codac2_ellipsoids' into codac2_ellipsoids

examples/ellipsoid_dev/CMakeLists.txt

@@ -0,0 +1,44 @@
+# ==================================================================
+#  codac / ellipsoid example - cmake configuration file
+# ==================================================================
+
+  cmake_minimum_required(VERSION 3.0.2)
+  project(ellipsoid_dev 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 Eigen3
+
+  # In case you installed Eigen3 in a local directory, you need
+  # to specify its path with the CMAKE_PREFIX_PATH option, e.g.
+  # set(CMAKE_PREFIX_PATH "~/eigen/build_install")
+
+  find_package(Eigen3 3.4 REQUIRED NO_MODULE)
+  message(STATUS "Found Eigen3 version ${Eigen3_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 "~/Documents/Code_these/codac/build_install")
+
+  find_package(CODAC REQUIRED)
+  message(STATUS "Found Codac version ${CODAC_VERSION}")
+
+# Compilation
+
+  add_executable(${PROJECT_NAME} main.cpp)
+  target_compile_options(${PROJECT_NAME} PUBLIC ${CODAC_CXX_FLAGS})
+  target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${CODAC_INCLUDE_DIRS})
+  target_link_libraries(${PROJECT_NAME} PUBLIC ${CODAC_LIBRARIES} Ibex::ibex Eigen3::Eigen)

examples/ellipsoid_dev/main.cpp

@@ -0,0 +1,123 @@
+#include <codac>
+
+using namespace std;
+using namespace codac2;
+
+int main() {
+    // ----------------------------------------------------------
+    // ellipsoid projections
+    // ----------------------------------------------------------
+
+    Vector mu4({1., 0., 0.});
+    Matrix G4({{1.,  0.5, 0.},
+               {0.5, 2.,  0.2},
+               {0.,  0.2, 3.}});
+    Ellipsoid e4(mu4, G4);
+
+    Matrix G5 = 0.7 * G4;
+    Ellipsoid e5(mu4, G5);
+
+    Matrix G6({{2., 0.,  0.5},
+               {0., 1.,  0.2},
+               {0., 0.2, 3.}});
+    Ellipsoid e6(mu4, G6);
+
+    Figure2D fig2("Projected ellipsoid xy", GraphicOutput::VIBES);
+    Figure2D fig3("Projected ellipsoid yz", GraphicOutput::VIBES);
+    Figure2D fig4("Projected ellipsoid xz", GraphicOutput::VIBES);
+
+    fig2.set_window_properties({700, 100}, {500, 500});
+    fig3.set_window_properties({1200, 100}, {500, 500});
+    fig4.set_window_properties({0, 100}, {500, 500});
+
+    fig2.set_axes(axis(0, {-3, 3}), axis(1, {-3, 3}));
+    fig3.set_axes(axis(1, {-3, 3}), axis(2, {-3, 3}));
+    fig4.set_axes(axis(0, {-3, 3}), axis(2, {-3, 3}));
+
+    fig2.draw_ellipsoid(e4, {Color::blue(), Color::blue(0.3)});
+    fig3.draw_ellipsoid(e4, {Color::blue(), Color::blue(0.3)});
+    fig4.draw_ellipsoid(e4, {Color::blue(), Color::blue(0.3)});
+
+    fig2.draw_ellipsoid(e5, {Color::red(), Color::red(0.3)});
+    fig3.draw_ellipsoid(e5, {Color::red(), Color::red(0.3)});
+    fig4.draw_ellipsoid(e5, {Color::red(), Color::red(0.3)});
+
+    fig2.draw_ellipsoid(e6, {Color::green(), Color::green(0.3)});
+    fig3.draw_ellipsoid(e6, {Color::green(), Color::green(0.3)});
+    fig4.draw_ellipsoid(e6, {Color::green(), Color::green(0.3)});
+
+    // ----------------------------------------------------------
+    // inclusion tests
+    // ----------------------------------------------------------
+
+    cout << "\nInclusion test e5 in e4: " << e5.is_concentric_subset(e4) << endl;
+
+    cout << "Inclusion test e4 in e5: " << e4.is_concentric_subset(e5) << endl;
+
+    cout << "Inclusion test e4 in e6: " << e6.is_concentric_subset(e4) << endl;
+
+    cout << "Inclusion test e5 in e6: " << e5.is_concentric_subset(e6) << endl;
+
+    // ----------------------------------------------------------
+    // test the non inclusion (e4 in e6)
+    // ----------------------------------------------------------
+
+    // step 1 identify the direction of non inclusion
+    Eigen::MatrixXd X = e4.G._e;
+    Eigen::MatrixXd Y = e6.G._e;
+    Eigen::MatrixXd Y_inv = Y.inverse();
+    Eigen::MatrixXd D = Eigen::MatrixXd::Identity(3,3)- X.transpose()*Y_inv.transpose()*Y_inv*X;
+    cout << "D = " << D << endl;
+
+    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(D);
+    Eigen::VectorXd eigenvalues = es.eigenvalues();
+
+    Eigen::VectorXd v;
+    double min_eigenvalue = MAXFLOAT;
+    for (int i = 0; i < eigenvalues.size(); i++){
+            if (eigenvalues(i) < min_eigenvalue){
+                min_eigenvalue = eigenvalues(i);
+                v = es.eigenvectors().col(i); // the vector associated to the smallest eigenvalue
+            }
+    }
+    cout << " the min eigenvalue of D is " << min_eigenvalue << endl;
+    cout << " the vector associated to the min eigenvalue is " << v << endl;
+
+    // select the less included point and draw it on the figures
+    Eigen::VectorXd x = e4.mu._e+e4.G._e*v;
+    cout << "the point " << x << " is the less included point" << endl;
+    fig2.draw_point(x, {Color::black()});
+    fig3.draw_point(x, {Color::black()});
+    fig4.draw_point(x, {Color::black()});
+
+//    IntervalVector x6(x);
+//    x6.inflate(0.001);
+//    fig2.draw_box(x6, {Color::black(), Color::black()});
+
+    // now let us prove by intervals that x is out of e6
+    // by proving (x-mu6).T@G6^(-T)@G6^(-1)@(x-mu6) > 1
+
+    IntervalVector w(x);
+    w -= IntervalVector(e6.mu);
+
+    IntervalMatrix Q(Y_inv);
+    Q = Q.transpose()*Q;
+
+    IntervalMatrix res = w.transpose()*Q*w;
+    cout << "the result of the inclusion test is " << res << endl;
+
+    if(res(0,0).lb() > 1.) {
+        cout << "the point is out of the ellipsoid" << endl;
+        cout << "the non inclusion is guaranteed" << endl;
+    }
+
+    // can we also prove that x is in e4? -> we should take the point not on the boundary or
+    // have guaranteed operations previously
+    Eigen::MatrixXd X_inv = X.inverse();
+    IntervalMatrix Qx(X_inv);
+    Qx = Qx.transpose()*Qx;
+
+    IntervalMatrix res2 = w.transpose()*Qx*w;
+    cout << "the result of the inclusion test in e4 is " << res2 << endl;
+
+}

src/core/domains/ellipsoid/codac2_Ellipsoid.cpp

@@ -88,6 +88,28 @@ namespace codac2 {
         return BoolInterval::TRUE;
     }
 
+    void Ellipsoid::projection2D(const Vector& d, const Vector& v, const Vector& u)
+    {
+        // from [Pope S. B. - Algorithms for Ellipsoids - 2008]
+        assert( d.size() == v.size());
+        assert( d.size() == u.size());
+        assert( d.size() == this->size());
+        assert( (v._e.transpose()*u._e).norm() == 0); // u & v orthogonal
+
+        // Normalized Projection matrix
+        // the plane (d,u,v) is also the affine plan {x|x=d+Tt} with T = [u,v]
+        Matrix T = Matrix(this->size(),2);
+        T.col(0) = v/v.norm();
+        T.col(1) = u/u.norm();
+
+        auto TTG = T._e.transpose() * this->G._e;
+        Eigen::BDCSVD<Eigen::MatrixXd> bdcsvd(TTG, Eigen::ComputeFullU);
+        Matrix U(bdcsvd.matrixU());
+        Matrix E((Eigen::MatrixXd) bdcsvd.singularValues().asDiagonal());
+        this->G = U._e * E._e;
+        this->mu = T._e.transpose() * (d._e + T._e * T._e.transpose() * (this->mu._e - d._e));
+    }
+
     Ellipsoid operator+(const Ellipsoid &e1, const Ellipsoid &e2) {
         assert_release(e1.size() == e2.size());
 

src/core/domains/ellipsoid/codac2_Ellipsoid.h

@@ -90,6 +90,14 @@ namespace codac2
        */
       BoolInterval is_concentric_subset(const Ellipsoid& e) const;
 
+      /**
+       * \brief Project the ellipsoid on the 2d plane (d,v,u)
+       * \param d a point on the plane
+       * \param v a vector of the plane
+       * \param u a other vector of the plane, orthogonal to v
+       */
+      void projection2D(const Vector& d, const Vector& v, const Vector& u);
+
     public:
 
       Vector mu; ///< midpoint vector

src/graphics/figures/codac2_Figure2D.cpp

@@ -203,22 +203,31 @@ void Figure2D::draw_ellipsoid(const Ellipsoid &e, const StyleProperties &s) {
         Vector mu_draw(2);
         // 2d projection of the ellipsoid
         if (e.size() > 2) {
-            // affine space of the projection
-            Vector d(Eigen::VectorXd::Zero(e.mu.nb_rows()));
-            Matrix T(Eigen::MatrixXd::Zero(e.G.nb_rows(), 2));
-            T(output_fig->i(), 0) = 1;
-            T(output_fig->j(), 1) = 1;
-
-            // project ellipsoid E(mu,Q) = {x in R^n | (x-mu).T*G.{-T}*G^{-1}*(x-mu)<1}
-            // on the affine plan A = {x|x=d+Tt} [Pope -2008]
-            // reduce the dimensions of mu and Q
-
-            auto TTG = T._e.transpose() * e.G._e;
-            Eigen::BDCSVD<Eigen::MatrixXd> bdcsvd(TTG, Eigen::ComputeFullU);
-            Matrix U(bdcsvd.matrixU());
-            Matrix E((Eigen::MatrixXd) bdcsvd.singularValues().asDiagonal());
-            G_draw = U._e * E._e;
-            mu_draw = T._e.transpose() * (d._e + T._e * T._e.transpose() * (e.mu._e - d._e));
+            Ellipsoid ep(e);
+            Vector v = Eigen::VectorXd::Zero(e.size());
+            v[output_fig->i()] = 1;
+            Vector u = Eigen::VectorXd::Zero(e.size());
+            u[output_fig->j()] = 1;
+            ep.projection2D(Eigen::VectorXd::Zero(e.size()),v,u);
+            G_draw = ep.G;
+            mu_draw = ep.mu;
+
+//            // affine space of the projection
+//            Vector d(Eigen::VectorXd::Zero(e.mu.nb_rows()));
+//            Matrix T(Eigen::MatrixXd::Zero(e.G.nb_rows(), 2));
+//            T(output_fig->i(), 0) = 1;
+//            T(output_fig->j(), 1) = 1;
+//
+//            // project ellipsoid E(mu,Q) = {x in R^n | (x-mu).T*G.{-T}*G^{-1}*(x-mu)<1}
+//            // on the affine plan A = {x|x=d+Tt} [Pope -2008]
+//            // reduce the dimensions of mu and Q
+//
+//            auto TTG = T._e.transpose() * e.G._e;
+//            Eigen::BDCSVD<Eigen::MatrixXd> bdcsvd(TTG, Eigen::ComputeFullU);
+//            Matrix U(bdcsvd.matrixU());
+//            Matrix E((Eigen::MatrixXd) bdcsvd.singularValues().asDiagonal());
+//            G_draw = U._e * E._e;
+//            mu_draw = T._e.transpose() * (d._e + T._e * T._e.transpose() * (e.mu._e - d._e));
         } else {
             G_draw = e.G;
             mu_draw = e.mu;