diff options
| author | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-12-07 09:39:53 +0000 |
|---|---|---|
| committer | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-12-07 09:39:53 +0000 |
| commit | 33c51358238382335caf892bbc24759c8aac59a0 (patch) | |
| tree | 499e768570391c30af22eac4443667604fa717d6 /src | |
| parent | da39f7cd8a0db5d7fa13c9c87f8fc3e038c10d01 (diff) | |
| parent | c8c2f91db880218bb7ab275fbadda53a23f88d35 (diff) | |
Changes piled up for months
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@932 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 0447901e608890eb607456fd12f3ea53547b8f10
Diffstat (limited to 'src')
46 files changed, 4367 insertions, 1214 deletions
diff --git a/src/Alpha_shapes/example/CMakeLists.txt b/src/Alpha_shapes/example/CMakeLists.txt index fb94ca05..753238a5 100644 --- a/src/Alpha_shapes/example/CMakeLists.txt +++ b/src/Alpha_shapes/example/CMakeLists.txt @@ -19,7 +19,8 @@ if(CGAL_FOUND) target_link_libraries(dtoffrw ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) - #add_definitions(-DDEBUG_TRACES) + # uncomment to display debug traces + # add_definitions(-DDEBUG_TRACES) add_executable ( stfromdt Simplex_tree_from_delaunay_triangulation.cpp ) target_link_libraries(stfromdt ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) else() diff --git a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp index a90704b6..b4c32321 100644 --- a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp +++ b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp @@ -20,7 +20,7 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#define BOOST_TEST_MODULE alpha_shapes +#define BOOST_TEST_MODULE alpha_shapes test #include <boost/test/included/unit_test.hpp> #include <boost/system/error_code.hpp> #include <boost/chrono/thread_clock.hpp> diff --git a/src/Alpha_shapes/test/CMakeLists.txt b/src/Alpha_shapes/test/CMakeLists.txt index a48c1a8f..e0d33827 100644 --- a/src/Alpha_shapes/test/CMakeLists.txt +++ b/src/Alpha_shapes/test/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 2.6) -project(GUDHIAlphaShapesTest) +project(GUDHIAlphaShapesUT) # need CGAL 4.6 # cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../.. @@ -15,10 +15,28 @@ if(CGAL_FOUND) include( ${EIGEN3_USE_FILE} ) include_directories (BEFORE "../../include") - add_definitions(-DDEBUG_TRACES) - add_executable ( AlphaShapesUnitTest Alpha_shapes_unit_test.cpp ) - target_link_libraries(AlphaShapesUnitTest ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) - add_test(AlphaShapesUnitTest ${CMAKE_CURRENT_BINARY_DIR}/AlphaShapesUnitTest) + if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") + endif() + if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") + endif() + + # uncomment to display debug traces + # add_definitions(-DDEBUG_TRACES) + add_executable ( AlphaShapesUT Alpha_shapes_unit_test.cpp ) + target_link_libraries(AlphaShapesUT ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + add_test(NAME AlphaShapesUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/AlphaShapesUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/AlphaShapesUT.xml --log_level=test_suite --report_level=no) + else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") @@ -28,4 +46,3 @@ if(CGAL_FOUND) endif () endif() -cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Alpha_shapes/include/gudhi") diff --git a/src/Alpha_shapes/test/README b/src/Alpha_shapes/test/README index 244a2b84..cddd46ca 100644 --- a/src/Alpha_shapes/test/README +++ b/src/Alpha_shapes/test/README @@ -1,12 +1,14 @@ To compile: *********** +cd /path-to-gudhi/ cmake . +cd /path-to-test/ make To launch with details: *********************** -./AlphaShapesUnitTest --report_level=detailed --log_level=all +./AlphaShapesUT --report_level=detailed --log_level=all ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h index 7e278209..73ad940b 100644 --- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h @@ -34,7 +34,7 @@ namespace Gudhi { namespace bottleneck { // Diagram_point is the type of the persistence diagram's points -typedef typename std::pair<double, double> Diagram_point; +typedef std::pair<double, double> Diagram_point; // Return the used index for encoding none of the points int null_point_index(); @@ -81,7 +81,7 @@ Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& dia swap(u, v); } -Persistence_diagrams_graph::Persistence_diagrams_graph::Persistence_diagrams_graph() +Persistence_diagrams_graph::Persistence_diagrams_graph() : u(), v() { } inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const { @@ -108,7 +108,7 @@ inline double Persistence_diagrams_graph::distance(int u_point_index, int v_poin return 0; Diagram_point p_u = get_u_point(u_point_index); Diagram_point p_v = get_v_point(v_point_index); - return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second)); + return (std::max)(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second)); } inline int Persistence_diagrams_graph::size() const { diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt index 7044372e..3dfd80cd 100644 --- a/src/Bottleneck/test/CMakeLists.txt +++ b/src/Bottleneck/test/CMakeLists.txt @@ -1,21 +1,25 @@ cmake_minimum_required(VERSION 2.6) -project(GUDHIBottleneckUnitTest) +project(GUDHIBottleneckUT) -if(NOT MSVC) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") endif() -add_executable ( BottleneckUnitTest bottleneck_unit_test.cpp ) -target_link_libraries(BottleneckUnitTest ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +add_executable ( BottleneckUT bottleneck_unit_test.cpp ) +target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests -add_test(BottleneckUnitTest ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUnitTest) - -if (LCOV_PATH) - # Lcov code coverage of unitary test - add_test(src/Bottleneck/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Bottleneck) -endif() +add_test(NAME BottleneckUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BottleneckUT.xml --log_level=test_suite --report_level=no) -cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Bottleneck/include/gudhi") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 362ff7a8..70fc9a45 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,25 +10,32 @@ if (NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE "Release") endif() if(MSVC) - SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4668 /wd4311 /wd4800 /wd4820 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018") + SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4668 /wd4311 /wd4800 /wd4820 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018") else() - list(APPEND CMAKE_CXX_FLAGS "-std=c++11") -endif() + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") +endif() -# BOOST ISSUE result_of vs C++11 -add_definitions(-DBOOST_RESULT_OF_USE_DECLTYPE) +set(Boost_USE_STATIC_LIBS ON) +set(Boost_USE_MULTITHREADED ON) +set(Boost_USE_STATIC_RUNTIME OFF) -find_package(Boost) +find_package(Boost) find_package(GMP) if(GMP_FOUND) - find_package(GMPXX) + find_package(GMPXX) endif() + find_package(CGAL) if(NOT Boost_FOUND) message(FATAL_ERROR "NOTICE: This demo requires Boost and will not be compiled.") else() - INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS}) + # BOOST ISSUE result_of vs C++11 + add_definitions(-DBOOST_RESULT_OF_USE_DECLTYPE) + # BOOST ISSUE with Libraries name resolution under Windows + add_definitions(-DBOOST_ALL_NO_LIB) + + INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS}) LINK_DIRECTORIES(${Boost_LIBRARY_DIRS}) include_directories(include/) @@ -41,4 +48,10 @@ else() add_subdirectory(example/Alpha_shapes) add_subdirectory(example/Bottleneck) + # GudhUI + add_subdirectory(GudhUI) + + # data points generator + add_subdirectory(data/points/generator) + endif() diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp index 5b178ff9..a41f65aa 100644 --- a/src/Contraction/example/Garland_heckbert.cpp +++ b/src/Contraction/example/Garland_heckbert.cpp @@ -7,7 +7,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA Sophia Antipolis-M�diterran�e (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -42,72 +42,70 @@ using namespace Gudhi; using namespace skbl; using namespace contraction; - -struct Geometry_trait{ - typedef Point_d Point; +struct Geometry_trait { + typedef Point_d Point; }; /** * The vertex stored in the complex contains a quadric. */ struct Garland_heckbert_traits : public Skeleton_blocker_simple_geometric_traits<Geometry_trait> { + public: -public: - struct Garland_heckbert_vertex : public Simple_geometric_vertex{ - Error_quadric<Geometry_trait::Point> quadric; - }; - typedef Garland_heckbert_vertex Graph_vertex; + struct Garland_heckbert_vertex : public Simple_geometric_vertex { + Error_quadric<Geometry_trait::Point> quadric; + }; + typedef Garland_heckbert_vertex Graph_vertex; }; typedef Skeleton_blocker_geometric_complex< Garland_heckbert_traits > Complex; typedef Edge_profile<Complex> EdgeProfile; typedef Skeleton_blocker_contractor<Complex> Complex_contractor; - /** * How the new vertex is placed after an edge collapse : here it is placed at * the point minimizing the cost of the quadric. */ -class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile>{ - Complex& complex_; -public: - typedef typename Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; - - GH_placement(Complex& complex):complex_(complex){} - - Placement_type operator()(const EdgeProfile& profile) const override{ - auto sum_quad(profile.v0().quadric); - sum_quad += profile.v1().quadric; - - boost::optional<Point> min_quadric_pt(sum_quad.min_cost()); - if (min_quadric_pt) - return Placement_type(*min_quadric_pt); - else - return profile.p0(); - } +class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> { + Complex& complex_; + public: + typedef Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; + + GH_placement(Complex& complex) : complex_(complex) { } + + Placement_type operator()(const EdgeProfile& profile) const override { + auto sum_quad(profile.v0().quadric); + sum_quad += profile.v1().quadric; + + boost::optional<Point> min_quadric_pt(sum_quad.min_cost()); + if (min_quadric_pt) + return Placement_type(*min_quadric_pt); + else + return profile.p0(); + } }; /** * How much cost an edge collapse : here the costs is given by a quadric * which expresses a squared distances with triangles planes. */ -class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile>{ - Complex& complex_; -public: - - typedef typename Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; - - GH_cost(Complex& complex):complex_(complex){} - - Cost_type operator()( EdgeProfile const& profile, boost::optional<Point> const& new_point ) const override { - Cost_type res; - if (new_point){ - auto sum_quad(profile.v0().quadric); - sum_quad += profile.v1().quadric; - res = sum_quad.cost(*new_point); - } - return res; - } +class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> { + Complex& complex_; + public: + + typedef Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; + + GH_cost(Complex& complex) : complex_(complex) { } + + Cost_type operator()(EdgeProfile const& profile, boost::optional<Point> const& new_point) const override { + Cost_type res; + if (new_point) { + auto sum_quad(profile.v0().quadric); + sum_quad += profile.v1().quadric; + res = sum_quad.cost(*new_point); + } + return res; + } }; /** @@ -115,79 +113,80 @@ public: * Here we initializes the quadrics of every vertex at the on_started call back * and we update them when contracting an edge (the quadric become the sum of both quadrics). */ -class GH_visitor: public Gudhi::contraction::Contraction_visitor<EdgeProfile> { - Complex& complex_; -public: - GH_visitor(Complex& complex):complex_(complex){} - - //Compute quadrics for every vertex v - //The quadric of v consists in the sum of quadric - //of every triangles passing through v weighted by its area - void on_started(Complex & complex) override{ - for(auto v : complex.vertex_range()){ - auto & quadric_v(complex[v].quadric); - for(auto t : complex.triangle_range(v)){ - auto t_it = t.begin(); - const auto& p0(complex.point(*t_it++)); - const auto& p1(complex.point(*t_it++)); - const auto& p2(complex.point(*t_it++)); - quadric_v+=Error_quadric<Point>(p0,p1,p2); - } - } - } - - /** - * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement. - */ - void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement) - override{ - profile.v0().quadric += profile.v1().quadric; - } +class GH_visitor : public Gudhi::contraction::Contraction_visitor<EdgeProfile> { + Complex& complex_; + public: + + GH_visitor(Complex& complex) : complex_(complex) { } + + //Compute quadrics for every vertex v + //The quadric of v consists in the sum of quadric + //of every triangles passing through v weighted by its area + + void on_started(Complex & complex) override { + for (auto v : complex.vertex_range()) { + auto & quadric_v(complex[v].quadric); + for (auto t : complex.triangle_range(v)) { + auto t_it = t.begin(); + const auto& p0(complex.point(*t_it++)); + const auto& p1(complex.point(*t_it++)); + const auto& p2(complex.point(*t_it++)); + quadric_v += Error_quadric<Point>(p0, p1, p2); + } + } + } + + /** + * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement. + */ + void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement) + override { + profile.v0().quadric += profile.v1().quadric; + } }; +int main(int argc, char *argv[]) { + if (argc != 4) { + std::cerr << "Usage " << argv[0] << " input.off output.off N to load the file input.off, contract N edges and save the result to output.off.\n"; + return EXIT_FAILURE; + } -int main(int argc, char *argv[]){ - if (argc!=4){ - std::cerr << "Usage "<<argv[0]<<" input.off output.off N to load the file input.off, contract N edges and save the result to output.off.\n"; - return EXIT_FAILURE; - } - - Complex complex; + Complex complex; - // load the points - Skeleton_blocker_off_reader<Complex> off_reader(argv[1],complex); - if(!off_reader.is_valid()){ - std::cerr << "Unable to read file:"<<argv[1]<<std::endl; - return EXIT_FAILURE; - } + // load the points + Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file:" << argv[1] << std::endl; + return EXIT_FAILURE; + } - std::cout << "Load complex with "<<complex.num_vertices()<<" vertices"<<std::endl; + std::cout << "Load complex with " << complex.num_vertices() << " vertices" << std::endl; - int num_contractions = atoi(argv[3]); + int num_contractions = atoi(argv[3]); - boost::timer::auto_cpu_timer t; + boost::timer::auto_cpu_timer t; - // constructs the contractor object with Garland Heckbert policies. - Complex_contractor contractor(complex, - new GH_cost(complex), - new GH_placement(complex), - contraction::make_link_valid_contraction<EdgeProfile>(), - new GH_visitor(complex) - ); + // constructs the contractor object with Garland Heckbert policies. + Complex_contractor contractor(complex, + new GH_cost(complex), + new GH_placement(complex), + contraction::make_link_valid_contraction<EdgeProfile>(), + new GH_visitor(complex) + ); - std::cout<<"Contract "<<num_contractions<<" edges"<<std::endl; - contractor.contract_edges(num_contractions); + std::cout << "Contract " << num_contractions << " edges" << std::endl; + contractor.contract_edges(num_contractions); - std::cout << "Final complex has "<< - complex.num_vertices()<<" vertices, "<< - complex.num_edges()<<" edges and" << - complex.num_triangles()<<" triangles."<<std::endl; + std::cout << "Final complex has " << + complex.num_vertices() << " vertices, " << + complex.num_edges() << " edges and" << + complex.num_triangles() << " triangles." << std::endl; - //write simplified complex - Skeleton_blocker_off_writer<Complex> off_writer(argv[2],complex); + //write simplified complex + Skeleton_blocker_off_writer<Complex> off_writer(argv[2], complex); - return EXIT_SUCCESS; + return EXIT_SUCCESS; } diff --git a/src/Contraction/example/Garland_heckbert/Error_quadric.h b/src/Contraction/example/Garland_heckbert/Error_quadric.h index 725a3a56..725a3a56 100755..100644 --- a/src/Contraction/example/Garland_heckbert/Error_quadric.h +++ b/src/Contraction/example/Garland_heckbert/Error_quadric.h diff --git a/src/GudhUI/CMakeLists.txt b/src/GudhUI/CMakeLists.txt index b62be8ff..ddbae969 100644 --- a/src/GudhUI/CMakeLists.txt +++ b/src/GudhUI/CMakeLists.txt @@ -1,82 +1,69 @@ cmake_minimum_required(VERSION 2.8) project(GudhUI) -#Specify Gudhi's path here -#for instance -#set(Gudhi_Path "C:/Users/dsalinas/Documents/Recherche/Code/c++/Gudhi_library_1.1.0/include") -set(Gudhi_Path "/home/dsalinas/Documents/CodeSVN/gudhi/src/include") -if(MSVC) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4668 /wd4311 /wd4800 /wd4820 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018 -frounding-math") -else() - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -Wall -pedantic -frounding-math") -endif() - -set(EXECUTABLE_OUTPUT_PATH bin/${CMAKE_BUILD_TYPE}) - -find_package(Boost REQUIRED COMPONENTS) find_package(CGAL COMPONENTS Qt4) find_package(Qt4) find_package(QGLViewer) find_package(OpenGL) +message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}") +message("CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG}") +message("CMAKE_CXX_FLAGS_RELEASE ${CMAKE_CXX_FLAGS_RELEASE}") if ( CGAL_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND ) - set( QT_USE_QTXML TRUE ) - set( QT_USE_QTMAIN TRUE ) - set( QT_USE_QTSCRIPT TRUE ) - set( QT_USE_QTOPENGL TRUE ) - SET(Boost_USE_STATIC_LIBS ON) - SET(Boost_USE_MULTITHREAD OFF) - - INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS}) - LINK_DIRECTORIES(${Boost_LIBRARY_DIRS}) - - include(${QT_USE_FILE}) - include(${CGAL_USE_FILE}) - - include_directories (${QGLVIEWER_INCLUDE_DIR}) - include_directories(.) - include_directories(${Gudhi_Path}) + set( QT_USE_QTXML TRUE ) + set( QT_USE_QTMAIN TRUE ) + set( QT_USE_QTSCRIPT TRUE ) + set( QT_USE_QTOPENGL TRUE ) + SET(Boost_USE_STATIC_LIBS ON) + SET(Boost_USE_MULTITHREAD OFF) + INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS}) + LINK_DIRECTORIES(${Boost_LIBRARY_DIRS}) + + include(${QT_USE_FILE}) + include(${CGAL_USE_FILE}) + include_directories (${QGLVIEWER_INCLUDE_DIR}) + include_directories(.) - # qt : ui file, created wih Qt Designer ############### - set(CMAKE_CURRENT_BINARY_DIR "${CMAKE_CURRENT_SOURCE_DIR}/gui") - qt4_wrap_ui( uis - gui/main_window.ui - gui/MenuEdgeContraction.ui - gui/KNearestNeighborsMenu.ui - gui/UniformNeighborsMenu.ui - gui/PersistenceMenu.ui - ) + # qt : ui file, created wih Qt Designer ############### + set(CMAKE_CURRENT_BINARY_DIR "${CMAKE_CURRENT_SOURCE_DIR}/gui") + qt4_wrap_ui( uis + gui/main_window.ui + gui/MenuEdgeContraction.ui + gui/KNearestNeighborsMenu.ui + gui/UniformNeighborsMenu.ui + gui/PersistenceMenu.ui + ) - set(CMAKE_CURRENT_BINARY_DIR "${CMAKE_CURRENT_SOURCE_DIR}/gui") - qt4_automoc( - gui/MainWindow.cpp - gui/Menu_k_nearest_neighbors.cpp - gui/Menu_uniform_neighbors.cpp - gui/Menu_edge_contraction.cpp - gui/Menu_persistence.cpp - ) + set(CMAKE_CURRENT_BINARY_DIR "${CMAKE_CURRENT_SOURCE_DIR}/gui") + qt4_automoc( + gui/MainWindow.cpp + gui/Menu_k_nearest_neighbors.cpp + gui/Menu_uniform_neighbors.cpp + gui/Menu_edge_contraction.cpp + gui/Menu_persistence.cpp + ) - set(CMAKE_CURRENT_BINARY_DIR "${CMAKE_CURRENT_SOURCE_DIR}/view") - qt4_automoc(view/Viewer_instructor.cpp - view/Viewer.cpp - ) - ##################################################################### + set(CMAKE_CURRENT_BINARY_DIR "${CMAKE_CURRENT_SOURCE_DIR}/view") + qt4_automoc(view/Viewer_instructor.cpp + view/Viewer.cpp + ) + ##################################################################### - add_executable ( GudhUI - gui/gudhui.cpp - gui/MainWindow.cpp - gui/Menu_k_nearest_neighbors.cpp - gui/Menu_uniform_neighbors.cpp - gui/Menu_edge_contraction.cpp - gui/Menu_persistence.cpp - view/Viewer_instructor.cpp - view/Viewer.cpp - ${uis} - ) + add_executable ( GudhUI + gui/gudhui.cpp + gui/MainWindow.cpp + gui/Menu_k_nearest_neighbors.cpp + gui/Menu_uniform_neighbors.cpp + gui/Menu_edge_contraction.cpp + gui/Menu_persistence.cpp + view/Viewer_instructor.cpp + view/Viewer.cpp + ${uis} + ) - target_link_libraries( GudhUI ${QT_LIBRARIES} ${QGLVIEWER_LIBRARIES} ) - target_link_libraries( GudhUI ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ) + target_link_libraries( GudhUI ${QT_LIBRARIES} ${QGLVIEWER_LIBRARIES} ) + target_link_libraries( GudhUI ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ) else() message(STATUS "NOTICE: This demo requires CGAL, the QGLViewer, OpenGL and Qt4, and will not be compiled.") diff --git a/src/GudhUI/model/Model.h b/src/GudhUI/model/Model.h index 6ac971d0..17a7d278 100644 --- a/src/GudhUI/model/Model.h +++ b/src/GudhUI/model/Model.h @@ -72,8 +72,6 @@ class Model{ public: Complex complex_; typedef Complex::Vertex_handle Vertex_handle; - typedef Complex::CVI CVI; - Model():complex_(){ } @@ -317,7 +315,8 @@ private: void run_chomp(){ save_complex_in_file_for_chomp(); std::cout << "Call CHOMP library\n"; - system("../src/utils/homsimpl chomp.sim"); + int returnValue = system("utils/homsimpl chomp.sim"); + std::cout << "CHOMP returns" << returnValue << std::endl; } void save_complex_in_file_for_chomp(){ diff --git a/src/GudhUI/view/Viewer.h b/src/GudhUI/view/Viewer.h index 5639aa56..5639aa56 100755..100644 --- a/src/GudhUI/view/Viewer.h +++ b/src/GudhUI/view/Viewer.h diff --git a/src/GudhUI/view/Viewer_instructor.h b/src/GudhUI/view/Viewer_instructor.h index 9a2a236b..9a2a236b 100755..100644 --- a/src/GudhUI/view/Viewer_instructor.h +++ b/src/GudhUI/view/Viewer_instructor.h diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 3697a419..9487cce6 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -9,7 +9,6 @@ if (NOT MSVC) add_executable(rips_persistence rips_persistence.cpp) target_link_libraries(rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) - add_test(rips_persistence_2 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 2 -m 100) add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 3 -m 100) add_executable(persistence_from_file persistence_from_file.cpp) @@ -23,14 +22,14 @@ if (NOT MSVC) add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) target_link_libraries(rips_multifield_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) - add_test(rips_multifield_persistence_2_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 3 -m 100) add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 71 -m 100) add_executable ( performance_rips_persistence performance_rips_persistence.cpp ) target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) if(CGAL_FOUND) - add_definitions(-DDEBUG_TRACES) + # uncomment to display debug traces + # add_definitions(-DDEBUG_TRACES) add_executable(alpha_shapes_persistence alpha_shapes_persistence.cpp) target_link_libraries(alpha_shapes_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) add_test(alpha_shapes_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/alpha_shapes_persistence ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 2 0.5) diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index 2505897e..297a8f98 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -70,8 +70,6 @@ int main (int argc, char * argv[]) st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree st.expansion( dim_max ); // expand the graph until dimension dim_max - std::cout << "st=" << st << std::endl; - // Sort the simplices in the order of the filtration st.initialize_filtration(); diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index c42e4be4..b0d68f09 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -238,8 +238,8 @@ class Persistent_cohomology { transverse_idx_(), // key -> row persistent_pairs_(), interval_length_policy(&cpx, 0), - column_pool_(new boost::object_pool<Column>()), // memory pools for the CAM - cell_pool_(new boost::object_pool<Cell>()) { + column_pool_(), // memory pools for the CAM + cell_pool_() { Simplex_key idx_fil = 0; for (auto & sh : cpx_->filtration_simplex_range()) { cpx_->assign_key(sh, idx_fil); @@ -273,9 +273,6 @@ class Persistent_cohomology { transverse_ref.second.row_->clear(); delete transverse_ref.second.row_; } -// Clear the memory pools - delete column_pool_; - delete cell_pool_; } private: @@ -528,8 +525,8 @@ class Persistent_cohomology { Arith_element charac) { Simplex_key key = cpx_->key(sigma); // Create a column containing only one cell, - Column * new_col = column_pool_->construct(Column(key)); - Cell * new_cell = cell_pool_->construct(Cell(key, x, new_col)); + Column * new_col = column_pool_.construct(Column(key)); + Cell * new_cell = cell_pool_.construct(Cell(key, x, new_col)); new_col->col_.push_back(*new_cell); // and insert it in the matrix, in constant time thanks to the hint cam_.end(). // Indeed *new_col has the biggest lexicographic value because key is the @@ -585,7 +582,7 @@ class Persistent_cohomology { if (curr_col->col_.empty()) { // If the column is null ds_repr_[curr_col->class_key_] = NULL; - column_pool_->free(curr_col); // delete curr_col; + column_pool_.destroy(curr_col); // delete curr_col; } else { // Find whether the column obtained is already in the CAM result_insert_cam = cam_.insert(*curr_col); @@ -602,7 +599,7 @@ class Persistent_cohomology { Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_); ds_repr_[key_tmp] = &(*(result_insert_cam.first)); result_insert_cam.first->class_key_ = key_tmp; - column_pool_->free(curr_col); // delete curr_col; + column_pool_.destroy(curr_col); // delete curr_col; } } } else { @@ -634,7 +631,7 @@ class Persistent_cohomology { ++target_it; } else { if (target_it->key_ > other_it->first) { - Cell * cell_tmp = cell_pool_->construct(Cell(other_it->first // key + Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first // key , coeff_field_.additive_identity(), &target)); cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); @@ -652,8 +649,7 @@ class Persistent_cohomology { Cell * tmp_cell_ptr = &(*tmp_it); target.col_.erase(tmp_it); // removed from column - coeff_field_.clear_coefficient(tmp_cell_ptr->coefficient_); - cell_pool_->free(tmp_cell_ptr); // delete from memory + cell_pool_.destroy(tmp_cell_ptr); // delete from memory } else { ++target_it; ++other_it; @@ -662,7 +658,7 @@ class Persistent_cohomology { } } while (other_it != other.end()) { - Cell * cell_tmp = cell_pool_->construct(Cell(other_it->first, coeff_field_.additive_identity(), &target)); + Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first, coeff_field_.additive_identity(), &target)); cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); target.col_.insert(target.col_.end(), *cell_tmp); @@ -767,8 +763,8 @@ class Persistent_cohomology { std::list<Persistent_interval> persistent_pairs_; length_interval interval_length_policy; - boost::object_pool<Column> * column_pool_; - boost::object_pool<Cell> * cell_pool_; + boost::object_pool<Column> column_pool_; + boost::object_pool<Cell> cell_pool_; }; /** @} */ // end defgroup persistent_cohomology diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h index 5d0c9f9f..2a4c8692 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h @@ -76,9 +76,6 @@ class Field_Zp { return plus_times_equal(0, y, (Element)w); } - void clear_coefficient(Element x) { - } - Element plus_equal(const Element& x, const Element& y) { return plus_times_equal(x, y, (Element)1); } @@ -88,7 +85,7 @@ class Field_Zp { return add_id_all; } /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ - const Element& multiplicative_identity(Element P = 0) const { + const Element& multiplicative_identity(Element = 0) const { return mult_id_all; } /** Returns the inverse in the field. Modifies P.*/ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h index b43bb16e..c6fd5282 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h @@ -77,6 +77,7 @@ class Multi_field { mpz_nextprime(tmp_prime, tmp_prime); curr_prime = mpz_get_ui(tmp_prime); } + mpz_clear(tmp_prime); // set m to primorial(bound_prime) prod_characteristics_ = 1; for (auto p : primes_) { @@ -102,10 +103,6 @@ class Multi_field { } } - void clear_coefficient(Element & x) { - mpz_clear(x.get_mpz_t()); - } - /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ const Element& additive_identity() const { return add_id_all; diff --git a/src/Persistent_cohomology/test/CMakeLists.txt b/src/Persistent_cohomology/test/CMakeLists.txt index 9dc19251..ed63a6ac 100644 --- a/src/Persistent_cohomology/test/CMakeLists.txt +++ b/src/Persistent_cohomology/test/CMakeLists.txt @@ -1,29 +1,39 @@ cmake_minimum_required(VERSION 2.6) -project(GUDHITestSimplexTree) +project(GUDHIPersistentCohomologyUT) -if(NOT MSVC) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") endif() -add_executable ( persistent_cohomology_unit_test persistent_cohomology_unit_test.cpp ) -target_link_libraries(persistent_cohomology_unit_test ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +add_executable ( PersistentCohomologyUT persistent_cohomology_unit_test.cpp ) +target_link_libraries(PersistentCohomologyUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests -add_test(persistent_cohomology_unit_test ${CMAKE_CURRENT_BINARY_DIR}/persistent_cohomology_unit_test ${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt) +add_test(NAME PersistentCohomologyUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/PersistentCohomologyUT + ${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/PersistentCohomologyUT.xml --log_level=test_suite --report_level=no) if(GMPXX_FOUND AND GMP_FOUND) - add_executable ( persistent_cohomology_unit_test_multi_field persistent_cohomology_unit_test_multi_field.cpp ) - target_link_libraries(persistent_cohomology_unit_test_multi_field ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + add_executable ( PersistentCohomologyMultiFieldUT persistent_cohomology_unit_test_multi_field.cpp ) + target_link_libraries(PersistentCohomologyMultiFieldUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) - # Unitary tests - add_test(persistent_cohomology_unit_test_multi_field ${CMAKE_CURRENT_BINARY_DIR}/persistent_cohomology_unit_test_multi_field ${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt) -endif() + # Unitary tests + add_test(NAME PersistentCohomologyMultiFieldUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/PersistentCohomologyMultiFieldUT + ${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/PersistentCohomologyMultiFieldUT.xml --log_level=test_suite --report_level=no) -if (LCOV_PATH) - # Lcov code coverage of unitary test - add_test(src/Persistent_cohomology/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Persistent_cohomology) endif() -cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/include/gudhi") diff --git a/src/Persistent_cohomology/test/README b/src/Persistent_cohomology/test/README index ddceac63..6c64b5fe 100644 --- a/src/Persistent_cohomology/test/README +++ b/src/Persistent_cohomology/test/README @@ -1,7 +1,9 @@ To compile: *********** +cd /path-to-gudhi/ cmake . +cd /path-to-test/ make To launch with details: @@ -9,13 +11,13 @@ To launch with details: SINGLE FIELD ------------ -./persistent_cohomology_unit_test simplex_tree_file_for_unit_test.txt --report_level=detailed --log_level=all +./PersistentCohomologyUT simplex_tree_file_for_unit_test.txt --report_level=detailed --log_level=all ==> echo $? returns 0 in case of success (non-zero otherwise) MULTI FIELD ----------- -./persistent_cohomology_unit_test_multi_field simplex_tree_file_for_multi_field_unit_test.txt --report_level=detailed --log_level=all +./PersistentCohomologyMultiFieldUT simplex_tree_file_for_multi_field_unit_test.txt --report_level=detailed --log_level=all ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp index 1e7a74a7..55bc7066 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp @@ -1,4 +1,4 @@ -#define BOOST_TEST_MODULE const_string test +#define BOOST_TEST_MODULE persistent_cohomology test #include <boost/test/included/unit_test.hpp> #include <boost/system/error_code.hpp> #include <boost/chrono/thread_clock.hpp> @@ -22,8 +22,7 @@ using namespace boost::unit_test; typedef Simplex_tree<> typeST; std::string test_rips_persistence(int coefficient, int min_persistence) { - // Check file name is given as parameter from CMakeLists.txt - BOOST_CHECK(framework::master_test_suite().argc >= 2); + // file name is given as parameter from CMakeLists.txt const std::string inputFile(framework::master_test_suite().argv[1]); std::ifstream simplex_tree_stream; diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp index e88add3a..18a4725e 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp @@ -1,4 +1,4 @@ -#define BOOST_TEST_MODULE const_string test +#define BOOST_TEST_MODULE persistent_cohomology_multi_field test #include <boost/test/included/unit_test.hpp> #include <boost/system/error_code.hpp> #include <boost/chrono/thread_clock.hpp> @@ -23,8 +23,7 @@ using namespace boost::unit_test; typedef Simplex_tree<> typeST; std::string test_rips_persistence(int min_coefficient, int max_coefficient, int min_persistence) { - // Check file name is given as parameter from CMakeLists.txt - BOOST_CHECK(framework::master_test_suite().argc >= 2); + // file name is given as parameter from CMakeLists.txt const std::string inputFile(framework::master_test_suite().argv[1]); std::ifstream simplex_tree_stream; @@ -38,9 +37,9 @@ std::string test_rips_persistence(int min_coefficient, int max_coefficient, int << " - filtration= " << st.filtration() << std::endl; // Check - BOOST_CHECK(st.num_simplices() == 6142604); + BOOST_CHECK(st.num_simplices() == 58); BOOST_CHECK(st.dimension() == 3); - BOOST_CHECK(st.filtration() == 0.249999); + BOOST_CHECK(st.filtration() == 0.4); // Sort the simplices in the order of the filtration st.initialize_filtration(); @@ -48,9 +47,9 @@ std::string test_rips_persistence(int min_coefficient, int max_coefficient, int // Compute the persistence diagram of the complex Persistent_cohomology<Simplex_tree<>, Multi_field> pcoh(st); - pcoh.init_coefficients(min_coefficient, max_coefficient); // initializes the coefficient field for homology + pcoh.init_coefficients(min_coefficient, max_coefficient); // initializes the coefficient field for homology // Check infinite rips - pcoh.compute_persistent_cohomology(min_persistence); // Minimal lifetime of homology feature to be recorded. + pcoh.compute_persistent_cohomology(min_persistence); // Minimal lifetime of homology feature to be recorded. std::ostringstream ossRips; pcoh.output_diagram(ossRips); @@ -60,68 +59,54 @@ std::string test_rips_persistence(int min_coefficient, int max_coefficient, int } void test_rips_persistence_in_dimension(int min_dimension, int max_dimension) { - std::string value0(" 0 0 inf"); - std::string value1(" 1 0.0702103 inf"); - std::string value2("2 1 0.0702103 inf"); - std::string value3("2 2 0.159992 inf"); + // there are 2 discontinued ensembles + std::string value0(" 0 0.25 inf"); + std::string value1(" 1 0.4 inf"); + // And a big hole - cut in 2 pieces after 0.3 + std::string value2(" 0 0.2 0.3"); + + // For dim <= 1 => + std::string value3(" 1 0.25 inf"); + std::string value4(" 2 0.25 inf"); + std::string value5(" 1 0.3 inf"); + std::string value6(" 2 0.3 inf"); + std::string value7(" 2 0.4 inf"); std::cout << "********************************************************************" << std::endl; std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD MIN_DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=0" << std::endl; - std::string str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 1); - - BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found - std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; - - std::cout << "********************************************************************" << std::endl; - std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=2" << std::endl; - - str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 2); - - BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + std::string str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, static_cast<Filtration_value> (0.0)); std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; - std::cout << "********************************************************************" << std::endl; - std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=3" << std::endl; - - str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 3); - BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found - std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; - std::cout << "********************************************************************" << std::endl; - std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=Inf" << std::endl; + if ((min_dimension < 2) && (max_dimension < 2)) { + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value4) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value5) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value6) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value7) != std::string::npos); // Check found + } else { + BOOST_CHECK(str_rips_persistence.find(value3) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value4) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value5) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value6) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value7) == std::string::npos); // Check not found + } - str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, std::numeric_limits<int>::max()); - - BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found - BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found - std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; } -BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_multi_field_dim_1_2 ) -{ - test_rips_persistence_in_dimension(1, 2); +BOOST_AUTO_TEST_CASE(rips_persistent_cohomology_multi_field_dim_1_2) { + test_rips_persistence_in_dimension(0, 1); } -BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_multi_field_dim_2_3 ) -{ - test_rips_persistence_in_dimension(2, 3); +BOOST_AUTO_TEST_CASE(rips_persistent_cohomology_multi_field_dim_2_3) { + test_rips_persistence_in_dimension(1, 3); } -BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_multi_field_dim_1_5 ) -{ +BOOST_AUTO_TEST_CASE(rips_persistent_cohomology_multi_field_dim_1_5) { test_rips_persistence_in_dimension(1, 5); } diff --git a/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt new file mode 100644 index 00000000..ed2c0c3d --- /dev/null +++ b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt @@ -0,0 +1,58 @@ +0 0 0.2 +0 3 0.2 +1 3 0 0.2 +0 6 0.2 +0 11 0.2 +1 11 6 0.2 +0 13 0.25 +0 14 0.25 +1 14 13 0.25 +0 15 0.25 +1 15 13 0.25 +1 15 14 0.25 +2 15 14 13 0.25 +0 1 0.3 +1 1 0 0.3 +0 2 0.3 +1 2 0 0.3 +1 2 1 0.3 +2 2 1 0 0.3 +0 4 0.3 +1 4 3 0.3 +0 5 0.3 +1 5 3 0.3 +1 5 4 0.3 +2 5 4 3 0.3 +0 9 0.3 +0 10 0.3 +1 10 2 0.3 +1 10 9 0.3 +1 11 9 0.3 +1 11 10 0.3 +2 11 10 9 0.3 +0 12 0.3 +1 12 2 0.3 +1 12 10 0.3 +2 12 10 2 0.3 +1 6 0 0.4 +1 6 1 0.4 +2 6 1 0 0.4 +0 7 0.4 +1 7 0 0.4 +1 7 1 0.4 +2 7 1 0 0.4 +1 7 6 0.4 +2 7 6 0 0.4 +2 7 6 1 0.4 +3 7 6 1 0 0.4 +0 8 0.4 +1 8 4 0.4 +1 8 5 0.4 +2 8 5 4 0.4 +1 9 4 0.4 +1 9 5 0.4 +2 9 5 4 0.4 +1 9 8 0.4 +2 9 8 4 0.4 +2 9 8 5 0.4 +3 9 8 5 4 0.4 diff --git a/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt.REMOVED.git-id b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt.REMOVED.git-id deleted file mode 100644 index 2dd38515..00000000 --- a/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt.REMOVED.git-id +++ /dev/null @@ -1 +0,0 @@ -ce87199d425b05f51c74cbf635870bfa5abbc7a1
\ No newline at end of file diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp index bde224f1..6d20e43e 100644 --- a/src/Simplex_tree/example/simple_simplex_tree.cpp +++ b/src/Simplex_tree/example/simple_simplex_tree.cpp @@ -28,7 +28,6 @@ using namespace Gudhi; typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; int main(int argc, char * const argv[]) { @@ -58,9 +57,8 @@ int main(int argc, char * const argv[]) { std::cout << " * INSERT 0" << std::endl; typeVectorVertex firstSimplexVector; firstSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); typePairSimplexBool returnValue = - simplexTree.insert_simplex(firstSimplex.first, firstSimplex.second); + simplexTree.insert_simplex(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 0 INSERTED" << std::endl; @@ -74,9 +72,8 @@ int main(int argc, char * const argv[]) { std::cout << " * INSERT 1" << std::endl; typeVectorVertex secondSimplexVector; secondSimplexVector.push_back(SECOND_VERTEX_HANDLE); - typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(secondSimplex.first, secondSimplex.second); + simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 1 INSERTED" << std::endl; @@ -91,9 +88,8 @@ int main(int argc, char * const argv[]) { typeVectorVertex thirdSimplexVector; thirdSimplexVector.push_back(FIRST_VERTEX_HANDLE); thirdSimplexVector.push_back(SECOND_VERTEX_HANDLE); - typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(thirdSimplex.first, thirdSimplex.second); + simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (0,1) INSERTED" << std::endl; @@ -107,9 +103,8 @@ int main(int argc, char * const argv[]) { std::cout << " * INSERT 2" << std::endl; typeVectorVertex fourthSimplexVector; fourthSimplexVector.push_back(THIRD_VERTEX_HANDLE); - typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(fourthSimplex.first, fourthSimplex.second); + simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 2 INSERTED" << std::endl; @@ -124,9 +119,8 @@ int main(int argc, char * const argv[]) { typeVectorVertex fifthSimplexVector; fifthSimplexVector.push_back(THIRD_VERTEX_HANDLE); fifthSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(fifthSimplex.first, fifthSimplex.second); + simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,0) INSERTED" << std::endl; @@ -141,9 +135,8 @@ int main(int argc, char * const argv[]) { typeVectorVertex sixthSimplexVector; sixthSimplexVector.push_back(THIRD_VERTEX_HANDLE); sixthSimplexVector.push_back(SECOND_VERTEX_HANDLE); - typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(sixthSimplex.first, sixthSimplex.second); + simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,1) INSERTED" << std::endl; @@ -159,9 +152,8 @@ int main(int argc, char * const argv[]) { seventhSimplexVector.push_back(THIRD_VERTEX_HANDLE); seventhSimplexVector.push_back(SECOND_VERTEX_HANDLE); seventhSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(seventhSimplex.first, seventhSimplex.second); + simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,1,0) INSERTED" << std::endl; @@ -175,9 +167,8 @@ int main(int argc, char * const argv[]) { std::cout << " * INSERT 3" << std::endl; typeVectorVertex eighthSimplexVector; eighthSimplexVector.push_back(FOURTH_VERTEX_HANDLE); - typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(eighthSimplex.first, eighthSimplex.second); + simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + 3 INSERTED" << std::endl; @@ -192,9 +183,8 @@ int main(int argc, char * const argv[]) { typeVectorVertex ninethSimplexVector; ninethSimplexVector.push_back(FOURTH_VERTEX_HANDLE); ninethSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(ninethSimplex.first, ninethSimplex.second); + simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (3,0) INSERTED" << std::endl; @@ -208,9 +198,8 @@ int main(int argc, char * const argv[]) { std::cout << " * INSERT 0 (already inserted)" << std::endl; typeVectorVertex tenthSimplexVector; tenthSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); // With a different filtration value returnValue = - simplexTree.insert_simplex(tenthSimplex.first, tenthSimplex.second); + simplexTree.insert_simplex(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); // With a different filtration value if (returnValue.second == true) { std::cout << " + 0 INSERTED" << std::endl; @@ -226,9 +215,8 @@ int main(int argc, char * const argv[]) { eleventhSimplexVector.push_back(THIRD_VERTEX_HANDLE); eleventhSimplexVector.push_back(SECOND_VERTEX_HANDLE); eleventhSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); returnValue = - simplexTree.insert_simplex(eleventhSimplex.first, eleventhSimplex.second); + simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); if (returnValue.second == true) { std::cout << " + (2,1,0) INSERTED" << std::endl; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index f95679cb..5d3753ca 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -326,7 +326,7 @@ class Simplex_tree { /** \brief Returns the filtration value of a simplex. * * Called on the null_simplex, returns INFINITY. */ - Filtration_value filtration(Simplex_handle sh) { + Filtration_value filtration(Simplex_handle sh) const { if (sh != null_simplex()) { return sh->second.filtration(); } else { @@ -334,34 +334,34 @@ class Simplex_tree { } // filtration(); } } /** \brief Returns an upper bound of the filtration values of the simplices. */ - Filtration_value filtration() { + Filtration_value filtration() const { return threshold_; } /** \brief Returns a Simplex_handle different from all Simplex_handles * associated to the simplices in the simplicial complex. * * One can call filtration(null_simplex()). */ - Simplex_handle null_simplex() { + Simplex_handle null_simplex() const { return Dictionary_it(NULL); } /** \brief Returns a key different for all keys associated to the * simplices of the simplicial complex. */ - Simplex_key null_key() { + Simplex_key null_key() const { return -1; } /** \brief Returns a Vertex_handle different from all Vertex_handles associated * to the vertices of the simplicial complex. */ - Vertex_handle null_vertex() { + Vertex_handle null_vertex() const { return null_vertex_; } /** \brief Returns the number of vertices in the complex. */ - size_t num_vertices() { + size_t num_vertices() const { return root_.members_.size(); } /** \brief Returns the number of simplices in the complex. * * Does not count the empty simplex. */ - const unsigned int& num_simplices() const { + unsigned int num_simplices() const { return num_simplices_; } @@ -378,13 +378,13 @@ class Simplex_tree { return dim - 1; } /** \brief Returns an upper bound on the dimension of the simplicial complex. */ - int dimension() { + int dimension() const { return dimension_; } /** \brief Returns true iff the node in the simplex tree pointed by * sh has children.*/ - bool has_children(Simplex_handle sh) { + bool has_children(Simplex_handle sh) const { return (sh->second.children()->parent() == sh->first); } @@ -564,7 +564,7 @@ class Simplex_tree { threshold_ = fil; } /** Set a number of simplices for the simplicial complex. */ - void set_num_simplices(const unsigned int& num_simplices) { + void set_num_simplices(unsigned int num_simplices) { num_simplices_ = num_simplices; } /** Set a dimension for the simplicial complex. */ @@ -772,10 +772,10 @@ class Simplex_tree { } /** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2) * and assigns the maximal possible Filtration_value to the Nodes. */ - void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection, - Dictionary_it begin1, Dictionary_it end1, - Dictionary_it begin2, Dictionary_it end2, - Filtration_value filtration) { + static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection, + Dictionary_it begin1, Dictionary_it end1, + Dictionary_it begin2, Dictionary_it end2, + Filtration_value filtration) { if (begin1 == end1 || begin2 == end2) return; // ----->> while (true) { @@ -802,8 +802,8 @@ class Simplex_tree { } } /** Maximum over 3 values.*/ - Filtration_value maximum(Filtration_value a, Filtration_value b, - Filtration_value c) { + static Filtration_value maximum(Filtration_value a, Filtration_value b, + Filtration_value c) { Filtration_value max = (a < b) ? b : a; return ((max < c) ? c : max); } diff --git a/src/Simplex_tree/test/CMakeLists.txt b/src/Simplex_tree/test/CMakeLists.txt index 02ef9d8b..b6a1c0b6 100644 --- a/src/Simplex_tree/test/CMakeLists.txt +++ b/src/Simplex_tree/test/CMakeLists.txt @@ -1,21 +1,25 @@ cmake_minimum_required(VERSION 2.6) -project(GUDHITestSimplexTree) +project(GUDHISimplexTreeUT) -if(NOT MSVC) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") endif() -add_executable ( simplex_tree_unit_test simplex_tree_unit_test.cpp ) -target_link_libraries(simplex_tree_unit_test ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +add_executable ( SimplexTreeUT simplex_tree_unit_test.cpp ) +target_link_libraries(SimplexTreeUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests -add_test(simplex_tree_unit_test ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_unit_test) - -if (LCOV_PATH) - # Lcov code coverage of unitary test - add_test(src/Simplex_tree/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Simplex_tree) -endif() +add_test(NAME SimplexTreeUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/SimplexTreeUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/SimplexTreeUT.xml --log_level=test_suite --report_level=no) -cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Simplex_tree/include/gudhi") diff --git a/src/Simplex_tree/test/README b/src/Simplex_tree/test/README index 620bcd5f..21c3d871 100644 --- a/src/Simplex_tree/test/README +++ b/src/Simplex_tree/test/README @@ -1,12 +1,14 @@ To compile: *********** +cd /path-to-gudhi/ cmake . +cd /path-to-test/ make To launch with details: *********************** -./simplex_tree_unit_test --report_level=detailed --log_level=all +./SimplexTreeUT --report_level=detailed --log_level=all ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index c0cfced1..6b0a1f3d 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -1,4 +1,4 @@ -#define BOOST_TEST_MODULE const_string test +#define BOOST_TEST_MODULE simplex_tree test #include <boost/test/included/unit_test.hpp> #include <boost/system/error_code.hpp> #include <boost/chrono/thread_clock.hpp> diff --git a/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp index 92fa17f3..126e32ec 100644 --- a/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp +++ b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp @@ -64,8 +64,10 @@ int main (int argc, char *argv[]){ // or edges, complex.num_vertices() and complex.num_edges() are // more appropriated! unsigned num_vertices = 0; - for(auto v : complex.vertex_range()) - ++num_vertices; + for(auto v : complex.vertex_range()) { + std::cout << "Vertex " << v <<std::endl; + ++num_vertices; + } // such loop can also be done directly with distance as iterators are STL compliant auto edges = complex.edge_range(); diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h index 049db6d5..289819b5 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h @@ -109,7 +109,7 @@ and point access in addition. \subsection Visitor The class Skeleton_blocker_complex has a visitor that is called when usual operations such as adding an edge or remove a vertex are called. -You may want to use this visitor to compute statistics or to update another data-structure (for instance this visitor is heavily used in the \ref contr package. +You may want to use this visitor to compute statistics or to update another data-structure (for instance this visitor is heavily used in the \ref contr package). diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h index 86a12d90..dd8d898e 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -223,7 +223,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex_hand for (auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it) for (auto v_it = u_it; ++v_it != sigma.end(); /**/) { - std::cout << "add edge" << *u_it << " " << *v_it << std::endl; + // std::cout << "add edge" << *u_it << " " << *v_it << std::endl; add_edge(*u_it, *v_it); } remove_blocker_include_in_simplex(sigma); diff --git a/src/Skeleton_blocker/test/CMakeLists.txt b/src/Skeleton_blocker/test/CMakeLists.txt index c69bfec7..e62600a2 100644 --- a/src/Skeleton_blocker/test/CMakeLists.txt +++ b/src/Skeleton_blocker/test/CMakeLists.txt @@ -1,12 +1,17 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIskbl) -if(NOT MSVC) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g") - - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") endif() add_executable(TestSkeletonBlockerComplex TestSkeletonBlockerComplex.cpp) @@ -17,10 +22,3 @@ add_test(TestSkeletonBlockerComplex ${CMAKE_CURRENT_BINARY_DIR}/TestSkeletonBloc add_test(TestSimplifiable ${CMAKE_CURRENT_BINARY_DIR}/TestSimplifiable) add_test(TestGeometricComplex ${CMAKE_CURRENT_BINARY_DIR}/TestGeometricComplex) - -if (LCOV_PATH) - # Lcov code coverage of unitary test - add_test(src/Skeleton_blocker/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Skeleton_blocker) -endif() - -cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Skeleton_blocker/include/gudhi") diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 23919b4a..ff372d16 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -1,16 +1,7 @@ cmake_minimum_required(VERSION 2.6) project(GUDHIWitnessComplex) -#cmake -DCGAL_DIR=~/GitDrive/CGAL/ ../../.. -#if (CGAL_FOUND) - #message(STATUS "CGAL version: ${CGAL_VERSION}.") - #include( ${CGAL_USE_FILE} ) - - #find_package(Eigen3 3.1.0) - #include( ${EIGEN3_USE_FILE} ) - - #INCLUDE_DIRECTORIES(${EIGEN3_INCLUDE_DIR}) - #INCLUDE_DIRECTORIES(${CGAL_INCLUDE_DIRS}) +# A simple example add_executable ( simple_witness_complex simple_witness_complex.cpp ) add_test(simple_witness_complex ${CMAKE_CURRENT_BINARY_DIR}/simple_witness_complex) @@ -21,7 +12,7 @@ project(GUDHIWitnessComplex) add_executable( witness_complex_from_off witness_complex_from_off.cpp ) add_executable( witness_complex_from_wl_matrix witness_complex_from_wl_matrix.cpp ) -#endif() + # An example with Simplex-tree using CGAL alpha_shapes_3 @@ -51,6 +42,7 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") include( ${EIGEN3_USE_FILE} ) + message(STATUS "Eigen3 use file: ${EIGEN3_USE_FILE}.") include_directories (BEFORE "../../include") add_executable ( witness_complex_knn_landmarks witness_complex_knn_landmarks.cpp ) diff --git a/src/Witness_complex/example/protected_sets/output_tikz.h b/src/Witness_complex/example/protected_sets/output_tikz.h new file mode 100644 index 00000000..edfd9a5f --- /dev/null +++ b/src/Witness_complex/example/protected_sets/output_tikz.h @@ -0,0 +1,67 @@ +#ifndef OUTPUT_TIKZ_H +#define OUTPUT_TIKZ_H + +#include <vector> +#include <string> +#include <algorithm> +#include <fstream> +#include <cmath> + +void write_tikz_plot(std::vector<FT> data, std::string filename) +{ + int n = data.size(); + FT vmax = *(std::max_element(data.begin(), data.end())); + //std::cout << std::log10(vmax) << " " << std::floor(std::log10(vmax)); + + FT order10 = pow(10,std::floor(std::log10(vmax))); + int digit = std::floor( vmax / order10) + 1; + if (digit == 4 || digit == 6) digit = 5; + if (digit > 6) digit = 10; + FT plot_max = digit*order10; + std::cout << plot_max << " " << vmax; + FT hstep = 10.0/(n-1); + FT wstep = 10.0 / plot_max; + + std::cout << "(eps_max-eps_min)/(N-48) = " << (vmax-*data.begin())/(data.size()-48) << "\n"; + std::ofstream ofs(filename, std::ofstream::out); + + ofs << + "\\documentclass{standalone}\n" << + "\\usepackage[utf8]{inputenc}\n" << + "\\usepackage{amsmath}\n" << + "\\usepackage{tikz}\n\n" << + "\\begin{document}\n" << + "\\begin{tikzpicture}\n"; + + ofs << "\\draw[->] (0,0) -- (0,11);" << std::endl << + "\\draw[->] (0,0) -- (11,0);" << std::endl << + "\\foreach \\i in {1,...,10}" << std::endl << + "\\draw (0,\\i) -- (-0.05,\\i);" << std::endl << + "\\foreach \\i in {1,...,10}" << std::endl << + "\\draw (\\i,0) -- (\\i,-0.05);" << std::endl << std::endl << + + "\\foreach \\i in {1,...,10}" << std::endl << + "\\draw[dashed] (-0.05,\\i) -- (11,\\i);" << std::endl << std::endl << + + "\\node at (-0.5,11) {$*$}; " << std::endl << + "\\node at (11,-0.5) {$*$}; " << std::endl << + "\\node at (-0.5,-0.5) {0}; " << std::endl << + "\\node at (-0.5,10) {" << plot_max << "}; " << std::endl << + "%\\node at (10,-0.5) {2}; " << std::endl; + + ofs << "\\draw[red] (0," << wstep*data[0] << ")"; + for (int i = 1; i < n; ++i) + ofs << " -- (" << hstep*i << "," << wstep*data[i] << ")"; + ofs << ";\n"; + + ofs << + "\\end{tikzpicture}\n" << + "\\end{document}"; + + ofs.close(); + + + +} + +#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets.h b/src/Witness_complex/example/protected_sets/protected_sets.h new file mode 100644 index 00000000..ec627808 --- /dev/null +++ b/src/Witness_complex/example/protected_sets/protected_sets.h @@ -0,0 +1,597 @@ +#ifndef PROTECTED_SETS_H +#define PROTECTED_SETS_H + +#include <algorithm> +#include <CGAL/Cartesian_d.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Euclidean_distance.h> +#include <CGAL/Kernel_d/Sphere_d.h> +#include <CGAL/Kernel_d/Hyperplane_d.h> +#include <CGAL/Kernel_d/Vector_d.h> + +#include <CGAL/point_generators_d.h> +#include <CGAL/constructions_d.h> + + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::Point_d Point_d; +typedef K::Vector_d Vector_d; +typedef K::Oriented_side_d Oriented_side_d; +typedef K::Has_on_positive_side_d Has_on_positive_side_d; +typedef K::Sphere_d Sphere_d; +typedef K::Hyperplane_d Hyperplane_d; + +typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; +typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; + +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; + +FT _sfty = pow(10,-14); + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// +// AUXILLARY FUNCTIONS +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well + * + * W is the initial point vector + * chosen_landmark is the index of the chosen point in W + * landmarks_ind is the vector of indices of already chosen points in W + * delaunay is the Delaunay triangulation + * landmark_count is the current number of chosen vertices + * torus is true iff you are working on a flat torus [-1,1]^d + * OUT: Vertex handle to the newly inserted point + */ +Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) +{ + if (!torus) + { + Delaunay_vertex v =delaunay.insert(W[chosen_landmark]); + landmarks_ind.push_back(chosen_landmark); + landmark_count++; + return v; + } + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + Delaunay_vertex v; + for (int i = 0; i < nb_cells; ++i) + { + std::vector<FT> point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + v = delaunay.insert(point); + } + landmarks_ind.push_back(chosen_landmark); + landmark_count++; + return v; + } +} + +/** Small check if the vertex v is in the full cell fc + */ + +bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) +{ + for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) + if (*v_it == v) + return true; + return false; +} + +/** Fill chosen point vector from indices with copies if you are working on a flat torus + * + * IN: W is the point vector + * OUT: landmarks is the output vector + * IN: landmarks_ind is the vector of indices + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) +{ + if (!torus) + for (unsigned j = 0; j < landmarks_ind.size(); ++j) + landmarks.push_back(W[landmarks_ind[j]]); + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + int nbL = landmarks_ind.size(); + // Fill landmarks + for (int i = 0; i < nb_cells-1; ++i) + for (int j = 0; j < nbL; ++j) + { + int cell_i = i; + Point_d point; + for (int l = 0; l < D; ++l) + { + point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1)); + cell_i /= 3; + } + landmarks.push_back(point); + } + } +} + +/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices + * + * IN: t is the Delaunay triangulation + * OUT: full_cells is the output vector + */ + +void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) +{ + // Store vertex indices in a map + int ind = 0; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (t.is_infinite(v_it)) + continue; + else + index_of_vertex[v_it] = ind++; + // Write full cells as vectors in full_cells + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + std::vector<int> cell; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + cell.push_back(index_of_vertex[*v_it]); + full_cells.push_back(cell); + } +} + +//////////////////////////////////////////////////////////////////////////////////////////////////////////// +// IS VIOLATED TEST +//////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Check if a newly created cell is protected from old vertices + * + * t is the Delaunay triangulation + * vertices is the vector containing the point to insert and a facet f in t + * v1 is the vertex of t, such that f and v1 form a simplex + * v2 is the vertex of t, such that f and v2 form another simplex + * delta is the protection constant + * power_protection is true iff the delta-power protection is used + */ + +bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) +{ + assert(vertices.size() == vertices[0].size() || + vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 + assert(v1 != v2); + if (vertices.size() == vertices[0].size() + 1) + // FINITE CASE + { + Sphere_d cs(vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + { + //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + */ + // Check if the simplex is thick enough + Hyperplane_d tau_h(vertices.begin()+1, vertices.end()); + Vector_d orth_tau = tau_h.orthogonal_vector(); + /* + p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1)); + */ + //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n"; + FT orth_length = sqrt(orth_tau.squared_length()); + K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; + // Compute the altitude + FT h = 0; + for (o_it = orth_tau.cartesian_begin(), + p_it = vertices.begin()->cartesian_begin(), + s_it = (vertices.begin()+1)->cartesian_begin(); + o_it != orth_tau.cartesian_end(); + ++o_it, ++p_it, ++s_it) + h += (*o_it)*(*p_it - *s_it)/orth_length; + h = fabs(h); + // Is the center inside the box? + bool inside_the_box = true; + for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) + if (*c_it > 1.0 || *c_it < -1.0) + { + inside_the_box = false; break; + } + if (inside_the_box && h/r < theta0) + return true; + if (!t.is_infinite(v1)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + if (!t.is_infinite(v2)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + else + // INFINITE CASE + { + Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); + while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) + v++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); + Vector_d orth_v = facet_plane.orthogonal_vector(); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + std::vector<FT> coords; + Point_d p = v_it->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!p_is_inside && p_delta_is_inside) + return true; + } + */ + if (!t.is_infinite(v1)) + { + std::vector<FT> coords; + Point_d p = v1->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + if (!t.is_infinite(v2)) + { + std::vector<FT> coords; + Point_d p = v2->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + } + return false; +} + +/** Auxillary recursive function to check if the point p violates the protection of the cell c and + * if there is a violation of an eventual new cell + * + * p is the point to insert + * t is the current triangulation + * c is the current cell (simplex) + * parent_cell is the parent cell (simplex) + * index is the index of the facet between c and parent_cell from parent_cell's point of view + * D is the dimension of the triangulation + * delta is the protection constant + * marked_cells is the vector of all visited cells containing p in their circumscribed ball + * power_protection is true iff you are working with delta-power protection + * + * OUT: true iff inserting p hasn't produced any violation so far + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0) +{ + Euclidean_distance ed; + std::vector<Point_d> vertices; + if (!t.is_infinite(c)) + { + // if the cell is finite, we look if the protection is violated + for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); + FT dist2 = ed.transformed_distance(center_cs, p); + // if the new point is inside the protection ball of a non conflicting simplex + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + // if the new point is inside the circumscribing ball : continue violation searching on neighbours + //if (dist2 < r*r) + //if (dist2 < (5*r+delta)*(5*r+delta)) + if (dist2 < r*r) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) + return true; + } + } + // if the new point is outside the protection sphere + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is guaranteed to be finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + else + { + // Inside of the convex hull is + side. Outside is - side. + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!t.is_infinite(*vh_it)) + vertices.push_back((*vh_it)->point()); + Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); + while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) + v_it++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); + //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); + Vector_d orth_v = facet_plane.orthogonal_vector(); + std::vector<FT> coords; + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + + // If we work with power protection, we just ignore any conflicts + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + //if the cell is infinite we look at the neighbours regardless + if (p_is_inside) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) + return true; + } + } + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is finite if the parent cell is finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + if (!t.is_infinite(parent_cell->vertex(i))) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + //c->tds_data().clear_visited(); + //marked_cells.pop_back(); + return false; +} + +/** Checks if inserting the point p in t will make conflicts + * + * p is the point to insert + * t is the current triangulation + * D is the dimension of triangulation + * delta is the protection constant + * power_protection is true iff you are working with delta-power protection + * OUT: true iff inserting p produces a violation of delta-protection. + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + std::vector<Full_cell_handle> marked_cells; + c = t.locate(p, lt, f, ft, v); + bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0); + for (Full_cell_handle fc : marked_cells) + fc->tds_data().clear(); + return violation_existing_cells; +} + +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// +//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////// +// LANDMARK CHOICE PROCEDURE +/////////////////////////////////////////////////////////////////////// + +/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. + * + * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d + * IN: nbP is the size of W + * OUT: landmarks is the output vector for the points + * OUT: landmarks_ind is the output vector for the indices of the selected points in W + * IN: delta is the constant of protection + * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT delta, std::vector<std::vector<int>>& full_cells, bool torus, bool power_protection, FT theta0) +{ + bool return_ = true; + unsigned D = W[0].size(); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + CGAL::Random rand; + int landmark_count = 0; + std::list<int> index_list; + // shuffle the list of indexes (via a vector) + { + std::vector<int> temp_vector; + for (int i = 0; i < nbP; ++i) + temp_vector.push_back(i); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); + for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + index_list.push_front(*it); + } + if (!torus) + for (unsigned pos1 = 0; pos1 < D+1; ++pos1) + { + std::vector<FT> point; + for (unsigned i = 0; i < pos1; ++i) + point.push_back(-1); + if (pos1 != D) + point.push_back(1); + for (unsigned i = pos1+1; i < D; ++i) + point.push_back(0); + assert(point.size() == D); + W[index_list.front()] = Point_d(point); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + } + else if (D == 2) + { + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 2; ++j) + { + W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0}); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j}); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + } + } + else + std::cout << "No torus starter available for dim>2\n"; + std::list<int>::iterator list_it = index_list.begin(); + while (list_it != index_list.end()) + { + if (!is_violating_protection(W[*list_it], t, D, delta, power_protection, theta0)) + { + // If no conflicts then insert in every copy of T^3 + + insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus); + if (return_) + { + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + index_list.erase(list_it++); + /* + // PIECE OF CODE FOR DEBUGGING PURPOSES + + Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + if (triangulation_is_protected(t, delta)) + { + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + { //THAT'S WHERE SOMETHING'S WRONG + t.remove(inserted_v); + landmarks_ind.pop_back(); + landmark_count--; + write_delaunay_mesh(t, W[*list_it], is2d); + is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore + } + */ + //std::cout << "index_list_size() = " << index_list.size() << "\n"; + } + else + { + list_it++; + //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; + } + //if (list_it != index_list.end()) + // write_delaunay_mesh(t, W[*list_it], is2d); + } + fill_landmarks(W, landmarks, landmarks_ind, torus); + fill_full_cell_vector(t, full_cells); + /* + if (triangulation_is_protected(t, delta)) + std::cout << "Triangulation is ok\n"; + else + { + std::cout << "Triangulation is BAD!! T_T ã—ãã—ã!\n"; + } + */ + //write_delaunay_mesh(t, W[0], is2d); + //std::cout << t << std::endl; +} + +#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp b/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp new file mode 100644 index 00000000..f3df3f1e --- /dev/null +++ b/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp @@ -0,0 +1,610 @@ +#ifndef PROTECTED_SETS_H +#define PROTECTED_SETS_H + +#include <algorithm> +#include <CGAL/Cartesian_d.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Euclidean_distance.h> +#include <CGAL/Kernel_d/Sphere_d.h> +#include <CGAL/Kernel_d/Hyperplane_d.h> +#include <CGAL/Kernel_d/Vector_d.h> + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::Point_d Point_d; +typedef K::Vector_d Vector_d; +typedef K::Oriented_side_d Oriented_side_d; +typedef K::Has_on_positive_side_d Has_on_positive_side_d; +typedef K::Sphere_d Sphere_d; +typedef K::Hyperplane_d Hyperplane_d; + +typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; +typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; + +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; + +FT _sfty = pow(10,-14); + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// +// AUXILLARY FUNCTIONS +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well + * + * W is the initial point vector + * chosen_landmark is the index of the chosen point in W + * landmarks_ind is the vector of indices of already chosen points in W + * delaunay is the Delaunay triangulation + * landmark_count is the current number of chosen vertices + * torus is true iff you are working on a flat torus [-1,1]^d + * OUT: Vertex handle to the newly inserted point + */ +Delaunay_vertex insert_delaunay_landmark_with_copies(Point_d& p, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) +{ + if (!torus) + { + Delaunay_vertex v =delaunay.insert(p); + landmark_count++; + return v; + } + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + Delaunay_vertex v; + for (int i = 0; i < nb_cells; ++i) + { + std::vector<FT> point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(p[l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + v = delaunay.insert(point); + } + landmark_count++; + return v; + } +} + +/** Small check if the vertex v is in the full cell fc + */ + +bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) +{ + for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) + if (*v_it == v) + return true; + return false; +} + +/** Fill chosen point vector from indices with copies if you are working on a flat torus + * + * IN: W is the point vector + * OUT: landmarks is the output vector + * IN: landmarks_ind is the vector of indices + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) +{ + if (!torus) + for (unsigned j = 0; j < landmarks_ind.size(); ++j) + landmarks.push_back(W[landmarks_ind[j]]); + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + int nbL = landmarks_ind.size(); + // Fill landmarks + for (int i = 0; i < nb_cells-1; ++i) + for (int j = 0; j < nbL; ++j) + { + int cell_i = i; + Point_d point; + for (int l = 0; l < D; ++l) + { + point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1)); + cell_i /= 3; + } + landmarks.push_back(point); + } + } +} + +/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices + * + * IN: t is the Delaunay triangulation + * OUT: full_cells is the output vector + */ + +void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) +{ + // Store vertex indices in a map + int ind = 0; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (t.is_infinite(v_it)) + continue; + else + index_of_vertex[v_it] = ind++; + // Write full cells as vectors in full_cells + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + std::vector<int> cell; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + cell.push_back(index_of_vertex[*v_it]); + full_cells.push_back(cell); + } +} + +//////////////////////////////////////////////////////////////////////////////////////////////////////////// +// IS VIOLATED TEST +//////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Check if a newly created cell is protected from old vertices + * + * t is the Delaunay triangulation + * vertices is the vector containing the point to insert and a facet f in t + * v1 is the vertex of t, such that f and v1 form a simplex + * v2 is the vertex of t, such that f and v2 form another simplex + * delta is the protection constant + * power_protection is true iff the delta-power protection is used + */ + +bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) +{ + assert(vertices.size() == vertices[0].size() || + vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 + assert(v1 != v2); + if (vertices.size() == vertices[0].size() + 1) + // FINITE CASE + { + Sphere_d cs(vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + { + //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + */ + // Check if the simplex is thick enough + Hyperplane_d tau_h(vertices.begin()+1, vertices.end()); + Vector_d orth_tau = tau_h.orthogonal_vector(); + /* + p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1)); + */ + //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n"; + FT orth_length = sqrt(orth_tau.squared_length()); + K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; + // Compute the altitude + FT h = 0; + for (o_it = orth_tau.cartesian_begin(), + p_it = vertices.begin()->cartesian_begin(), + s_it = (vertices.begin()+1)->cartesian_begin(); + o_it != orth_tau.cartesian_end(); + ++o_it, ++p_it, ++s_it) + h += (*o_it)*(*p_it - *s_it)/orth_length; + h = fabs(h); + // Is the center inside the box? + bool inside_the_box = true; + for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) + if (*c_it > 1.0 || *c_it < -1.0) + { + inside_the_box = false; break; + } + if (inside_the_box && h/r < theta0) + return true; + if (!t.is_infinite(v1)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + if (!t.is_infinite(v2)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + else + // INFINITE CASE + { + Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); + while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) + v++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); + Vector_d orth_v = facet_plane.orthogonal_vector(); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + std::vector<FT> coords; + Point_d p = v_it->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!p_is_inside && p_delta_is_inside) + return true; + } + */ + if (!t.is_infinite(v1)) + { + std::vector<FT> coords; + Point_d p = v1->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + if (!t.is_infinite(v2)) + { + std::vector<FT> coords; + Point_d p = v2->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + } + return false; +} + +/** Auxillary recursive function to check if the point p violates the protection of the cell c and + * if there is a violation of an eventual new cell + * + * p is the point to insert + * t is the current triangulation + * c is the current cell (simplex) + * parent_cell is the parent cell (simplex) + * index is the index of the facet between c and parent_cell from parent_cell's point of view + * D is the dimension of the triangulation + * delta is the protection constant + * marked_cells is the vector of all visited cells containing p in their circumscribed ball + * power_protection is true iff you are working with delta-power protection + * + * OUT: true iff inserting p hasn't produced any violation so far + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0) +{ + Euclidean_distance ed; + std::vector<Point_d> vertices; + if (!t.is_infinite(c)) + { + // if the cell is finite, we look if the protection is violated + for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); + FT dist2 = ed.transformed_distance(center_cs, p); + // if the new point is inside the protection ball of a non conflicting simplex + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + // if the new point is inside the circumscribing ball : continue violation searching on neighbours + //if (dist2 < r*r) + //if (dist2 < (5*r+delta)*(5*r+delta)) + if (dist2 < r*r) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) + return true; + } + } + // if the new point is outside the protection sphere + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is guaranteed to be finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + else + { + // Inside of the convex hull is + side. Outside is - side. + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!t.is_infinite(*vh_it)) + vertices.push_back((*vh_it)->point()); + Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); + while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) + v_it++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); + //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); + Vector_d orth_v = facet_plane.orthogonal_vector(); + std::vector<FT> coords; + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + + // If we work with power protection, we just ignore any conflicts + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + //if the cell is infinite we look at the neighbours regardless + if (p_is_inside) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) + return true; + } + } + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is finite if the parent cell is finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + if (!t.is_infinite(parent_cell->vertex(i))) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + //c->tds_data().clear_visited(); + //marked_cells.pop_back(); + return false; +} + +/** Checks if inserting the point p in t will make conflicts + * + * p is the point to insert + * t is the current triangulation + * D is the dimension of triangulation + * delta is the protection constant + * power_protection is true iff you are working with delta-power protection + * OUT: true iff inserting p produces a violation of delta-protection. + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + std::vector<Full_cell_handle> marked_cells; + c = t.locate(p, lt, f, ft, v); + bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0); + for (Full_cell_handle fc : marked_cells) + fc->tds_data().clear(); + return violation_existing_cells; +} + +////////////////////////////////////////////////////////////////////// +// INITIALIZATION +////////////////////////////////////////////////////////////////////// + +void initialize(Search_Tree& W, Delaunay& t, int D, int width, bool torus) +{ + if (!torus) + std::cout << "Non-toric case is not supported\n"; + else + { + if (D == 2) + { + FT stepx = 2.0/width; + FT stepy = sqrt(3)/width; + for (int i = 0; i < width; ++i) + for (int j = 0; j < floor(2*width/sqrt(3)); ++j) + { + insert_delaunay_landmark_with_copies(Point_d(step*i,)) + } + } + else (D == 3) + { + + } + else std::cout << "T^d with d>3 not supported"; + } +} + +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// +//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// + +/////////////////////////////////////////////////////////////////////// +// LANDMARK CHOICE PROCEDURE AS IN PAPER +/////////////////////////////////////////////////////////////////////// + +/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. + * + * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d + * IN: nbP is the size of W + * OUT: landmarks is the output vector for the points + * OUT: landmarks_ind is the output vector for the indices of the selected points in W + * IN: delta is the constant of protection + * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +template<class Search_Tree> +void protected_delaunay_refinement(Search_Tree& W, int nbP, Point_Vector& landmarks, FT delta, bool torus, bool power_protection, FT theta0) +{ + bool return_ = true; + unsigned D = W[0].size(); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + CGAL::Random rand; + int landmark_count = 0; + //std::list<int> index_list; + // shuffle the list of indexes (via a vector) + // { + // std::vector<int> temp_vector; + // for (int i = 0; i < nbP; ++i) + // temp_vector.push_back(i); + // unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + // std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + // //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); + // for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + // index_list.push_front(*it); + // } + if (torus) + if (D == 2) + // \T^2 + { + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 2; ++j) + { + W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0}); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j}); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + } + } + else if (D == 3) + { + + } + //std::cout << "No torus starter available for dim>2\n"; + std::list<int>::iterator list_it = index_list.begin(); + while (list_it != index_list.end()) + { + if (!is_violating_protection(W[*list_it], t, D, delta, power_protection, theta0)) + { + // If no conflicts then insert in every copy of T^3 + + insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus); + if (return_) + { + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + index_list.erase(list_it++); + /* + // PIECE OF CODE FOR DEBUGGING PURPOSES + + Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + if (triangulation_is_protected(t, delta)) + { + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + { //THAT'S WHERE SOMETHING'S WRONG + t.remove(inserted_v); + landmarks_ind.pop_back(); + landmark_count--; + write_delaunay_mesh(t, W[*list_it], is2d); + is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore + } + */ + //std::cout << "index_list_size() = " << index_list.size() << "\n"; + } + else + { + list_it++; + //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; + } + //if (list_it != index_list.end()) + // write_delaunay_mesh(t, W[*list_it], is2d); + } + fill_landmarks(W, landmarks, landmarks_ind, torus); + fill_full_cell_vector(t, full_cells); + /* + if (triangulation_is_protected(t, delta)) + std::cout << "Triangulation is ok\n"; + else + { + std::cout << "Triangulation is BAD!! T_T ã—ãã—ã!\n"; + } + */ + //write_delaunay_mesh(t, W[0], is2d); + //std::cout << t << std::endl; +} + +#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.h b/src/Witness_complex/example/protected_sets/protected_sets_paper.h new file mode 100644 index 00000000..61fcc75b --- /dev/null +++ b/src/Witness_complex/example/protected_sets/protected_sets_paper.h @@ -0,0 +1,917 @@ +#ifndef PROTECTED_SETS_H +#define PROTECTED_SETS_H + +#include <algorithm> +#include <CGAL/Cartesian_d.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Euclidean_distance.h> +#include <CGAL/Kernel_d/Sphere_d.h> +#include <CGAL/Kernel_d/Hyperplane_d.h> +#include <CGAL/Kernel_d/Vector_d.h> + +#include <CGAL/Orthogonal_k_neighbor_search.h> +#include <CGAL/Kd_tree.h> +#include <CGAL/Fuzzy_sphere.h> + +#include <boost/heap/fibonacci_heap.hpp> +#include <boost/heap/policies.hpp> + +#include "output_tikz.h" + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::Point_d Point_d; +typedef K::Line_d Line_d; +typedef K::Vector_d Vector_d; +typedef K::Oriented_side_d Oriented_side_d; +typedef K::Has_on_positive_side_d Has_on_positive_side_d; +typedef K::Sphere_d Sphere_d; +typedef K::Hyperplane_d Hyperplane_d; + +typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; +typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; + +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; + +typedef CGAL::Search_traits_adapter< + std::ptrdiff_t, Point_d*, Traits_base> STraits; +//typedef K TreeTraits; +//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; +//typedef CGAL::Kd_tree<STraits> Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; + + +FT _sfty = pow(10,-14); + +bool experiment1, experiment2 = false; + +/* Experiment 1: epsilon as function on time **********************/ +std::vector<FT> eps_vector; + +/* Experiment 2: R/epsilon on delta *******************************/ +std::vector<FT> epsratio_vector; + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// +// AUXILLARY FUNCTIONS +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well + * + * W is the initial point vector + * chosen_landmark is the index of the chosen point in W + * landmarks_ind is the vector of indices of already chosen points in W + * delaunay is the Delaunay triangulation + * landmark_count is the current number of chosen vertices + * torus is true iff you are working on a flat torus [-1,1]^d + * OUT: Vertex handle to the newly inserted point + */ +Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) +{ + if (!torus) + { + Delaunay_vertex v =delaunay.insert(W[chosen_landmark]); + landmarks_ind.push_back(chosen_landmark); + landmark_count++; + return v; + } + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + Delaunay_vertex v; + for (int i = 0; i < nb_cells; ++i) + { + std::vector<FT> point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + if (i == nb_cells/2) + v = delaunay.insert(point); //v = center point + else + delaunay.insert(point); + } + landmarks_ind.push_back(chosen_landmark); + landmark_count++; + return v; + } +} + +/** Small check if the vertex v is in the full cell fc + */ + +bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) +{ + for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) + if (*v_it == v) + return true; + return false; +} + +/** Fill chosen point vector from indices with copies if you are working on a flat torus + * + * IN: W is the point vector + * OUT: landmarks is the output vector + * IN: landmarks_ind is the vector of indices + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) +{ + if (!torus) + for (unsigned j = 0; j < landmarks_ind.size(); ++j) + landmarks.push_back(W[landmarks_ind[j]]); + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + int nbL = landmarks_ind.size(); + // Fill landmarks + for (int i = 0; i < nb_cells-1; ++i) + for (int j = 0; j < nbL; ++j) + { + int cell_i = i; + Point_d point; + for (int l = 0; l < D; ++l) + { + point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1)); + cell_i /= 3; + } + landmarks.push_back(point); + } + } +} + +/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices + * + * IN: t is the Delaunay triangulation + * OUT: full_cells is the output vector + */ + +void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) +{ + // Store vertex indices in a map + int ind = 0; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (t.is_infinite(v_it)) + continue; + else + index_of_vertex[v_it] = ind++; + // Write full cells as vectors in full_cells + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + std::vector<int> cell; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + cell.push_back(index_of_vertex[*v_it]); + full_cells.push_back(cell); + } +} + +bool sphere_intersects_cube(Point_d& c, FT r) +{ + bool in_cube = true; + // int i = 0, D = p.size(); + for (auto xi = c.cartesian_begin(); xi != c.cartesian_end(); ++xi) + // if ((*xi < 1.0 || *xi > -1.0) && + // (*xi-r < 1.0 || *xi-r > -1.0) && + // (*xi+r < 1.0 || *xi+r > -1.0)) + + if ((*xi-r < -1.0 && *xi+r < -1.0) || + (*xi-r > 1.0 && *xi+r > 1.0 )) + { + in_cube = false; break; + } + return in_cube; +} + +/** Recursive function for checking if the simplex is good, + * meaning it does not contain a k-face, which is not theta0^(k-1) thick + */ + +bool is_theta0_good(std::vector<Point_d>& vertices, FT theta0) +{ + if (theta0 > 1) + { + std::cout << "Warning! theta0 is set > 1\n"; + return false; + } + int D = vertices.size()-1; + if (D <= 1) + return true; // Edges are always good + //******** Circumscribed sphere + Euclidean_distance ed; + Sphere_d cs(vertices.begin(), vertices.end()); + FT r = sqrt(cs.squared_radius()); + for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it) + { + std::vector<Point_d> facet; + for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it) + if (f_it != v_it) + facet.push_back(*f_it); + // Compute the altitude + + if (vertices[0].size() == 3 && D == 2) + { + //Vector_d l = facet[0] - facet[1]; + FT orth_length2 = ed.transformed_distance(facet[0],facet[1]); + K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it; + FT h = 0; + // Scalar product = <sp,l> + FT scalar = 0; + for (p_it = v_it->cartesian_begin(), + s_it = facet[0].cartesian_begin(), + l_it = facet[1].cartesian_begin(); + p_it != v_it->cartesian_end(); + ++l_it, ++p_it, ++s_it) + scalar += (*l_it - *s_it)*(*p_it - *s_it); + // Gram-Schmidt for one vector + for (p_it = v_it->cartesian_begin(), + s_it = facet[0].cartesian_begin(), + l_it = facet[1].cartesian_begin(); + p_it != v_it->cartesian_end(); + ++l_it, ++p_it, ++s_it) + { + FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2; + h += hx*hx; + } + h = sqrt(h); + + if (h/(2*r) < pow(theta0, D-1)) + return false; + if (!is_theta0_good(facet, theta0)) + return false; + } + else + { + Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it); + Vector_d orth_tau = tau_h.orthogonal_vector(); + FT orth_length = sqrt(orth_tau.squared_length()); + K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; + FT h = 0; + for (o_it = orth_tau.cartesian_begin(), + p_it = v_it->cartesian_begin(), + s_it = (facet.begin())->cartesian_begin(); + o_it != orth_tau.cartesian_end(); + ++o_it, ++p_it, ++s_it) + h += (*o_it)*(*p_it - *s_it)/orth_length; + h = fabs(h); + if (h/(2*r) < pow(theta0, D-1)) + return false; + if (!is_theta0_good(facet, theta0)) + return false; + } + } + return true; +} + + +//////////////////////////////////////////////////////////////////////////////////////////////////////////// +// IS VIOLATED TEST +//////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Check if a newly created cell is protected from old vertices + * + * t is the Delaunay triangulation + * vertices is the vector containing the point to insert and a facet f in t + * v1 is the vertex of t, such that f and v1 form a simplex + * v2 is the vertex of t, such that f and v2 form another simplex + * delta is the protection constant + * power_protection is true iff the delta-power protection is used + */ + +bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) +{ + assert(vertices.size() == vertices[0].size() || + vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 + assert(v1 != v2); + if (vertices.size() == vertices[0].size() + 1) + // FINITE CASE + { + Sphere_d cs(vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + { + //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + */ + // Check if the simplex is theta0-good + if (!is_theta0_good(vertices, theta0)) + return true; + // Is the center inside the box? (only Euclidean case) + // if (!torus) + // { + // bool inside_the_box = true; + // for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) + // if (*c_it > 1.0 || *c_it < -1.0) + // { + // inside_the_box = false; break; + // } + // if (inside_the_box && h/r < theta0) + // return true; + // } + // Check the two vertices (if not infinite) + if (!t.is_infinite(v1)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + if (!t.is_infinite(v2)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + else + // INFINITE CASE + { + Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); + while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) + v++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); + Vector_d orth_v = facet_plane.orthogonal_vector(); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + std::vector<FT> coords; + Point_d p = v_it->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!p_is_inside && p_delta_is_inside) + return true; + } + */ + if (!t.is_infinite(v1)) + { + std::vector<FT> coords; + Point_d p = v1->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + if (!t.is_infinite(v2)) + { + std::vector<FT> coords; + Point_d p = v2->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + } + return false; +} + +/** Auxillary recursive function to check if the point p violates the protection of the cell c and + * if there is a violation of an eventual new cell + * + * p is the point to insert + * t is the current triangulation + * c is the current cell (simplex) + * parent_cell is the parent cell (simplex) + * index is the index of the facet between c and parent_cell from parent_cell's point of view + * D is the dimension of the triangulation + * delta is the protection constant + * marked_cells is the vector of all visited cells containing p in their circumscribed ball + * power_protection is true iff you are working with delta-power protection + * + * OUT: true iff inserting p hasn't produced any violation so far + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0) +{ + Euclidean_distance ed; + std::vector<Point_d> vertices; + if (!t.is_infinite(c)) + { + // if the cell is finite, we look if the protection is violated + for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); + FT dist2 = ed.transformed_distance(center_cs, p); + // if the new point is inside the protection ball of a non conflicting simplex + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + // if the new point is inside the circumscribing ball : continue violation searching on neighbours + //if (dist2 < r*r) + //if (dist2 < (5*r+delta)*(5*r+delta)) + if (dist2 < r*r) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) + return true; + } + } + // if the new point is outside the protection sphere + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is guaranteed to be finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + else + { + // Inside of the convex hull is + side. Outside is - side. + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!t.is_infinite(*vh_it)) + vertices.push_back((*vh_it)->point()); + Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); + while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) + v_it++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); + //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); + Vector_d orth_v = facet_plane.orthogonal_vector(); + std::vector<FT> coords; + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + + // If we work with power protection, we just ignore any conflicts + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + //if the cell is infinite we look at the neighbours regardless + if (p_is_inside) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0)) + return true; + } + } + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is finite if the parent cell is finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + if (!t.is_infinite(parent_cell->vertex(i))) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + //c->tds_data().clear_visited(); + //marked_cells.pop_back(); + return false; +} + +/** Checks if inserting the point p in t will make conflicts + * + * p is the point to insert + * t is the current triangulation + * D is the dimension of triangulation + * delta is the protection constant + * power_protection is true iff you are working with delta-power protection + * OUT: true iff inserting p produces a violation of delta-protection. + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + std::vector<Full_cell_handle> marked_cells; + //c = t.locate(p, lt, f, ft, v); + c = t.locate(p); + bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0); + for (Full_cell_handle fc : marked_cells) + fc->tds_data().clear(); + return violation_existing_cells; +} + + +//////////////////////////////////////////////////////////////////////// +// INITIALIZATION +//////////////////////////////////////////////////////////////////////// + +// Query for a sphere near a cite in all copies of a torus +// OUT points_inside +void torus_search(Tree& treeW, int D, Point_d cite, FT r, std::vector<int>& points_inside) +{ + int nb_cells = pow(3, D); + Delaunay_vertex v; + for (int i = 0; i < nb_cells; ++i) + { + std::vector<FT> cite_copy; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + cite_copy.push_back(cite[l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + Fuzzy_sphere fs(cite_copy, r, 0, treeW.traits()); + treeW.search(std::insert_iterator<std::vector<int>>(points_inside, points_inside.end()), fs); + } +} + + +void initialize_torus(Point_Vector& W, Tree& treeW, Delaunay_triangulation& t, FT epsilon, std::vector<int>& landmarks_ind, int& landmark_count) +{ + int D = W[0].size(); + if (D == 2) + { + int xw = 6, yw = 4; + // Triangular lattice close to regular triangles h=0.866a ~ 0.875a : 48p + for (int i = 0; i < xw; ++i) + for (int j = 0; j < yw; ++j) + { + Point_d cite1(std::vector<FT>{2.0/xw*i, 1.0/yw*j}); + std::vector<int> points_inside; + torus_search(treeW, D, cite1, epsilon, points_inside); + assert(points_inside.size() > 0); + insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), + landmarks_ind, t, landmark_count, true); + Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 1.0/yw*(j+0.5)}); + points_inside.clear(); + torus_search(treeW, D, cite2, epsilon, points_inside); + assert(points_inside.size() > 0); + insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), + landmarks_ind, t, landmark_count, true); + } + } + else if (D == 3) + { + int wd = 3; + // Body-centered cubic lattice : 54p + for (int i = 0; i < wd; ++i) + for (int j = 0; j < wd; ++j) + for (int k = 0; k < wd; ++k) + { + Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k}); + std::vector<int> points_inside; + torus_search(treeW, D, cite1, epsilon, points_inside); + assert(points_inside.size() > 0); + insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), + landmarks_ind, t, landmark_count, true); + Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)}); + points_inside.clear(); + torus_search(treeW, D, cite2, epsilon, points_inside); + assert(points_inside.size() > 0); + insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), + landmarks_ind, t, landmark_count, true); + } + } +} + +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// +//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// + +// Struct for R_max_heap elements + +struct R_max_handle +{ + FT value; + Point_d center; + + R_max_handle(FT value_, Point_d c): value(value_), center(c) + {} +}; + +struct R_max_compare +{ + bool operator()(const R_max_handle& rmh1, const R_max_handle& rmh2) const + { + return rmh1.value < rmh2.value; + } +}; + +// typedef boost::heap::fibonacci_heap<R_max_handle, boost::heap::compare<R_max_compare>> Heap; + +// void make_heap(Delaunay_triangulation& t, Heap& R_max_heap) +// { +// R_max_heap.clear(); +// for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) +// { +// if (t.is_infinite(fc_it)) +// continue; +// Point_Vector vertices; +// for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) +// vertices.push_back((*fc_v_it)->point()); +// Sphere_d cs( vertices.begin(), vertices.end()); +// Point_d csc = cs.center(); +// FT r = sqrt(cs.squared_radius()); +// // A ball is in the heap, if it intersects the cube +// bool accepted = sphere_intersects_cube(csc, sqrt(r)); +// if (!accepted) +// continue; +// R_max_heap.push(R_max_handle(r, fc_it, csc)); +// } +// } + +////////////////////////////////////////////////////////////////////////////////////////////////////////// +// SAMPLING RADIUS +////////////////////////////////////////////////////////////////////////////////////////////////////////// + +R_max_handle sampling_radius(Delaunay_triangulation& t) +{ + FT epsilon2 = 0; + Point_d final_center; + Point_d control_point; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin())); + if (epsilon2 < r2) + { + epsilon2 = r2; + final_center = csc; + control_point = (*vertices.begin()); + } + } + return R_max_handle(sqrt(epsilon2), final_center); +} + +/////////////////////////////////////////////////////////////////////// +// LANDMARK CHOICE PROCEDURE +/////////////////////////////////////////////////////////////////////// + +/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. + * + * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d + * IN: nbP is the size of W + * OUT: landmarks is the output vector for the points + * OUT: landmarks_ind is the output vector for the indices of the selected points in W + * IN: delta is the constant of protection + * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void protected_delaunay(Point_Vector& W, + //Point_Vector& landmarks, + std::vector<int>& landmarks_ind, + FT delta, + FT epsilon, + FT alpha, + FT theta0, + //std::vector<std::vector<int>>& full_cells, + bool torus, + bool power_protection + ) +{ + //bool return_ = true; + unsigned D = W[0].size(); + int nbP = W.size(); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + CGAL::Random rand; + int landmark_count = 0; + std::list<int> index_list; + //****************** Kd Tree W + STraits traits(&(W[0])); + Tree treeW(boost::counting_iterator<std::ptrdiff_t>(0), + boost::counting_iterator<std::ptrdiff_t>(nbP), + typename Tree::Splitter(), + traits); + // shuffle the list of indexes (via a vector) + { + std::vector<int> temp_vector; + for (int i = 0; i < nbP; ++i) + temp_vector.push_back(i); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); + for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + index_list.push_front(*it); + } + //******************** Initialize point set + if (!torus) + for (unsigned pos1 = 0; pos1 < D+1; ++pos1) + { + std::vector<FT> point; + for (unsigned i = 0; i < pos1; ++i) + point.push_back(-1); + if (pos1 != D) + point.push_back(1); + for (unsigned i = pos1+1; i < D; ++i) + point.push_back(0); + assert(point.size() == D); + W[index_list.front()] = Point_d(point); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + } + else + initialize_torus(W, treeW, t, epsilon, landmarks_ind, landmark_count); + //std::cout << "Size of treeW: " << treeW.size() << "\n"; + //std::cout << "Size of t: " << t.number_of_vertices() << "\n"; + //******************* Initialize heap for R_max + //Heap R_max_heap; + //make_heap(t, R_max_heap); + + + R_max_handle rh = sampling_radius(t); + FT epsilon0 = rh.value; + if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); + //******************** Iterative algorithm + std::vector<int> candidate_points; + torus_search(treeW, D, + rh.center, + alpha*rh.value, + candidate_points); + std::list<int>::iterator list_it; + std::vector<int>::iterator cp_it = candidate_points.begin(); + while (cp_it != candidate_points.end()) + { + if (!is_violating_protection(W[*cp_it], t, D, delta, power_protection, theta0)) + { + insert_delaunay_landmark_with_copies(W, *cp_it, landmarks_ind, t, landmark_count, torus); + //make_heap(t, R_max_heap); + rh = sampling_radius(t); + if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); + //std::cout << "rhvalue = " << rh.value << "\n"; + //std::cout << "D = " << + candidate_points.clear(); + torus_search(treeW, D, + rh.center, + alpha*rh.value, + candidate_points); + /* + // PIECE OF CODE FOR DEBUGGING PURPOSES + + Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + if (triangulation_is_protected(t, delta)) + { + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + { //THAT'S WHERE SOMETHING'S WRONG + t.remove(inserted_v); + landmarks_ind.pop_back(); + landmark_count--; + write_delaunay_mesh(t, W[*list_it], is2d); + is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore + } + */ + //std::cout << "index_list_size() = " << index_list.size() << "\n"; + } + else + { + cp_it++; + //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; + } + //if (list_it != index_list.end()) + // write_delaunay_mesh(t, W[*list_it], is2d); + } + if (experiment2) epsratio_vector.push_back(rh.value/epsilon0); + std::cout << "The iteration ended when cp_count = " << candidate_points.size() << "\n"; + std::cout << "alphaRmax = " << alpha*rh.value << "\n"; + std::cout << "epsilon' = " << rh.value << "\n"; + std::cout << "nbL = " << landmarks_ind.size() << "\n"; + //fill_landmarks(W, landmarks, landmarks_ind, torus); + //fill_full_cell_vector(t, full_cells); + /* + if (triangulation_is_protected(t, delta)) + std::cout << "Triangulation is ok\n"; + else + { + std::cout << "Triangulation is BAD!! T_T ã—ãã—ã!\n"; + } + */ + //write_delaunay_mesh(t, W[0], is2d); + //std::cout << t << std::endl; +} + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Series of experiments +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + +void start_experiments(Point_Vector& W, FT theta0, std::vector<int>& landmarks_ind, FT epsilon) +{ + // Experiment 1 + experiment1 = true; + protected_delaunay(W, landmarks_ind, 0.1*epsilon, epsilon, 0.5, 0, true, true); + write_tikz_plot(eps_vector,"epstime.tikz"); + experiment1 = false; + + // Experiment 2 + // experiment2 = true; + // for (FT delta = 0; delta < epsilon; delta += 0.1*epsilon) + // protected_delaunay(W, landmarks_ind, delta, epsilon, 0.5, 0, true, true); + // write_tikz_plot(epsratio_vector,"epsratio_delta.tikz"); + // experiment2 = false; + +} + +#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper2.h b/src/Witness_complex/example/protected_sets/protected_sets_paper2.h new file mode 100644 index 00000000..04b5e3bc --- /dev/null +++ b/src/Witness_complex/example/protected_sets/protected_sets_paper2.h @@ -0,0 +1,1384 @@ +#ifndef PROTECTED_SETS_H +#define PROTECTED_SETS_H + +#include <algorithm> +#include <CGAL/Cartesian_d.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Euclidean_distance.h> +#include <CGAL/Kernel_d/Sphere_d.h> +#include <CGAL/Kernel_d/Hyperplane_d.h> +#include <CGAL/Kernel_d/Vector_d.h> + +#include <CGAL/Orthogonal_k_neighbor_search.h> +#include <CGAL/Kd_tree.h> +#include <CGAL/Fuzzy_sphere.h> + +#include <boost/heap/fibonacci_heap.hpp> +#include <boost/heap/policies.hpp> + +#include "output_tikz.h" +#include "../output.h" +#include "../generators.h" + +#include <CGAL/point_generators_d.h> + + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::Point_d Point_d; +typedef K::Line_d Line_d; +typedef K::Vector_d Vector_d; +typedef K::Oriented_side_d Oriented_side_d; +typedef K::Has_on_positive_side_d Has_on_positive_side_d; +typedef K::Sphere_d Sphere_d; +typedef K::Hyperplane_d Hyperplane_d; + +typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; +typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; + +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; + +typedef CGAL::Search_traits_adapter< + std::ptrdiff_t, Point_d*, Traits_base> STraits; +//typedef K TreeTraits; +//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; +//typedef CGAL::Kd_tree<STraits> Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; + +typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; + + +FT _sfty = pow(10,-14); + +bool experiment1, experiment2, experiment3, experiment5 = false; + +/* Experiment 1: epsilon as function on time **********************/ +std::vector<FT> eps_vector; + +/* Experiment 2: R/epsilon on alpha *******************************/ +std::vector<FT> epsratio_vector; +std::vector<FT> epsslope_vector; + +/* Experiment 3: theta on delta ***********************************/ +std::vector<FT> thetamin_vector; FT curr_theta; +std::vector<FT> gammamin_vector; + +/* Statistical data ***********************************************/ +int refused_case1, refused_case2, refused_bad, refused_centers1, refused_centers2; + +void initialize_statistics() +{ + refused_case1 = 0; + refused_case2 = 0; + refused_bad = 0; + refused_centers1 = 0; + refused_centers2 = 0; +} + +void print_statistics() +{ + std::cout << " * Old simplex not protected: " << refused_case1 << "\n"; + std::cout << " * New simplex not protected: " << refused_case2 << "\n"; + std::cout << " * New simplex not good: " << refused_bad << "\n"; + std::cout << " * New-old centers too close: " << refused_centers1 << "\n"; + std::cout << " * New-new centers too close: " << refused_centers2 << "\n"; +} + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// +// AUXILLARY FUNCTIONS +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well + * + * W is the initial point vector + * chosen_landmark is the index of the chosen point in W + * landmarks_ind is the vector of indices of already chosen points in W + * delaunay is the Delaunay triangulation + * landmark_count is the current number of chosen vertices + * torus is true iff you are working on a flat torus [-1,1]^d + * OUT: Vertex handle to the newly inserted point + */ +Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) +{ + if (!torus) + { + Delaunay_vertex v =delaunay.insert(W[chosen_landmark]); + landmarks_ind.push_back(chosen_landmark); + landmark_count++; + return v; + } + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + Delaunay_vertex v; + for (int i = 0; i < nb_cells; ++i) + { + std::vector<FT> point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + if (i == nb_cells/2) + v = delaunay.insert(point); //v = center point + else + delaunay.insert(point); + } + landmarks_ind.push_back(chosen_landmark); + landmark_count++; + return v; + } +} + +/** Small check if the vertex v is in the full cell fc + */ + +bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) +{ + for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) + if (*v_it == v) + return true; + return false; +} + +/** Fill chosen point vector from indices with copies if you are working on a flat torus + * + * IN: W is the point vector + * OUT: landmarks is the output vector + * IN: landmarks_ind is the vector of indices + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) +{ + if (!torus) + for (unsigned j = 0; j < landmarks_ind.size(); ++j) + landmarks.push_back(W[landmarks_ind[j]]); + else + { + int D = W[0].size(); + int nb_cells = pow(3, D); + int nbL = landmarks_ind.size(); + // Fill landmarks + for (int i = 0; i < nb_cells-1; ++i) + for (int j = 0; j < nbL; ++j) + { + int cell_i = i; + Point_d point; + for (int l = 0; l < D; ++l) + { + point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1)); + cell_i /= 3; + } + landmarks.push_back(point); + } + } +} + +/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices + * + * IN: t is the Delaunay triangulation + * OUT: full_cells is the output vector + */ + +void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) +{ + // Store vertex indices in a map + int ind = 0; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (t.is_infinite(v_it)) + continue; + else + index_of_vertex[v_it] = ind++; + // Write full cells as vectors in full_cells + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + std::vector<int> cell; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + cell.push_back(index_of_vertex[*v_it]); + full_cells.push_back(cell); + } +} + +bool sphere_intersects_cube(Point_d& c, FT r) +{ + bool in_cube = true; + // int i = 0, D = p.size(); + for (auto xi = c.cartesian_begin(); xi != c.cartesian_end(); ++xi) + // if ((*xi < 1.0 || *xi > -1.0) && + // (*xi-r < 1.0 || *xi-r > -1.0) && + // (*xi+r < 1.0 || *xi+r > -1.0)) + + if ((*xi-r < -1.0 && *xi+r < -1.0) || + (*xi-r > 1.0 && *xi+r > 1.0 )) + { + in_cube = false; break; + } + return in_cube; +} + +/** Recursive function for checking if the simplex is good, + * meaning it does not contain a k-face, which is not theta0^(k-1) thick + */ + +bool is_theta0_good(std::vector<Point_d>& vertices, FT theta0) +{ + if (theta0 > 1) + { + std::cout << "Warning! theta0 is set > 1\n"; + return false; + } + int D = vertices.size()-1; + if (D <= 1) + return true; // Edges are always good + //******** Circumscribed sphere + Euclidean_distance ed; + Sphere_d cs(vertices.begin(), vertices.end()); + FT r = sqrt(cs.squared_radius()); + for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it) + { + std::vector<Point_d> facet; + for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it) + if (f_it != v_it) + facet.push_back(*f_it); + // Compute the altitude + + if (vertices[0].size() == 3 && D == 2) + { + //Vector_d l = facet[0] - facet[1]; + FT orth_length2 = ed.transformed_distance(facet[0],facet[1]); + K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it; + FT h = 0; + // Scalar product = <sp,l> + FT scalar = 0; + for (p_it = v_it->cartesian_begin(), + s_it = facet[0].cartesian_begin(), + l_it = facet[1].cartesian_begin(); + p_it != v_it->cartesian_end(); + ++l_it, ++p_it, ++s_it) + scalar += (*l_it - *s_it)*(*p_it - *s_it); + // Gram-Schmidt for one vector + for (p_it = v_it->cartesian_begin(), + s_it = facet[0].cartesian_begin(), + l_it = facet[1].cartesian_begin(); + p_it != v_it->cartesian_end(); + ++l_it, ++p_it, ++s_it) + { + FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2; + h += hx*hx; + } + h = sqrt(h); + + if (h/(2*r) < pow(theta0, D-1)) + return false; + if (!is_theta0_good(facet, theta0)) + return false; + } + else + { + Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it); + Vector_d orth_tau = tau_h.orthogonal_vector(); + FT orth_length = sqrt(orth_tau.squared_length()); + K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; + FT h = 0; + for (o_it = orth_tau.cartesian_begin(), + p_it = v_it->cartesian_begin(), + s_it = (facet.begin())->cartesian_begin(); + o_it != orth_tau.cartesian_end(); + ++o_it, ++p_it, ++s_it) + h += (*o_it)*(*p_it - *s_it)/orth_length; + h = fabs(h); + if (experiment3 && thetamin_vector[thetamin_vector.size()-1] > pow(h/(2*r), 1.0/(D-1))) + { + thetamin_vector[thetamin_vector.size()-1] = pow(h/(2*r), 1.0/(D-1)); + //std::cout << "theta=" << h/(2*r) << ", "; + } + if (h/(2*r) < pow(theta0, D-1)) + return false; + if (!is_theta0_good(facet, theta0)) + return false; + } + } + return true; +} + +/** Recursive function for checking the goodness of a simplex, + * meaning it does not contain a k-face, which is not theta0^(k-1) thick + */ + +FT theta(std::vector<Point_d>& vertices) +{ + FT curr_value = 1.0; + int D = vertices.size()-1; + if (D <= 1) + return 1; // Edges are always good + //******** Circumscribed sphere + Euclidean_distance ed; + Sphere_d cs(vertices.begin(), vertices.end()); + FT r = sqrt(cs.squared_radius()); + for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it) + { + std::vector<Point_d> facet; + for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it) + if (f_it != v_it) + facet.push_back(*f_it); + // Compute the altitude + curr_value = std::min(curr_value, theta(facet)); // Check the corresponding facet + if (vertices[0].size() == 3 && D == 2) + { + //Vector_d l = facet[0] - facet[1]; + FT orth_length2 = ed.transformed_distance(facet[0],facet[1]); + K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it; + FT h = 0; + // Scalar product = <sp,l> + FT scalar = 0; + for (p_it = v_it->cartesian_begin(), + s_it = facet[0].cartesian_begin(), + l_it = facet[1].cartesian_begin(); + p_it != v_it->cartesian_end(); + ++l_it, ++p_it, ++s_it) + scalar += (*l_it - *s_it)*(*p_it - *s_it); + // Gram-Schmidt for one vector + for (p_it = v_it->cartesian_begin(), + s_it = facet[0].cartesian_begin(), + l_it = facet[1].cartesian_begin(); + p_it != v_it->cartesian_end(); + ++l_it, ++p_it, ++s_it) + { + FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2; + h += hx*hx; + } + h = sqrt(h); + curr_value = std::min(curr_value, std::pow(h/(2*r), 1.0/(D-1))); + } + else + { + Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it); + Vector_d orth_tau = tau_h.orthogonal_vector(); + FT orth_length = sqrt(orth_tau.squared_length()); + K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; + FT h = 0; + for (o_it = orth_tau.cartesian_begin(), + p_it = v_it->cartesian_begin(), + s_it = (facet.begin())->cartesian_begin(); + o_it != orth_tau.cartesian_end(); + ++o_it, ++p_it, ++s_it) + h += (*o_it)*(*p_it - *s_it)/orth_length; + h = fabs(h); + curr_value = std::min(curr_value, pow(h/(2*r), 1.0/(D-1))); + } + } + return curr_value; +} + +// Doubling in a way 1->2->5->10 +void double_round(int& i) +{ + FT order10 = pow(10,std::floor(std::log10(i))); + int digit = std::floor( i / order10); + std::cout << digit; + if (digit == 1) + i *= 2; + else if (digit == 2) + i = 5*i/2; + else if (digit == 5) + i *= 2; + else + std::cout << "digit not correct. digit = " << digit << std::endl; +} + +//////////////////////////////////////////////////////////////////////////////////////////////////////////// +// IS VIOLATED TEST +//////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** Check if a newly created cell is protected from old vertices + * + * t is the Delaunay triangulation + * vertices is the vector containing the point to insert and a facet f in t + * v1 is the vertex of t, such that f and v1 form a simplex + * v2 is the vertex of t, such that f and v2 form another simplex + * delta is the protection constant + * power_protection is true iff the delta-power protection is used + */ + +bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta0, bool power_protection, FT theta0, FT gamma0) +{ + assert(vertices.size() == vertices[0].size() || + vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 + assert(v1 != v2); + if (vertices.size() == vertices[0].size() + 1) + // FINITE CASE + { + Sphere_d cs(vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + { + //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) + return true; + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) + return true; + } + } + */ + // Is the center inside the box? (only Euclidean case) + // if (!torus) + // { + // bool inside_the_box = true; + // for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) + // if (*c_it > 1.0 || *c_it < -1.0) + // { + // inside_the_box = false; break; + // } + // if (inside_the_box && h/r < theta0) + // return true; + // } + // Check the two vertices (if not infinite) + if (!t.is_infinite(v1)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0)) + { refused_case2++; return true;} + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+r*r*delta0*delta0) + { refused_case2++; return true;} + // Check if the centers are not too close + std::vector<Point_d> sigma(vertices); + sigma[0] = v1->point(); + Sphere_d cs_sigma(sigma.begin(), sigma.end()); + Point_d csc_sigma = cs_sigma.center(); + FT r_sigma = sqrt(cs_sigma.squared_radius()); + FT dcc = sqrt(Euclidean_distance().transformed_distance(center_cs, csc_sigma)); + if (experiment3 && dcc/r < gammamin_vector[gammamin_vector.size()-1]) + gammamin_vector[gammamin_vector.size()-1] = dcc/r; + if (experiment3 && dcc/r_sigma < gammamin_vector[gammamin_vector.size()-1]) + gammamin_vector[gammamin_vector.size()-1] = dcc/r_sigma; + if (dcc < r*gamma0 || dcc < r_sigma*gamma0) + { refused_centers1++; return true; } + } + if (!t.is_infinite(v2)) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0)) + { refused_case2++; return true;} + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+r*r*delta0*delta0) + { refused_case2++; return true;} + // Check if the centers are not too close + std::vector<Point_d> sigma(vertices); + sigma[0] = v2->point(); + Sphere_d cs_sigma(sigma.begin(), sigma.end()); + Point_d csc_sigma = cs_sigma.center(); + FT r_sigma = sqrt(cs_sigma.squared_radius()); + FT dcc = sqrt(Euclidean_distance().transformed_distance(center_cs, csc_sigma)); + if (experiment3 && dcc/r < gammamin_vector[gammamin_vector.size()-1]) + gammamin_vector[gammamin_vector.size()-1] = dcc/r; + if (experiment3 && dcc/r_sigma < gammamin_vector[gammamin_vector.size()-1]) + gammamin_vector[gammamin_vector.size()-1] = dcc/r_sigma; + if (dcc < r*gamma0 || dcc < r_sigma*gamma0) + { refused_centers1++; return true; } + } + // Check if the simplex is theta0-good + if (!is_theta0_good(vertices, theta0)) + { refused_bad++; return true;} + + } + else + // INFINITE CASE + { + Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); + while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) + v++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); + Vector_d orth_v = facet_plane.orthogonal_vector(); + /* + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + std::vector<FT> coords; + Point_d p = v_it->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!p_is_inside && p_delta_is_inside) + return true; + } + */ + if (!t.is_infinite(v1)) + { + std::vector<FT> coords; + Point_d p = v1->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta0 / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + if (!t.is_infinite(v2)) + { + std::vector<FT> coords; + Point_d p = v2->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta0 / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + } + } + return false; +} + +/** Auxillary recursive function to check if the point p violates the protection of the cell c and + * if there is a violation of an eventual new cell + * + * p is the point to insert + * t is the current triangulation + * c is the current cell (simplex) + * parent_cell is the parent cell (simplex) + * index is the index of the facet between c and parent_cell from parent_cell's point of view + * D is the dimension of the triangulation + * delta is the protection constant + * marked_cells is the vector of all visited cells containing p in their circumscribed ball + * power_protection is true iff you are working with delta-power protection + * + * OUT: true iff inserting p hasn't produced any violation so far + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta0, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0, FT gamma0) +{ + Euclidean_distance ed; + std::vector<Point_d> vertices; + if (!t.is_infinite(c)) + { + // if the cell is finite, we look if the protection is violated + for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); + FT dist2 = ed.transformed_distance(center_cs, p); + // if the new point is inside the protection ball of a non conflicting simplex + if (!power_protection) + if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0)) + { refused_case1++; return true;} + if (power_protection) + if (dist2 >= r*r-_sfty && dist2 <= r*r+delta0*delta0*r*r) + { refused_case1++; return true;} + // if the new point is inside the circumscribing ball : continue violation searching on neighbours + //if (dist2 < r*r) + //if (dist2 < (5*r+delta)*(5*r+delta)) + if (dist2 < r*r) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta0, marked_cells, power_protection, theta0, gamma0)) + return true; + } + } + // if the new point is outside the protection sphere + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is guaranteed to be finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta0, power_protection, theta0, gamma0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + else + { + // Inside of the convex hull is + side. Outside is - side. + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!t.is_infinite(*vh_it)) + vertices.push_back((*vh_it)->point()); + Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); + while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) + v_it++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); + //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); + Vector_d orth_v = facet_plane.orthogonal_vector(); + std::vector<FT> coords; + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta0 / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + + // If we work with power protection, we just ignore any conflicts + if (!power_protection && !p_is_inside && p_delta_is_inside) + return true; + //if the cell is infinite we look at the neighbours regardless + if (p_is_inside) + { + c->tds_data().mark_visited(); + marked_cells.push_back(c); + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta0, marked_cells, power_protection, theta0, gamma0)) + return true; + } + } + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is finite if the parent cell is finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + if (!t.is_infinite(parent_cell->vertex(i))) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check = t.infinite_vertex(); + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta0, power_protection, theta0, gamma0)) + //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) + return true; + } + } + //c->tds_data().clear_visited(); + //marked_cells.pop_back(); + return false; +} + +/** Checks if inserting the point p in t will make conflicts + * + * p is the point to insert + * t is the current triangulation + * D is the dimension of triangulation + * delta is the protection constant + * power_protection is true iff you are working with delta-power protection + * OUT: true iff inserting p produces a violation of delta-protection. + */ + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta0, bool power_protection, FT theta0, FT gamma0) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + std::vector<Full_cell_handle> marked_cells; + //c = t.locate(p, lt, f, ft, v); + c = t.locate(p); + bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta0, marked_cells, power_protection, theta0, gamma0); + for (Full_cell_handle fc : marked_cells) + fc->tds_data().clear(); + return violation_existing_cells; +} + + +//////////////////////////////////////////////////////////////////////// +// INITIALIZATION +//////////////////////////////////////////////////////////////////////// + +// Query for a sphere near a cite in all copies of a torus +// OUT points_inside +void torus_search(Tree& treeW, int D, Point_d cite, FT r, std::vector<int>& points_inside) +{ + int nb_cells = pow(3, D); + Delaunay_vertex v; + for (int i = 0; i < nb_cells; ++i) + { + std::vector<FT> cite_copy; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + cite_copy.push_back(cite[l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + Fuzzy_sphere fs(cite_copy, r, 0, treeW.traits()); + treeW.search(std::insert_iterator<std::vector<int>>(points_inside, points_inside.end()), fs); + } +} + + +void initialize_torus(Point_Vector& W, Tree& treeW, Delaunay_triangulation& t, FT epsilon, std::vector<int>& landmarks_ind, int& landmark_count, std::vector<bool>& point_taken) +{ + initialize_statistics(); + int D = W[0].size(); + if (D == 2) + { + int xw = 6, yw = 4; + // Triangular lattice close to regular triangles h=0.866a ~ 0.875a : 48p + for (int i = 0; i < xw; ++i) + for (int j = 0; j < yw; ++j) + { + Point_d cite1(std::vector<FT>{2.0/xw*i, 2.0/yw*j}); + std::vector<int> points_inside; + torus_search(treeW, D, cite1, epsilon, points_inside); + //std::cout << "i=" << i << ", j=" << j << " "; print_vector(points_inside); std::cout << "\n"; + std::vector<int>::iterator p_it = points_inside.begin(); + while (p_it != points_inside.end() && point_taken[*p_it]) + ++p_it; + assert(p_it != points_inside.end()); + //W[*p_it] = cite1; // debug purpose + insert_delaunay_landmark_with_copies(W, *p_it, + landmarks_ind, t, landmark_count, true); + point_taken[*p_it] = true; + + Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 2.0/yw*(j+0.5)}); + points_inside.clear(); + torus_search(treeW, D, cite2, epsilon, points_inside); + //std::cout << "i=" << i << ", j=" << j << " "; print_vector(points_inside); std::cout << "\n"; + p_it = points_inside.begin(); + while (p_it != points_inside.end() && point_taken[*p_it]) + ++p_it; + assert(p_it != points_inside.end()); + //W[*p_it] = cite2; // debug purpose + insert_delaunay_landmark_with_copies(W, *p_it, + landmarks_ind, t, landmark_count, true); + point_taken[*p_it] = true; + } + } + else if (D == 3) + { + int wd = 3; + // Body-centered cubic lattice : 54p + for (int i = 0; i < wd; ++i) + for (int j = 0; j < wd; ++j) + for (int k = 0; k < wd; ++k) + { + Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k}); + std::vector<int> points_inside; + torus_search(treeW, D, cite1, epsilon, points_inside); + std::vector<int>::iterator p_it = points_inside.begin(); + while (p_it != points_inside.end() && point_taken[*p_it]) + ++p_it; + assert(p_it != points_inside.end()); + insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), + landmarks_ind, t, landmark_count, true); + point_taken[*p_it] = true; + + Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)}); + points_inside.clear(); + torus_search(treeW, D, cite2, epsilon, points_inside); + p_it = points_inside.begin(); + while (p_it != points_inside.end() && point_taken[*p_it]) + ++p_it; + assert(p_it != points_inside.end()); + insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), + landmarks_ind, t, landmark_count, true); + point_taken[*p_it] = true; + } + } + //write_mesh +} + +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// +//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// +/////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////// + +// Struct for R_max_heap elements + +struct R_max_handle +{ + FT value; + Point_d center; + + R_max_handle(FT value_, Point_d c): value(value_), center(c) + {} +}; + +struct R_max_compare +{ + bool operator()(const R_max_handle& rmh1, const R_max_handle& rmh2) const + { + return rmh1.value < rmh2.value; + } +}; + +// typedef boost::heap::fibonacci_heap<R_max_handle, boost::heap::compare<R_max_compare>> Heap; + +// void make_heap(Delaunay_triangulation& t, Heap& R_max_heap) +// { +// R_max_heap.clear(); +// for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) +// { +// if (t.is_infinite(fc_it)) +// continue; +// Point_Vector vertices; +// for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) +// vertices.push_back((*fc_v_it)->point()); +// Sphere_d cs( vertices.begin(), vertices.end()); +// Point_d csc = cs.center(); +// FT r = sqrt(cs.squared_radius()); +// // A ball is in the heap, if it intersects the cube +// bool accepted = sphere_intersects_cube(csc, sqrt(r)); +// if (!accepted) +// continue; +// R_max_heap.push(R_max_handle(r, fc_it, csc)); +// } +// } + +////////////////////////////////////////////////////////////////////////////////////////////////////////// +// SAMPLING RADIUS +////////////////////////////////////////////////////////////////////////////////////////////////////////// + +R_max_handle sampling_radius(Delaunay_triangulation& t) +{ + FT epsilon2 = 0; + Point_d final_center; + Point_d control_point; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin())); + if (epsilon2 < r2) + { + epsilon2 = r2; + final_center = csc; + control_point = (*vertices.begin()); + } + } + return R_max_handle(sqrt(epsilon2), final_center); +} + +FT sampling_fatness(Delaunay_triangulation& t) +{ + FT curr_theta = 1.0; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + FT theta_f = theta(vertices); + curr_theta = std::min(curr_theta, theta_f); + //std::cout << "theta(sigma) = " << theta_f << "\n"; + } + return curr_theta; +} + +// Generate an epsilon sample for a given epsilon +void generate_epsilon_sample_torus(Point_Vector& W, FT epsilon, int dim, Delaunay_triangulation& t) +{ + W.clear(); + t.clear(); + int point_count = 0; + std::vector<int> point_ind; + // std::vector<FT> coords; + FT curr_eps = 2*dim; + // Initialize + // for (int i = 0; i < dim; ++i) + // coords.push_back(-1); + // R_max_handle rmh(2*sqrt(dim), Point_d(coords)); + // int N = dim; std::floor(std::pow(1/epsilon,dim)); + // std::cout << N << "\n"; + typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; + Random_cube_iterator rp(dim, 1.0); + W.push_back(*rp++); + insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); + curr_eps = sampling_radius(t).value; + while (curr_eps > epsilon) + { + + W.push_back(*rp++); + insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); + + Point_d c = sampling_radius(t).center; + W.push_back(c); + insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); + curr_eps = sampling_radius(t).value; + + std::cout << "curr_eps = " << curr_eps << "\n"; + } + // Iterate and insert in a torus + // while (rmh.value > epsilon) + // { + // W.push_back(rmh.center); + // insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); + // rmh = sampling_radius(t); + // //std::cout << rmh.value; + // } +} + +/////////////////////////////////////////////////////////////////////// +// LANDMARK CHOICE PROCEDURE +/////////////////////////////////////////////////////////////////////// + +/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. + * + * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d + * IN: nbP is the size of W + * OUT: landmarks is the output vector for the points + * OUT: landmarks_ind is the output vector for the indices of the selected points in W + * IN: delta is the constant of protection + * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation + * IN: torus is true iff you are working on a flat torus [-1,1]^d + */ + +void protected_delaunay(Point_Vector& W, + //Point_Vector& landmarks, + std::vector<int>& landmarks_ind, + FT alpha, + FT epsilon, + FT delta0, + FT theta0, + FT gamma0, + //std::vector<std::vector<int>>& full_cells, + bool torus, + bool power_protection + ) +{ + //bool return_ = true; + unsigned D = W[0].size(); + int nbP = W.size(); + //FT beta = 1/(1-alpha); + //FT Ad = pow((4*alpha + 8*beta)/alpha, D); + //FT theta0 = 1/Ad; + //FT delta0 = pow(1/Ad,D); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + std::vector<bool> point_taken(nbP,false); + CGAL::Random rand; + int landmark_count = 0; + std::list<int> index_list; + //****************** Kd Tree W + STraits traits(&(W[0])); + Tree treeW(boost::counting_iterator<std::ptrdiff_t>(0), + boost::counting_iterator<std::ptrdiff_t>(nbP), + typename Tree::Splitter(), + traits); + // shuffle the list of indexes (via a vector) + { + std::vector<int> temp_vector; + for (int i = 0; i < nbP; ++i) + temp_vector.push_back(i); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); + for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + index_list.push_front(*it); + } + //******************** Initialize point set + if (!torus) + for (unsigned pos1 = 0; pos1 < D+1; ++pos1) + { + std::vector<FT> point; + for (unsigned i = 0; i < pos1; ++i) + point.push_back(-1); + if (pos1 != D) + point.push_back(1); + for (unsigned i = pos1+1; i < D; ++i) + point.push_back(0); + assert(point.size() == D); + W[index_list.front()] = Point_d(point); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); + index_list.pop_front(); + } + else + initialize_torus(W, treeW, t, epsilon, landmarks_ind, landmark_count, point_taken); + //std::cout << "Size of treeW: " << treeW.size() << "\n"; + //std::cout << "Size of t: " << t.number_of_vertices() << "\n"; + //******************* Initialize heap for R_max + //Heap R_max_heap; + //make_heap(t, R_max_heap); + + + R_max_handle rh = sampling_radius(t); + FT epsilon0 = rh.value; + if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); + //******************** Iterative algorithm + std::vector<int> candidate_points; + torus_search(treeW, D, + rh.center, + alpha*rh.value, + candidate_points); + std::list<int>::iterator list_it; + std::vector<int>::iterator cp_it = candidate_points.begin(); + while (cp_it != candidate_points.end()) + { + if (!point_taken[*cp_it] && !is_violating_protection(W[*cp_it], t, D, delta0, power_protection, theta0, gamma0)) + { + Delaunay_vertex v = insert_delaunay_landmark_with_copies(W, *cp_it, landmarks_ind, t, landmark_count, torus); + { + // Simple check if the new cells don't have centers too close one to another + std::vector<Full_cell_handle> inc_cells; + std::back_insert_iterator<std::vector<Full_cell_handle>> out(inc_cells); + t.tds().incident_full_cells(v, out); + + std::vector<Sphere_d> spheres; + for (auto i_it = inc_cells.begin(); i_it != inc_cells.end(); ++i_it) + { + std::vector<Point_d> vertices; + for (auto v_it = (*i_it)->vertices_begin(); v_it != (*i_it)->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + spheres.push_back(Sphere_d(vertices.begin(), vertices.end())); + } + for (auto s_it = spheres.begin(); s_it != spheres.end(); ++s_it) + for (auto t_it = s_it+1; t_it != spheres.end(); ++t_it) + { + FT ddc2 = ed.transformed_distance(s_it->center(),t_it->center()); + if (ddc2 < gamma0*gamma0*s_it->squared_radius() || + ddc2 < gamma0*gamma0*t_it->squared_radius()) + { refused_centers2++; } + } + } + + //std::cout << *cp_it << ",\n"; + //make_heap(t, R_max_heap); + point_taken[*cp_it] = true; + rh = sampling_radius(t); + if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); + //std::cout << "rhvalue = " << rh.value << "\n"; + //std::cout << "D = " << + candidate_points.clear(); + torus_search(treeW, D, + rh.center, + alpha*rh.value, + candidate_points); + cp_it = candidate_points.begin(); + /* + // PIECE OF CODE FOR DEBUGGING PURPOSES + + Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + if (triangulation_is_protected(t, delta)) + { + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + { //THAT'S WHERE SOMETHING'S WRONG + t.remove(inserted_v); + landmarks_ind.pop_back(); + landmark_count--; + write_delaunay_mesh(t, W[*list_it], is2d); + is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore + } + */ + //std::cout << "index_list_size() = " << index_list.size() << "\n"; + } + else + { + cp_it++; + //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; + } + //if (list_it != index_list.end()) + // write_delaunay_mesh(t, W[*list_it], is2d); + } + + if (experiment2) epsratio_vector.push_back(rh.value/epsilon0); + if (experiment2) epsslope_vector.push_back( (pow(1/rh.value,D)-pow(1/epsilon0,D))/(landmarks_ind.size() - 48) ); + std::cout << "The iteration ended when cp_count = " << candidate_points.size() << "\n"; + std::cout << "alphaRmax = " << alpha*rh.value << "\n"; + std::cout << "epsilon' = " << rh.value << "\n"; + std::cout << "nbL = " << landmarks_ind.size() << "\n"; + print_statistics(); + //print_vector(landmarks_ind); std::cout << std::endl; + //std::sort(landmarks_ind.begin(), landmarks_ind.end()); + print_vector(landmarks_ind); std::cout << std::endl; + if (experiment3) thetamin_vector[thetamin_vector.size()-1] = sampling_fatness(t); + std::cout << "theta = " << sampling_fatness(t) << "\n"; + //fill_landmarks(W, landmarks, landmarks_ind, torus); + //fill_full_cell_vector(t, full_cells); + /* + if (triangulation_is_protected(t, delta)) + std::cout << "Triangulation is ok\n"; + else + { + std::cout << "Triangulation is BAD!! T_T ã—ãã—ã!\n"; + } + */ + write_delaunay_mesh(t, W[0], true); + //std::cout << t << std::endl; +} + +void run_experiment5(Point_Vector& W, + int D, + FT alpha, + FT epsilon, + FT delta0, + FT theta0, + FT gamma0, + //std::vector<std::vector<int>>& full_cells, + bool torus, + bool power_protection + ) +{ + // INITIALIZATION + Delaunay_triangulation t(D); + std::vector<int> landmarks_ind; + int landmark_count = 0; + initialize_statistics(); + if (D == 2) + { + int xw = 6, yw = 4; + // Triangular lattice close to regular triangles h=0.866a ~ 0.875a : 48p + for (int i = 0; i < xw; ++i) + for (int j = 0; j < yw; ++j) + { + Point_d cite1(std::vector<FT>{2.0/xw*i, 2.0/yw*j}); + W.push_back(cite1); // debug purpose + insert_delaunay_landmark_with_copies(W, W.size()-1, + landmarks_ind, t, landmark_count, true); + + Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 2.0/yw*(j+0.5)}); + W.push_back(cite2); // debug purpose + insert_delaunay_landmark_with_copies(W, W.size()-1, + landmarks_ind, t, landmark_count, true); + } + } + else if (D == 3) + { + int wd = 3; + // Body-centered cubic lattice : 54p + for (int i = 0; i < wd; ++i) + for (int j = 0; j < wd; ++j) + for (int k = 0; k < wd; ++k) + { + Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k}); + W.push_back(cite1); // debug purpose + insert_delaunay_landmark_with_copies(W, W.size()-1, + landmarks_ind, t, landmark_count, true); + + Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)}); + W.push_back(cite2); // debug purpose + insert_delaunay_landmark_with_copies(W, W.size()-1, + landmarks_ind, t, landmark_count, true); + } + } + + // ITERATIONS + R_max_handle rh = sampling_radius(t); + Point_d rp = *(Random_point_iterator(D, alpha*rh.value)); + int death_count = 0; + std::cout << "death count " << death_count << " rp = " << rp << "\n"; + while (death_count < 100) + { + std::vector<FT> coords; + for (auto c_it = rh.center.cartesian_begin(), + r_it = rp.cartesian_begin(); + c_it != rh.center.cartesian_end(); + ++c_it, ++r_it) + coords.push_back(*c_it + *r_it); + Point_d new_p(coords); + if (!is_violating_protection(new_p, t, D, delta0, power_protection, theta0, gamma0)) + { + W.push_back(new_p); + insert_delaunay_landmark_with_copies(W, W.size()-1, landmarks_ind, t, landmark_count, torus); + rh = sampling_radius(t); + rp = *(Random_point_iterator(D, alpha*rh.value)); + death_count = 0; + std::cout << "death count " << death_count << " rp = " << rp << "\n"; + } + else + { + rp = *(Random_point_iterator(D, alpha*rh.value)); + death_count++; + std::cout << "death count " << death_count << " rp = " << rp << "\n"; + } + //Point_d new_p = (*rp++) + Vector_d; + } +} + +/////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Series of experiments +/////////////////////////////////////////////////////////////////////////////////////////////////////////// + +void start_experiments(Point_Vector& W, FT alpha, std::vector<int>& landmarks_ind, FT epsilon) +{ + int experiment_no = 1; + FT delta0 = 0.1; + FT theta0 = 0.1; + FT gamma0 = 0.01; + std::string suffix; + //std::cout << "よã†ã“ãジプシー我ãŒç¥žç§˜ã®éƒ¨å±‹ã¸:\n"; + while (experiment_no != 0) + { + std::cout << "Enter experiment no (0 to exit): "; + std::cin >> experiment_no; + switch (experiment_no) + { + case 1: + // Experiment 1 + experiment1 = true; + eps_vector = {}; + std::cout << "Enter delta0: "; std::cin >> delta0; + std::cout << "Enter theta0: "; std::cin >> theta0; + std::cout << "Enter gamma0: "; std::cin >> gamma0; + protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true); + write_tikz_plot(eps_vector,"epstime.tikz"); + experiment1 = false; + break; + + case 2: + // Experiment 2 + suffix = ""; + experiment2 = true; + epsratio_vector = {0}; + epsslope_vector = {0}; + std::cout << "File name suffix: "; + std::cin >> suffix; + for (FT alpha = 0.01; alpha < 0.999; alpha += 0.01) + { + landmarks_ind.clear(); + std::cout << "Test for alpha = " << alpha << "\n"; + protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true); + } + write_tikz_plot(epsratio_vector,"epsratio_alpha." + suffix + ".tex"); + write_tikz_plot(epsslope_vector,"epsslope_alpha." + suffix + ".tex"); + experiment2 = false; + break; + + case 3: + // Experiment 3 + experiment3 = true; + thetamin_vector = {}; + gammamin_vector = {}; + theta0 = 0; + gamma0 = 0; + for (FT delta0 = 0; delta0 < 0.999; delta0 += 0.05) + { + landmarks_ind.clear(); + thetamin_vector.push_back(1.0); //0.7489 fatness of the initialization + gammamin_vector.push_back(10); + std::cout << "Test for delta0 = " << delta0 << "\n"; + protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true); + } + write_tikz_plot(thetamin_vector,"thetamin_delta.tex"); + write_tikz_plot(gammamin_vector,"gammamin_delta.tex"); + experiment3 = false; + break; + + // case 4: + // // Experiment 4 + // { + // int dim; + // std::cout << "Enter dimension: "; + // std::cin >> dim; + // Delaunay_triangulation t(dim); + // // for (FT eps = 0.7; eps < 1.1; eps += 0.1) + // // { + // // generate_epsilon_sample_torus(W, eps, dim, t); + // // for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + // // { + // // if (t.is_infinite(v_it)) + // // continue; + // // bool in_cube = true; + // // for (auto xi = v_it->cartesian_begin(); xi != v_it->cartesian_end(); ++xi) + // // if (*xi > 1.0 || *xi < -1.0) + // // { + // // in_cube = false; break; + // // } + // // if (!in_cube) + // // continue; + // // for (auto t.tds().incident_full_cells()) + // // } + // // std::cout << "eps = " << eps << ", real epsilon = " << sampling_radius(t).value << "\n"; + // // } + // // } + // break; + + + case 5: + // Experiment 5 + experiment5 = true; + // std::cout << "Enter dimension: "; + // std::cin >> dim; + + landmarks_ind.clear(); + W.clear(); + run_experiment5(W, alpha, epsilon, delta0, theta0, gamma0, true, true); + experiment5 = false; + break; + } + + } + +} + +#endif diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp index a9a2959b..e448c55d 100644 --- a/src/Witness_complex/example/witness_complex_cube.cpp +++ b/src/Witness_complex/example/witness_complex_cube.cpp @@ -20,6 +20,11 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +// Avoiding the max arity issue with CGAL +#ifndef BOOST_PARAMETER_MAX_ARITY +# define BOOST_PARAMETER_MAX_ARITY 12 +#endif + #include <iostream> #include <fstream> #include <ctime> @@ -37,6 +42,10 @@ #include "gudhi/Witness_complex.h" #include "gudhi/reader_utils.h" #include "Torus_distance.h" +#include "generators.h" +#include "output.h" +//#include "protected_sets/protected_sets.h" +#include "protected_sets/protected_sets_paper2.h" #include <CGAL/Cartesian_d.h> #include <CGAL/Search_traits.h> @@ -106,8 +115,6 @@ typedef std::vector<Point_d> Point_Vector; //typedef K::Equal_d Equal_d; //typedef CGAL::Random_points_in_cube_d<CGAL::Point_d<CGAL::Cartesian_d<FT> > > Random_cube_iterator; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; typedef Delaunay_triangulation::Facet Facet; @@ -117,449 +124,84 @@ typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; typedef K::Sphere_d Sphere_d; typedef K::Hyperplane_d Hyperplane_d; +/*////////////////////////////////////// + * GLOBAL VARIABLES ******************** + *////////////////////////////////////// -bool toric=false; - - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , Point_Vector & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - Point_d p(point.begin(), point.end()); - if (point.size() != 1) - points.push_back(p); - } - in_file.close(); -} - -void generate_points_grid(Point_Vector& W, int width, int D) -{ - int nb_points = 1; - for (int i = 0; i < D; ++i) - nb_points *= width; - for (int i = 0; i < nb_points; ++i) - { - std::vector<double> point; - int cell_i = i; - for (int l = 0; l < D; ++l) - { - point.push_back(0.01*(cell_i%width)); - cell_i /= width; - } - W.push_back(point); - } -} - -void generate_points_random_box(Point_Vector& W, int nbP, int dim) -{ - /* - Random_cube_iterator rp(dim, 1); - for (int i = 0; i < nbP; i++) - { - std::vector<double> point; - for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it) - point.push_back(*it); - W.push_back(Point_d(point)); - rp++; - } - */ - Random_cube_iterator rp(dim, 1.0); - for (int i = 0; i < nbP; i++) - { - W.push_back(*rp++); - } -} - - -void write_wl( std::string file_name, std::vector< std::vector <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} +//NA bool toric=false; +bool power_protection = true; +bool grid_points = true; +bool is2d = true; +//FT _sfty = pow(10,-14); +bool torus = false; -void write_points( std::string file_name, std::vector< Point_d > & points) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : points) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); -} - -void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); -} - - -void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count) -{ - delaunay.insert(W[chosen_landmark]); - landmarks_ind.push_back(chosen_landmark); - landmark_count++; -} - -bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) -{ - for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) - if (*v_it == v) - return true; - return false; -} - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, bool is_infinite, const Point_d& p, FT delta) -{ - if (!is_infinite) - // FINITE CASE - { - Sphere_d cs(vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) - { - //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); - if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); - //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) - if (dist2 >= r*r && dist2 <= r*r+delta*delta) - return true; - } - } - } - else - // INFINITE CASE - { - Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); - while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) - v++; - Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); - Vector_d orth_v = facet_plane.orthogonal_vector(); - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) - if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) - { - std::vector<FT> coords; - Point_d p = v_it->point(); - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - if (!p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells) -{ - Euclidean_distance ed; - std::vector<Point_d> vertices; - if (!t.is_infinite(c)) - { - // if the cell is finite, we look if the protection is violated - for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) - vertices.push_back((*v_it)->point()); - Sphere_d cs( vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); - FT dist2 = ed.transformed_distance(center_cs, p); - // if the new point is inside the protection ball of a non conflicting simplex - //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) - if (dist2 >= r*r && dist2 <= r*r+delta*delta) - return true; - c->tds_data().mark_visited(); - marked_cells.push_back(c); - // if the new point is inside the circumscribing ball : continue violation searching on neughbours - if (dist2 < r*r) - for (int i = 0; i < D+1; ++i) - { - Full_cell_handle next_c = c->neighbor(i); - if (next_c->tds_data().is_clear() && - is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells)) - return true; - } - // if the new point is outside the protection sphere - else - { - // facet f is on the border of the conflict zone : check protection of simplex {p,f} - // the new simplex is guaranteed to be finite - vertices.clear(); vertices.push_back(p); - for (int i = 0; i < D+1; ++i) - if (i != index) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check; - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!vertex_is_in_full_cell(*vh_it, parent_cell)) - { - vertex_to_check = *vh_it; break; - } - if (new_cell_is_violated(t, vertices, false, vertex_to_check->point(), delta)) - return true; - } - } - else - { - // Inside of the convex hull is + side. Outside is - side. - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!t.is_infinite(*vh_it)) - vertices.push_back((*vh_it)->point()); - Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); - while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) - v_it++; - Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); - //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); - Vector_d orth_v = facet_plane.orthogonal_vector(); - std::vector<FT> coords; - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - if (!p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - c->tds_data().mark_visited(); - marked_cells.push_back(c); - if (p_is_inside) - for (int i = 0; i < D+1; ++i) - { - Full_cell_handle next_c = c->neighbor(i); - if (next_c->tds_data().is_clear() && - is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells)) - return true; - } - else - { - // facet f is on the border of the conflict zone : check protection of simplex {p,f} - // the new simplex is finite if the parent cell is finite - vertices.clear(); vertices.push_back(p); - bool new_simplex_is_finite = false; - for (int i = 0; i < D+1; ++i) - if (i != index) - { - if (t.is_infinite(parent_cell->vertex(i))) - new_simplex_is_finite = true; - else - vertices.push_back(parent_cell->vertex(i)->point()); - } - Delaunay_vertex vertex_to_check; - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!vertex_is_in_full_cell(*vh_it, parent_cell)) - { - vertex_to_check = *vh_it; break; - } - if (new_cell_is_violated(t, vertices, new_simplex_is_finite, vertex_to_check->point(), delta)) - return true; - } - } - return false; -} - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) -{ - Euclidean_distance ed; - Delaunay_triangulation::Vertex_handle v; - Delaunay_triangulation::Face f(t.current_dimension()); - Delaunay_triangulation::Facet ft; - Delaunay_triangulation::Full_cell_handle c; - Delaunay_triangulation::Locate_type lt; - std::vector<Full_cell_handle> marked_cells; - c = t.locate(p, lt, f, ft, v); - bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - -bool old_is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) +bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) { + std::cout << "Start protection verification\n"; Euclidean_distance ed; - Delaunay_triangulation::Vertex_handle v; - Delaunay_triangulation::Face f(t.current_dimension()); - Delaunay_triangulation::Facet ft; - Delaunay_triangulation::Full_cell_handle c; - Delaunay_triangulation::Locate_type lt; - c = t.locate(p, lt, f, ft, v); - for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) - if (!t.is_infinite(fc_it)) - { - std::vector<Point_d> vertices; - for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) - vertices.push_back((*v_it)->point()); - Sphere_d cs( vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); - FT dist2 = ed.transformed_distance(center_cs, p); - //if the new point is inside the protection ball of a non conflicting simplex - if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) - return true; - } - t.insert(p, c); - return false; -} - -void write_delaunay_mesh(Delaunay_triangulation& t, const Point_d& p) -{ - std::ofstream ofs ("delaunay.mesh", std::ofstream::out); - int nbV = t.number_of_vertices()+1; - ofs << "MeshVersionFormatted 1\nDimension 2\n"; - ofs << "Vertices\n" << nbV << "\n"; - int ind = 1; //index of a vertex + // Fill the map Vertices -> Numbers std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + int ind = 0; for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) { if (t.is_infinite(v_it)) continue; - for (auto coord = v_it->point().cartesian_begin(); coord != v_it->point().cartesian_end(); ++coord) - ofs << *coord << " "; - ofs << "508\n"; index_of_vertex[v_it] = ind++; } - for (auto coord = p.cartesian_begin(); coord != p.cartesian_end(); ++coord) - ofs << *coord << " "; - ofs << "208\n"; - /* - int nbFacets = 0; - for (auto ft_it = t.finite_facets_begin(); ft_it != t.finite_facets_end(); ++ft_it) - nbFacets++; - ofs << "\nEdges\n" << nbFacets << "\n\n"; - for (auto ft_it = t.facets_begin(); ft_it != t.facets_end(); ++ft_it) - { - if (t.is_infinite(ft_it)) - continue; - for (auto vh_it = ft_it->vertices_begin(); vh_it != ft_it->vertices_end(); ++vh_it) - ofs << index_of_vertex[*vh_it] << " "; - } - */ - ofs << "Triangles " << t.number_of_finite_full_cells()+1 << "\n"; - for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) - { - if (t.is_infinite(fc_it)) - continue; - for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) - ofs << index_of_vertex[*vh_it] << " "; - ofs << "508\n"; - } - ofs << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; - ofs << "End\n"; - ofs.close(); -} - -bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) -{ - // Verification part - Euclidean_distance ed; for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) if (!t.is_infinite(fc_it)) - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) + { + std::vector<Point_d> vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(0)->point())); + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) //check if vertex belongs to the face - if (!vertex_is_in_full_cell(v_it, fc_it)) - { - std::vector<Point_d> vertices; - for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) - vertices.push_back((*fc_v_it)->point()); - Sphere_d cs( vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(0)->point())); - FT dist2 = ed.transformed_distance(center_cs, v_it->point()); - //if the new point is inside the protection ball of a non conflicting simplex - //std::cout << "Dist^2 = " << dist2 << " (r+delta)*(r+delta) = " << (r+delta)*(r+delta) << " r^2 = " << r*r <<"\n"; - //if (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r) - if (dist2 <= r*r+delta*delta && dist2 >= r*r) - { - write_delaunay_mesh(t, v_it->point()); - std::cout << "Problematic vertex " << *v_it << " "; - std::cout << "Problematic cell " << *fc_it << "\n"; - std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", r^2+delta^2 = " << r*r+delta*delta << "\n"; - return false; - } - } - + if (!vertex_is_in_full_cell(v_it, fc_it)) + { + FT dist2 = ed.transformed_distance(center_cs, v_it->point()); + //if the new point is inside the protection ball of a non conflicting simplex + //std::cout << "Dist^2 = " << dist2 << " (r+delta)*(r+delta) = " << (r+delta)*(r+delta) << " r^2 = " << r*r <<"\n"; + if (!power_protection) + if (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r) + { + write_delaunay_mesh(t, v_it->point(), is2d); + // Output the problems + std::cout << "Problematic vertex " << index_of_vertex[v_it] << " "; + std::cout << "Problematic cell "; + for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) + if (!t.is_infinite(*vh_it)) + std::cout << index_of_vertex[*vh_it] << " "; + std::cout << "\n"; + std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", (r+delta)^2 = " << (r+delta)*(r+delta) << "\n"; + return false; + } + if (power_protection) + if (dist2 <= r*r+delta*delta && dist2 >= r*r) + { + write_delaunay_mesh(t, v_it->point(), is2d); + std::cout << "Problematic vertex " << *v_it << " "; + std::cout << "Problematic cell " << *fc_it << "\n"; + std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", r^2+delta^2 = " << r*r+delta*delta << "\n"; + return false; + } + } + } return true; } -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind) -{ - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); -} - -void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) -{ - // Store vertex indices in a map - int ind = 0; //index of a vertex - std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (t.is_infinite(v_it)) - continue; - else - index_of_vertex[v_it] = ind++; - // Write full cells as vectors in full_cells - for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) - { - if (t.is_infinite(fc_it)) - continue; - std::vector<int> cell; - for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) - cell.push_back(index_of_vertex[*v_it]); - full_cells.push_back(cell); - } -} +////////////////////////////////////////////////////////////////////////////////////////////////////////// +// SAMPLING RADIUS +////////////////////////////////////////////////////////////////////////////////////////////////////////// -FT sampling_radius(Delaunay_triangulation& t) +FT sampling_radius(Delaunay_triangulation& t, FT epsilon0) { - FT epsilon2 = 4.0; + FT epsilon2 = 0; + Point_d control_point; for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) { if (t.is_infinite(fc_it)) @@ -578,134 +220,106 @@ FT sampling_radius(Delaunay_triangulation& t) if (!in_cube) continue; FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin())); - if (epsilon2 > r2) - epsilon2 = r2; + if (epsilon2 < r2) + { + epsilon2 = r2; + control_point = (*vertices.begin()); + } + } + if (epsilon2 < epsilon0*epsilon0) + { + std::cout << "ACHTUNG! E' < E\n"; + std::cout << "eps = " << epsilon0 << " eps' = " << sqrt(epsilon2) << "\n"; + write_delaunay_mesh(t, control_point, is2d); } return sqrt(epsilon2); } -FT point_sampling_radius_by_delaunay(Point_Vector& points) +FT point_sampling_radius_by_delaunay(Point_Vector& points, FT epsilon0) { Delaunay_triangulation t(points[0].size()); t.insert(points.begin(), points.end()); - return sampling_radius(t); + return sampling_radius(t, epsilon0); } -void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT delta, std::vector<std::vector<int>>& full_cells) +// A little script to make a tikz histogram of epsilon distribution +// Returns the average epsilon +FT epsilon_histogram(Delaunay_triangulation& t, int n) { - unsigned D = W[0].size(); - Torus_distance td; - Euclidean_distance ed; - Delaunay_triangulation t(D); - CGAL::Random rand; - int landmark_count = 0; - std::list<int> index_list; - // shuffle the list of indexes (via a vector) - { - std::vector<int> temp_vector; - for (int i = 0; i < nbP; ++i) - temp_vector.push_back(i); - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); - //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); - for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) - index_list.push_front(*it); - } - for (unsigned pos1 = 0; pos1 < D+1; ++pos1) - { - std::vector<FT> point; - for (unsigned i = 0; i < pos1; ++i) - point.push_back(-1); - if (pos1 != D) - point.push_back(1); - for (unsigned i = pos1+1; i < D; ++i) - point.push_back(0); - assert(point.size() == D); - W[index_list.front()] = Point_d(point); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - } - // add the first D+1 vertices to form one finite cell - /* - for (int i = 0; i <= D+1; ++i) - { - t.insert(W[index_list.front()]); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - } - */ - /* - { - std::vector<FT> coords; - for (int i = 0; i < D; ++i) - coords.push_back(-1); - W[index_list.front()] = Point_d(coords); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - for (int i = 0; i < D; ++i) - { - coords.clear(); - for (int j = 0; j < D; ++j) - if (i == j) - coords.push_back(1); - else - coords.push_back(-1); - W[index_list.front()] = Point_d(coords); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - } - } - */ - //std::cout << t; - //assert(t.number_of_vertices() == D+1); - //assert(landmarks_ind.size() == D+1); - //assert(W[landmarks_ind[0]][0] == 0); - // add other vertices if they don't violate protection - std::list<int>::iterator list_it = index_list.begin(); - while (list_it != index_list.end()) + FT epsilon_max = 0; //sampling_radius(t,0); + FT sum_epsilon = 0; + int count_simplices = 0; + std::vector<int> histo(n+1, 0); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) { - if (!is_violating_protection(W[*list_it], t, D, delta)) - { - // If no conflicts then insert in every copy of T^3 - is_violating_protection(W[*list_it], t, D, delta); - insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - index_list.erase(list_it); - list_it = index_list.begin(); - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - list_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //write_delaunay_mesh(t, W[*list_it]); + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + FT r = sqrt(Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin()))); + if (r > epsilon_max) + epsilon_max = r; + sum_epsilon += r; + count_simplices++; + histo[floor(r/epsilon_max*n)]++; } - fill_landmarks(W, landmarks, landmarks_ind); - fill_full_cell_vector(t, full_cells); - if (triangulation_is_protected(t, delta)) - std::cout << "Triangulation is ok\n"; - else - std::cout << "Triangulation is BAD!! T_T ã—ãã—ã!\n"; - write_delaunay_mesh(t, Point_d(std::vector<FT>({0,0}))); - //std::cout << t << std::endl; + std::ofstream ofs ("histogram.tikz", std::ofstream::out); + FT barwidth = 20.0/n; + int max_value = *(std::max_element(histo.begin(), histo.end())); + std::cout << max_value << std::endl; + FT ten_power = pow(10, ceil(log10(max_value))); + FT max_histo = ten_power; + if (max_value/ten_power < 2) + max_histo = 0.2*ten_power; + if (max_value/ten_power < 5) + max_histo = 0.5*ten_power; + std::cout << ceil(log10(max_value)) << std::endl << max_histo << std::endl; + FT unitht = max_histo/10.0; + + ofs << "\\draw[->] (0,0) -- (0,11);\n" << + "\\draw[->] (0,0) -- (21,0);\n" << + "\\foreach \\i in {1,...,10}\n" << + "\\draw (0,\\i) -- (-0.1,\\i);\n" << + "\\foreach \\i in {1,...,20}\n" << + "\\draw (\\i,0) -- (\\i,-0.1);\n" << + + "\\node at (-1,11) {$\\epsilon$};\n" << + "\\node at (22,-1) {$\\epsilon/\\epsilon_{max}$};\n" << + "\\node at (-0.5,-0.5) {0};\n" << + "\\node at (-0.5,10) {" << max_histo << "};\n" << + "\\node at (20,-0.5) {1};\n"; + + + for (int i = 0; i < n; ++i) + ofs << "\\draw (" << barwidth*i << "," << histo[i]/unitht << ") -- (" + << barwidth*(i+1) << "," << histo[i]/unitht << ") -- (" + << barwidth*(i+1) << ",0) -- (" << barwidth*i << ",0) -- cycle;\n"; + + ofs.close(); + + //return sum_epsilon/count_simplices; + return epsilon_max; } -template <typename T> -void print_vector(std::vector<T> v) +FT epsilon_histogram_by_delaunay(Point_Vector& points, int n) { - std::cout << "["; - if (!v.empty()) - { - std::cout << *(v.begin()); - for (auto it = v.begin()+1; it != v.end(); ++it) - { - std::cout << ","; - std::cout << *it; - } - } - std::cout << "]"; + Delaunay_triangulation t(points[0].size()); + t.insert(points.begin(), points.end()); + return epsilon_histogram(t, n); } + int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind, std::vector<std::vector<int>>& full_cells) { //******************** Preface: origin point @@ -764,7 +378,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std } } std::string out_file = "wl_result"; - write_wl(out_file,WL); + //write_wl(out_file,WL); //******************** Constructng a witness complex std::cout << "Entered witness complex construction\n"; @@ -787,7 +401,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std not_in << " are not.\n"; //******************** Making a set of bad link landmarks - /* + std::cout << "Entered bad links\n"; std::set< int > perturbL; int count_badlinks = 0; @@ -814,7 +428,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std if (count_bad[i] != 0) std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; - */ + //*********************** Perturb bad link landmarks /* for (auto u: perturbL) @@ -848,16 +462,19 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std ofs.close(); } - write_edges("landmarks/edges", witnessComplex, landmarks); + //write_edges("landmarks/edges", witnessComplex, landmarks); /* return count_badlinks; */ return 0; } - int main (int argc, char * const argv[]) { + power_protection = true;//false; + grid_points = false;//true; + torus = true; + if (argc != 4) { std::cerr << "Usage: " << argv[0] @@ -866,40 +483,98 @@ int main (int argc, char * const argv[]) } int nbP = atoi(argv[1]); int dim = atoi(argv[2]); - double delta = atof(argv[3]); + double theta0 = atof(argv[3]); + //double delta = atof(argv[3]); + + is2d = (dim == 2); std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; - generate_points_random_box(point_vector, nbP, dim); - FT epsilon = point_sampling_radius_by_delaunay(point_vector); + if (grid_points) + { + generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim, torus); + nbP = (int)pow((int)pow(nbP, 1.0/dim), dim); + } + else + generate_points_random_box(point_vector, nbP, dim); + FT epsilon = point_sampling_radius_by_delaunay(point_vector, 0); + //FT epsilon = epsilon_histogram_by_delaunay(point_vector,50); std::cout << "Initial epsilon = " << epsilon << std::endl; Point_Vector L; std::vector<int> chosen_landmarks; //write_points("landmarks/initial_pointset",point_vector); //write_points("landmarks/initial_landmarks",L); CGAL::Timer timer; + + int n = 1; + std::vector<FT> values(n,0); + std::vector<FT> time(n,0); + + //FT step = 0.001; + //FT delta = 0.01*epsilon; + //FT alpha = 0.5; + //FT step = atof(argv[3]); + + start_experiments(point_vector, theta0, chosen_landmarks, epsilon); + + // for (int i = 0; i < n; i++) + // //for (int i = 0; bl > 0; i++) + // { + // //std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + // //double delta = pow(10, -(1.0*i)/2); + // //delta = step*i*epsilon; + // //theta0 = step*i; + // std::cout << "delta/epsilon = " << delta/epsilon << std::endl; + // std::cout << "theta0 = " << theta0 << std::endl; + // // Averaging the result + // int sum_values = 0; + // int nb_iterations = 1; + // std::vector<std::vector<int>> full_cells; + // for (int i = 0; i < nb_iterations; ++i) + // { + // //L = {}; + // chosen_landmarks = {}; + // //full_cells = {}; + // //timer.start(); + // //protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, epsilon, alpha, theta0, full_cells, torus, power_protection); + // protected_delaunay(point_vector, chosen_landmarks, delta, epsilon, alpha, theta0, torus, power_protection); + // //timer.stop(); + // sum_values += chosen_landmarks.size(); + // } + // //FT epsilon2 = point_sampling_radius_by_delaunay(L, epsilon); + // //std::cout << "Final epsilon = " << epsilon2 << ". Ratio = " << epsilon2/epsilon << std::endl; + // //write_points("landmarks/initial_landmarks",L); + // //std::cout << "delta/epsilon' = " << delta/epsilon2 << std::endl; + // FT nbL = (sum_values*1.0)/nb_iterations; + // //values[i] = pow((1.0*nbL)/nbP, -1.0/dim); + // values[i] = (1.0*nbL)/nbP; + // std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl; + // //landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); + // time[i] = timer.time(); + // timer.reset(); + // //write_points("landmarks/landmarks0",L); + // } + + // // OUTPUT A PLOT + // FT hstep = 20.0/(n-1); + // FT wstep = 10.0; + + // std::ofstream ofs("N'Nplot.tikz", std::ofstream::out); + // ofs << "\\draw[red] (0," << wstep*values[0] << ")"; + // for (int i = 1; i < n; ++i) + // ofs << " -- (" << hstep*i << "," << wstep*values[i] << ")"; + // ofs << ";\n"; + // ofs.close(); /* - for (int i = 0; i < 11; i++) - //for (int i = 0; bl > 0; i++) - { - //std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; - double delta = pow(10, -(1.0*i)/2); - std::cout << "delta = " << delta << std::endl; - L = {}; chosen_landmarks = {}; - std::vector<std::vector<int>> full_cells; - timer.start(); - landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, full_cells); - timer.stop(); - FT epsilon2 = point_sampling_radius_by_delaunay(L); - std::cout << "Final epsilon = " << epsilon2 << ". Ratio = " << epsilon/epsilon2 << std::endl; - write_points("landmarks/initial_landmarks",L); - int nbL = chosen_landmarks.size(); - std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl; - landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); - timer.reset(); - //write_points("landmarks/landmarks0",L); - } - */ + wstep = 0.1; + ofs = std::ofstream("time.tikz", std::ofstream::out); + ofs << "\\draw[red] (0," << wstep*time[0] << ")"; + for (int i = 1; i < n; ++i) + ofs << " -- (" << hstep*i << "," << wstep*time[i] << ")"; + ofs << ";\n"; + ofs.close(); + + std::vector<std::vector<int>> full_cells; timer.start(); landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, full_cells); @@ -909,6 +584,7 @@ int main (int argc, char * const argv[]) write_points("landmarks/initial_landmarks",L); int nbL = chosen_landmarks.size(); std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl; - landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); + //landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); timer.reset(); + */ } diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index 69ef5fbf..49383154 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -776,8 +776,8 @@ int main (int argc, char * const argv[]) std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; //read_points_cust(file_name, point_vector); - generate_points_random_box(point_vector, nbP, dim); - //generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim); + //generate_points_random_box(point_vector, nbP, dim); + generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim); //nbP = (int)(pow((int)pow(nbP, 1.0/dim), dim)); /* for (auto &p: point_vector) diff --git a/src/Witness_complex/example/witness_complex_knn_landmarks.cpp b/src/Witness_complex/example/witness_complex_knn_landmarks.cpp index e4a1c324..c45bc0c1 100644 --- a/src/Witness_complex/example/witness_complex_knn_landmarks.cpp +++ b/src/Witness_complex/example/witness_complex_knn_landmarks.cpp @@ -32,6 +32,8 @@ //#include "gudhi/graph_simplicial_complex.h" #include "gudhi/Witness_complex.h" #include "gudhi/reader_utils.h" +#include "generators.h" +#include "output.h" //#include <boost/filesystem.hpp> //#include <CGAL/Delaunay_triangulation.h> @@ -73,60 +75,6 @@ typedef K_neighbor_search::iterator KNS_range; typedef boost::container::flat_map<int, int> Point_etiquette_map; typedef std::vector<Point_d> Point_Vector; -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , Point_Vector & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - Point_d p(point.begin(), point.end()); - if (point.size() != 1) - points.push_back(p); - } - in_file.close(); -} - -/* -void read_points_to_tree (std::string file_name, Tree& tree) -{ - //I assume here that tree is empty - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector<double> coords; - std::istringstream iss( line ); - while(iss >> x) { coords.push_back(x); } - if (coords.size() != 1) - { - Point_d point(coords.begin(), coords.end()); - tree.insert(point); - } - } - in_file.close(); -} -*/ /** Function that chooses landmarks from W and place it in the kd-tree L. * Note: nbL hould be removed if the code moves to Witness_complex @@ -184,19 +132,6 @@ void d_nearest_landmarks(Point_Vector &W, Tree &L, Point_etiquette_map &L_i, std } } - -void write_wl( std::string file_name, std::vector< std::vector <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - int main (int argc, char * const argv[]) { if (argc != 3) @@ -270,6 +205,6 @@ int main (int argc, char * const argv[]) out_file = "output/"+file_name+"_"+argv[2]+".badlinks"; std::ofstream ofs2(out_file, std::ofstream::out); - witnessComplex.write_bad_links(ofs2); + //witnessComplex.write_bad_links(ofs2); ofs2.close(); } diff --git a/src/Witness_complex/example/witness_complex_protected_delaunay.cpp b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp index 2f795a5f..77a167a5 100644 --- a/src/Witness_complex/example/witness_complex_protected_delaunay.cpp +++ b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp @@ -268,6 +268,15 @@ void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, landmark_count++; } + + + +//////////////////////////////////////////////////////////////////////// +// OLD CODE VVVVVVVV +//////////////////////////////////////////////////////////////////////// + + +/* bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) { Euclidean_distance ed; @@ -592,3 +601,4 @@ int main (int argc, char * const argv[]) } } +*/ diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index 550c9392..bf3015fa 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -35,6 +35,8 @@ //#include "gudhi/graph_simplicial_complex.h" #include "gudhi/Witness_complex.h" #include "gudhi/reader_utils.h" +#include "generators.h" +#include "output.h" //#include <boost/filesystem.hpp> //#include <CGAL/Delaunay_triangulation.h> @@ -94,101 +96,9 @@ typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; typedef std::vector<Point_d> Point_Vector; //typedef K::Equal_d Equal_d; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; bool toric=false; -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , Point_Vector & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - Point_d p(point.begin(), point.end()); - if (point.size() != 1) - points.push_back(p); - } - in_file.close(); -} - -void generate_points_grid(Point_Vector& W, int width, int D) -{ - -} - -void generate_points_random_box(Point_Vector& W, int nbP, int dim) -{ - Random_cube_iterator rp(dim, 1); - for (int i = 0; i < nbP; i++) - { - W.push_back(*rp++); - } -} - -/* NOT TORUS RELATED - */ -void generate_points_sphere(Point_Vector& W, int nbP, int dim) -{ - CGAL::Random_points_on_sphere_d<Point_d> rp(dim,1); - for (int i = 0; i < nbP; i++) - W.push_back(*rp++); -} -/* -void read_points_to_tree (std::string file_name, Tree& tree) -{ - //I assume here that tree is empty - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector<double> coords; - std::istringstream iss( line ); - while(iss >> x) { coords.push_back(x); } - if (coords.size() != 1) - { - Point_d point(coords.begin(), coords.end()); - tree.insert(point); - } - } - in_file.close(); -} -*/ - -void write_wl( std::string file_name, std::vector< std::vector <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - - std::vector<Point_d> convert_to_torus(std::vector< Point_d>& points) { std::vector< Point_d > points_torus; @@ -205,82 +115,6 @@ std::vector<Point_d> convert_to_torus(std::vector< Point_d>& points) return points_torus; } -void write_points_torus( std::string file_name, std::vector< Point_d > & points) -{ - std::ofstream ofs (file_name, std::ofstream::out); - std::vector<Point_d> points_torus = convert_to_torus(points); - for (auto w : points_torus) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); -} - -void write_points( std::string file_name, std::vector< Point_d > & points) -{ - if (toric) write_points_torus(file_name, points); - else - { - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : points) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); - } -} - - -void write_edges_torus(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - Point_Vector l_torus = convert_to_torus(landmarks); - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = l_torus[u].cartesian_begin(); it != l_torus[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = l_torus[v].cartesian_begin(); it != l_torus[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); -} - -void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - if (toric) write_edges_torus(file_name, witness_complex, landmarks); - else - { - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); - } -} - - /** Function that chooses landmarks from W and place it in the kd-tree L. * Note: nbL hould be removed if the code moves to Witness_complex */ @@ -356,6 +190,7 @@ void landmark_choice_600cell(Point_Vector&W, int nbP, int nbL, Point_Vector& lan int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector<int>& landmarks_ind) { //********************Preface: origin point + clock_t start, end; int D = W[0].size(); std::vector<FT> orig_vector; for (int i=0; i<D; i++) @@ -383,6 +218,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< */ std::cout << "Enter (D+1) nearest landmarks\n"; //std::cout << "Size of the tree is " << L.size() << std::endl; + start = clock(); for (int i = 0; i < nbP; i++) { //std::cout << "Entered witness number " << i << std::endl; @@ -416,7 +252,9 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< } } //std::cout << "\n"; - + end = clock(); + std::cout << "Landmark choice for " << nbL << " landmarks took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; std::string out_file = "wl_result"; write_wl(out_file,WL); @@ -424,14 +262,19 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::cout << "Entered witness complex construction\n"; Witness_complex<> witnessComplex; witnessComplex.setNbL(nbL); + start = clock(); witnessComplex.witness_complex(WL); + // + end = clock(); + std::cout << "Howdy world! The process took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + //witnessComplex.witness_complex(WL); /* if (witnessComplex.is_witness_complex(WL)) std::cout << "!!YES. IT IS A WITNESS COMPLEX!!\n"; else std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n"; */ - */ //******************** Making a set of bad link landmarks std::cout << "Entered bad links\n"; std::set< int > perturbL; @@ -575,8 +418,8 @@ int main (int argc, char * const argv[]) */ } int bl = nbL, curr_min = bl; - write_points("landmarks/initial_pointset",point_vector); - write_points("landmarks/initial_landmarks",L); + //write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_landmarks",L); for (int i = 0; bl > 0; i++) //for (int i = 0; i < 1; i++) @@ -585,7 +428,7 @@ int main (int argc, char * const argv[]) bl=landmark_perturbation(point_vector, L, chosen_landmarks); if (bl < curr_min) curr_min=bl; - write_points("landmarks/landmarks0",L); + //write_points("landmarks/landmarks0",L); } //end = clock(); diff --git a/src/cmake/modules/FindQGLViewer.cmake b/src/cmake/modules/FindQGLViewer.cmake new file mode 100644 index 00000000..65723d67 --- /dev/null +++ b/src/cmake/modules/FindQGLViewer.cmake @@ -0,0 +1,61 @@ +# - Try to find QGLViewer +# Once done this will define +# +# QGLVIEWER_FOUND - system has QGLViewer +# QGLVIEWER_INCLUDE_DIR - the QGLViewer include directory +# QGLVIEWER_LIBRARIES - Link these to use QGLViewer +# QGLVIEWER_DEFINITIONS - Compiler switches required for using QGLViewer +# + +find_path(QGLVIEWER_INCLUDE_DIR + NAMES QGLViewer/qglviewer.h + PATHS /usr/include + /usr/local/include + ENV QGLVIEWERROOT + ) + +find_library(QGLVIEWER_LIBRARY_RELEASE + NAMES qglviewer-qt4 qglviewer QGLViewer QGLViewer2 + PATHS /usr/lib + /usr/local/lib + ENV QGLVIEWERROOT + ENV LD_LIBRARY_PATH + ENV LIBRARY_PATH + PATH_SUFFIXES QGLViewer QGLViewer/release + ) + +find_library(QGLVIEWER_LIBRARY_DEBUG + NAMES dqglviewer dQGLViewer dQGLViewer2 QGLViewerd2 + PATHS /usr/lib + /usr/local/lib + ENV QGLVIEWERROOT + ENV LD_LIBRARY_PATH + ENV LIBRARY_PATH + PATH_SUFFIXES QGLViewer QGLViewer/debug + ) + +if(QGLVIEWER_LIBRARY_RELEASE) + if(QGLVIEWER_LIBRARY_DEBUG) + set(QGLVIEWER_LIBRARIES_ optimized ${QGLVIEWER_LIBRARY_RELEASE} debug ${QGLVIEWER_LIBRARY_DEBUG}) + else() + set(QGLVIEWER_LIBRARIES_ ${QGLVIEWER_LIBRARY_RELEASE}) + endif() + + set(QGLVIEWER_LIBRARIES ${QGLVIEWER_LIBRARIES_} CACHE FILEPATH "The QGLViewer library") + +endif() + +IF(QGLVIEWER_INCLUDE_DIR AND QGLVIEWER_LIBRARIES) + SET(QGLVIEWER_FOUND TRUE) +ENDIF(QGLVIEWER_INCLUDE_DIR AND QGLVIEWER_LIBRARIES) + +IF(QGLVIEWER_FOUND) + IF(NOT QGLViewer_FIND_QUIETLY) + MESSAGE(STATUS "Found QGLViewer: ${QGLVIEWER_LIBRARIES}") + ENDIF(NOT QGLViewer_FIND_QUIETLY) +ELSE(QGLVIEWER_FOUND) + IF(QGLViewer_FIND_REQUIRED) + MESSAGE(FATAL_ERROR "Could not find QGLViewer") + ENDIF(QGLViewer_FIND_REQUIRED) +ENDIF(QGLVIEWER_FOUND) + diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 83908905..315aa0ac 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -6,7 +6,8 @@ The Gudhi library (Geometric Understanding in Higher Dimensions) is a generic C++ library for topological analysis of high-dimensional data whose goal is to provide robust, efficient, flexible and easy to use implementations of -state-of-the-art algorithms and data structures for computational topology. +state-of-the-art algorithms and data structures for computational topology. +This library is part of the <a href="https://project.inria.fr/gudhi/">Gudhi project</a>. The current release of the library allows to use several data-structures for simplicial complexes : simplex tree, Hasse diagram or skeleton-blocker. Several operations can then be done on top of these @@ -17,48 +18,52 @@ We refer to \cite gudhilibrary_ICMS14 for a detailed description of the design of the library. +\section installation Gudhi installation + +As Gudhi is a header only library, there is no need to install the library. + +Examples of Gudhi headers inclusion can be found in \ref demos. + \section compiling Compiling -The library uses c++11 and requires Boost with version 1.48.0 or more recent : http://www.boost.org/. -It is a multiplaform library and compiles on Linux, Mac OSX and Visual Studio 2013. +The library uses c++11 and requires <a href="http://www.boost.org/">Boost</a> with version 1.48.0 or more recent. +It is a multi-platform library and compiles on Linux, Mac OSX and Visual Studio 2013. \subsection gmp GMP: The multi-field persistent homology algorithm requires GMP which is a free library for arbitrary-precision -arithmetic, operating on signed integers, rational numbers, and floating point numbers -The following examples require The GNU Multiple Precision Arithmetic Library (GMP) http://gmplib.org/ +arithmetic, operating on signed integers, rational numbers, and floating point numbers. + +The following examples require the <a href="http://gmplib.org/">GNU Multiple Precision Arithmetic Library</a> (GMP) and will not be built if GMP is not installed: - Persistent_cohomology/rips_multifield_persistence - Simplex_tree/simplex_tree_from_alpha_shapes_3 -Having GMP version 4.2 or higher installed is recommended. This library can be obtained from http://gmplib.org/ +Having GMP version 4.2 or higher installed is recommended. \subsection cgal CGAL: CGAL is a C++ library which provides easy access to efficient and reliable geometric algorithms. -The following example requires CGAL https://www.cgal.org/ and will not be built if CGAL is not installed: +The following example requires the <a href="http://www.cgal.org/">Computational Geometry Algorithms Library</a> (CGAL) +and will not be built if CGAL is not installed: - Simplex_tree/simplex_tree_from_alpha_shapes_3 Having CGAL version 4.4 or higher installed is recommended. The procedure to install this library according to your operating system is detailed here http://doc.cgal.org/latest/Manual/installation.html -\section demos Demos and Examples +\subsection demos Demos and examples -To build the library, run the following in a terminal: +To build the demos and libraries, run the following commands in a terminal: \verbatim cd /path-to-gudhi/ mkdir build cd build/ -cmake -DCMAKE_BUILD_TYPE=Release .. +cmake .. make \endverbatim - - - - \details \copyright GNU General Public License v3. |