diff --git a/examples/ellipsoid_dev/CMakeLists.txt b/examples/ellipsoid_dev/CMakeLists.txt new file mode 100644 index 00000000..a9eeab77 --- /dev/null +++ b/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) diff --git a/examples/ellipsoid_dev/main.cpp b/examples/ellipsoid_dev/main.cpp new file mode 100644 index 00000000..f8554893 --- /dev/null +++ b/examples/ellipsoid_dev/main.cpp @@ -0,0 +1,123 @@ +#include + +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 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; + +} diff --git a/src/core/domains/ellipsoid/codac2_Ellipsoid.cpp b/src/core/domains/ellipsoid/codac2_Ellipsoid.cpp index 1cf11247..04d0afe7 100644 --- a/src/core/domains/ellipsoid/codac2_Ellipsoid.cpp +++ b/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 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()); diff --git a/src/core/domains/ellipsoid/codac2_Ellipsoid.h b/src/core/domains/ellipsoid/codac2_Ellipsoid.h index d839b8b5..f2003619 100644 --- a/src/core/domains/ellipsoid/codac2_Ellipsoid.h +++ b/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 diff --git a/src/graphics/figures/codac2_Figure2D.cpp b/src/graphics/figures/codac2_Figure2D.cpp index 3df8eafa..c9c78181 100644 --- a/src/graphics/figures/codac2_Figure2D.cpp +++ b/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 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 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;