diff options
17 files changed, 1529 insertions, 14 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index edf7e63f..29ab2a55 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,7 +18,7 @@ if(MSVC) # Turn off some VC++ warnings SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4668 /wd4311 /wd4800 /wd4820 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018") else() - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wsign-compare") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}") endif() @@ -27,10 +27,10 @@ set(Boost_USE_STATIC_LIBS ON) set(Boost_USE_MULTITHREADED ON) set(Boost_USE_STATIC_RUNTIME OFF) -find_package(GMP) -if(GMP_FOUND) - find_package(GMPXX) -endif() +#find_package(GMP) +#if(GMP_FOUND) + #find_package(GMPXX) +#endif() find_package(CGAL) @@ -39,16 +39,28 @@ set(TBB_FIND_QUIETLY ON) find_package(TBB) # Required programs for unitary tests purpose -FIND_PROGRAM( GCOVR_PATH gcovr ) -if (GCOVR_PATH) - message("gcovr found in ${GCOVR_PATH}") +FIND_PROGRAM( LCOV_PATH lcov ) +if (LCOV_PATH) + message("lcov found in ${LCOV_PATH}") endif() -# Required programs for unitary tests purpose -FIND_PROGRAM( GPROF_PATH gprof ) -if (GPROF_PATH) - message("gprof found in ${GPROF_PATH}") + +FIND_PROGRAM( PYTHON_PATH python ) +if (PYTHON_PATH) + message("python found in ${PYTHON_PATH}") endif() +# Function to add_test cpplint on each header file of the Gudhi module +function(cpplint_add_tests the_directory) + if (PYTHON_PATH) + # Cpplint tests on coding style + file(GLOB files "${the_directory}/*.h" "${the_directory}/*/*.h") + foreach(filename ${files}) + message(${filename}) + add_test("${filename}.cpplint" ${CMAKE_SOURCE_DIR}/scripts/check_google_style.sh ${filename} ${CMAKE_SOURCE_DIR}/scripts/cpplint.py) + endforeach() + endif() +endfunction(cpplint_add_tests) + if(NOT Boost_FOUND) message(FATAL_ERROR "NOTICE: This demo requires Boost and will not be compiled.") @@ -71,6 +83,7 @@ else() include_directories(src/Persistent_cohomology/include/) include_directories(src/Simplex_tree/include/) include_directories(src/Skeleton_blocker/include/) + include_directories(src/Witness_complex/include/) add_subdirectory(src/Simplex_tree/test) add_subdirectory(src/Simplex_tree/example) @@ -79,6 +92,13 @@ else() add_subdirectory(src/Skeleton_blocker/test) add_subdirectory(src/Skeleton_blocker/example) add_subdirectory(src/Contraction/example) + # add_subdirectory(src/Hasse_complex/example) + add_subdirectory(src/Witness_complex/test) + add_subdirectory(src/Witness_complex/example) + # add_subdirectory(src/Alpha_shapes/example) + # add_subdirectory(src/Alpha_shapes/test) + add_subdirectory(src/Bottleneck/example) + add_subdirectory(src/Bottleneck/test) # data points generator add_subdirectory(data/points/generator) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3c38483c..5e79e0b5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -52,6 +52,10 @@ else() add_subdirectory(example/Persistent_cohomology) add_subdirectory(example/Skeleton_blocker) add_subdirectory(example/Contraction) + add_subdirectory(example/Hasse_complex) + add_subdirectort(example/Witness_complex) + # add_subdirectory(example/Alpha_shapes) + add_subdirectory(example/Bottleneck) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Doxyfile b/src/Doxyfile index faa0d3fe..6585e50c 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -834,8 +834,8 @@ EXAMPLE_RECURSIVE = NO IMAGE_PATH = doc/Skeleton_blocker/ \ doc/common/ \ - doc/Contraction/ - + doc/Contraction/ \ + doc/Witness_complex/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/Witness_complex/doc/bench_Cy8.png b/src/Witness_complex/doc/bench_Cy8.png Binary files differnew file mode 100644 index 00000000..d9045294 --- /dev/null +++ b/src/Witness_complex/doc/bench_Cy8.png diff --git a/src/Witness_complex/doc/bench_sphere.png b/src/Witness_complex/doc/bench_sphere.png Binary files differnew file mode 100644 index 00000000..ba6bb381 --- /dev/null +++ b/src/Witness_complex/doc/bench_sphere.png diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt new file mode 100644 index 00000000..b5770ba6 --- /dev/null +++ b/src/Witness_complex/example/CMakeLists.txt @@ -0,0 +1,29 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIWitnessComplex) + +# A simple example + add_executable( witness_complex_from_file witness_complex_from_file.cpp ) + add_test( witness_complex_from_bunny &{CMAKE_CURRENT_BINARY_DIR}/witness_complex_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + 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_sphere witness_complex_sphere.cpp ) + target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for witness_complex_sphere example.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile witness_complex_sphere example. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Witness_complex/example/generators.h b/src/Witness_complex/example/generators.h new file mode 100644 index 00000000..0d42cda2 --- /dev/null +++ b/src/Witness_complex/example/generators.h @@ -0,0 +1,134 @@ +#ifndef GENERATORS_H +#define GENERATORS_H + +#include <fstream> + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::FT FT; +typedef K::Point_d Point_d; +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; + +/** + * \brief Rock age method of reading off file + * + */ +inline void +off_reader_cust ( std::string file_name , std::vector<Point_d> & 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; + // Line OFF. No need in it + if (!getline(in_file, line)) + { + std::cerr << "No line OFF\n"; + return; + } + // Line with 3 numbers. No need + if (!getline(in_file, line)) + { + std::cerr << "No line with 3 numbers\n"; + return; + } + // Reading points + while( getline ( in_file , line ) ) + { + std::vector< double > point; + std::istringstream iss( line ); + while(iss >> x) { point.push_back(x); } + points.push_back(Point_d(point)); + } + in_file.close(); +} + +/** + * \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(); +} + +/** \brief Generate points on a grid in a cube of side 2 + * having {+-1}^D as vertices and insert them in W. + * The grid has "width" points on each side. + * If torus is true then it is supposed that the cube represents + * a flat torus, hence the opposite borders are associated. + * The points on border in this case are not placed twice. + */ +void generate_points_grid(Point_Vector& W, int width, int D, bool torus) +{ + 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) + { + if (torus) + point.push_back(-1+(2.0/(width-1))*(cell_i%width)); + else + point.push_back(-1+(2.0/width)*(cell_i%width)); + //attention: the bottom and the right are covered too! + cell_i /= width; + } + W.push_back(point); + } +} + +/** \brief Generate nbP points uniformly in a cube of side 2 + * having {+-1}^dim as its vertices and insert them in W. + */ +void generate_points_random_box(Point_Vector& W, int nbP, int dim) +{ + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + +/** \brief Generate nbP points uniformly on a (dim-1)-sphere + * and insert them in W. + */ +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++); +} + + +#endif diff --git a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp new file mode 100644 index 00000000..067321ce --- /dev/null +++ b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp @@ -0,0 +1,461 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <fstream> +#include <ctime> +#include <utility> +#include <algorithm> +#include <set> +#include <queue> +#include <iterator> + +#include <sys/types.h> +#include <sys/stat.h> +//#include <stdlib.h> + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Relaxed_witness_complex.h" +#include "gudhi/reader_utils.h" +#include "gudhi/Collapse/Collapse.h" +//#include <boost/filesystem.hpp> + +//#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Cartesian_d.h> +#include <CGAL/Search_traits.h> +#include <CGAL/Search_traits_adapter.h> +#include <CGAL/property_map.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Orthogonal_incremental_neighbor_search.h> +#include <CGAL/Kd_tree.h> +#include <CGAL/Euclidean_distance.h> + +#include <CGAL/Kernel_d/Vector_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/constructions_d.h> +#include <CGAL/Fuzzy_sphere.h> +#include <CGAL/Random.h> + + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::FT FT; +typedef K::Point_d Point_d; +typedef CGAL::Search_traits< + FT, Point_d, + typename K::Cartesian_const_iterator_d, + typename K::Construct_cartesian_const_iterator_d> Traits_base; +typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; +//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + +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_incremental_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> Neighbor_search; +typedef Neighbor_search::Tree Tree; +typedef Neighbor_search::Distance Distance; +typedef Neighbor_search::iterator KNS_iterator; +typedef Neighbor_search::iterator KNS_range; +typedef boost::container::flat_map<int, int> Point_etiquette_map; +typedef CGAL::Kd_tree<STraits> Tree2; + +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_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 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(); +} + +void write_rl( std::string file_name, std::vector< std::vector <std::vector<int>::iterator> > & rl) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : rl) + { + 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; + for (auto p: points) + { + FT theta = M_PI*p[0]; + FT phi = M_PI*p[1]; + std::vector<FT> p_torus; + p_torus.push_back((1+0.2*cos(theta))*cos(phi)); + p_torus.push_back((1+0.2*cos(theta))*sin(phi)); + p_torus.push_back(0.2*sin(theta)); + points_torus.push_back(Point_d(p_torus)); + } + 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 + */ +void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + //std::vector<Point_d> landmarks; + int chosen_landmark; + //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false); + Point_d* p; + CGAL::Random rand; + for (int i = 0; i < nbL; i++) + { + // while (!res.second) + // { + do chosen_landmark = rand.get_int(0,nbP); + while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end()); + //rand++; + //std::cout << "Chose " << chosen_landmark << std::endl; + p = &W[chosen_landmark]; + //L_i.emplace(chosen_landmark,i); + // } + landmarks.push_back(*p); + landmarks_ind.push_back(chosen_landmark); + //std::cout << "Added landmark " << chosen_landmark << std::endl; + } + } + + +void landmarks_to_witness_complex(Point_Vector &W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT alpha) +{ + //********************Preface: origin point + unsigned D = W[0].size(); + std::vector<FT> orig_vector; + for (unsigned i = 0; i < D; i++) + orig_vector.push_back(0); + Point_d origin(orig_vector); + //Distance dist; + //dist.transformed_distance(0,1); + //******************** Constructing a WL matrix + int nbP = W.size(); + int nbL = landmarks.size(); + STraits traits(&(landmarks[0])); + Euclidean_distance ed; + std::vector< std::vector <int> > WL(nbP); + std::vector< std::vector< typename std::vector<int>::iterator > > ope_limits(nbP); + Tree L(boost::counting_iterator<std::ptrdiff_t>(0), + boost::counting_iterator<std::ptrdiff_t>(nbL), + typename Tree::Splitter(), + traits); + + std::cout << "Enter (D+1) nearest landmarks\n"; + //std::cout << "Size of the tree is " << L.size() << std::endl; + for (int i = 0; i < nbP; i++) + { + //std::cout << "Entered witness number " << i << std::endl; + Point_d& w = W[i]; + std::queue< typename std::vector<int>::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark + Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0]))); + Neighbor_search::iterator search_it = search.begin(); + + //Incremental search and filling WL + while (WL[i].size() < D) + WL[i].push_back((search_it++)->first); + FT dtow = ed.transformed_distance(w, landmarks[WL[i][D-1]]); + while (search_it->second < dtow + alpha) + WL[i].push_back((search_it++)->first); + + //Filling the (1+epsilon)-limits table + for (std::vector<int>::iterator wl_it = WL[i].begin(); wl_it != WL[i].end(); ++wl_it) + { + ope_queue.push(wl_it); + FT d_to_curr_l = ed.transformed_distance(w, landmarks[*wl_it]); + //std::cout << "d_to_curr_l=" << d_to_curr_l << std::endl; + //std::cout << "d_to_front+alpha=" << d_to_curr_l << std::endl; + while (d_to_curr_l > alpha + ed.transformed_distance(w, landmarks[*(ope_queue.front())])) + { + ope_limits[i].push_back(wl_it); + ope_queue.pop(); + } + } + while (ope_queue.size() > 0) + { + ope_limits[i].push_back(WL[i].end()); + ope_queue.pop(); + } + //std::cout << "Safely constructed a point\n"; + ////Search D+1 nearest neighbours from the tree of landmarks L + /* + if (w[0]>0.95) + std::cout << i << std::endl; + */ + //K_neighbor_search search(L, w, D, FT(0), true, + // CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])) ); + //std::cout << "Safely found nearest landmarks\n"; + /* + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) + { + //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; + //Point_etiquette_map::iterator itm = L_i.find(it->first); + //assert(itm != L_i.end()); + //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; + WL[i].push_back(it->first); + //std::cout << "ITFIRST " << it->first << std::endl; + //std::cout << i << " " << it->first << ": " << it->second << std::endl; + } + */ + } + //std::cout << "\n"; + + //std::string out_file = "wl_result"; + write_wl("wl_result",WL); + write_rl("rl_result",ope_limits); + + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.relaxed_witness_complex(WL, ope_limits); + char buffer[100]; + int i = sprintf(buffer,"stree_result.txt"); + + if (i >= 0) + { + std::string out_file = (std::string)buffer; + std::ofstream ofs (out_file, std::ofstream::out); + witnessComplex.st_to_file(ofs); + ofs.close(); + } + write_edges("landmarks/edges", witnessComplex, landmarks); + std::cout << Distance().transformed_distance(Point_d(std::vector<double>({0.1,0.1})), Point_d(std::vector<double>({1.9,1.9}))) << std::endl; +} + + +int main (int argc, char * const argv[]) +{ + + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] + << " nbP nbL dim alpha\n"; + return 0; + } + /* + boost::filesystem::path p; + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + + int nbP = atoi(argv[1]); + int nbL = atoi(argv[2]); + int dim = atoi(argv[3]); + double alpha = atof(argv[4]); + //clock_t start, end; + //Construct the Simplex Tree + Witness_complex<> witnessComplex; + + std::cout << "Let the carnage begin!\n"; + Point_Vector point_vector; + //read_points_cust(file_name, point_vector); + generate_points_sphere(point_vector, nbP, dim); + /* + for (auto &p: point_vector) + { + assert(std::count(point_vector.begin(),point_vector.end(),p) == 1); + } + */ + //std::cout << "Successfully read the points\n"; + //witnessComplex.setNbL(nbL); + Point_Vector L; + std::vector<int> chosen_landmarks; + landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); + //start = clock(); + + write_points("landmarks/initial_pointset",point_vector); + write_points("landmarks/initial_landmarks",L); + + landmarks_to_witness_complex(point_vector, L, chosen_landmarks, alpha); + //end = clock(); + + /* + std::cout << "Landmark choice took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + start = clock(); + witnessComplex.witness_complex(WL); + // + end = clock(); + std::cout << "Howdy world! The process took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + */ + + /* + out_file = "output/"+file_name+"_"+argv[2]+".stree"; + std::ofstream ofs (out_file, std::ofstream::out); + witnessComplex.st_to_file(ofs); + ofs.close(); + + out_file = "output/"+file_name+"_"+argv[2]+".badlinks"; + std::ofstream ofs2(out_file, std::ofstream::out); + witnessComplex.write_bad_links(ofs2); + ofs2.close(); + */ +} diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp new file mode 100644 index 00000000..6add4e0a --- /dev/null +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -0,0 +1,118 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <fstream> +#include <ctime> + +#include <sys/types.h> +#include <sys/stat.h> + +#include "gudhi/Simplex_tree.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/Landmark_choice_by_random_point.h" +#include "gudhi/reader_utils.h" + +using namespace Gudhi; + +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< std::vector <double> > Point_Vector; + +typedef Witness_complex< Simplex_tree<> > WitnessComplex; + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , std::vector< std::vector< double > > & 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); } + if (point.size() != 1) + points.push_back(point); + } + in_file.close(); +} + +/** Write a gnuplot readable file. + * Data range is a random access range of pairs (arg, value) + */ +template < typename Data_range > +void write_data( Data_range & data, std::string filename ) +{ + std::ofstream ofs(filename, std::ofstream::out); + for (auto entry: data) + ofs << entry.first << ", " << entry.second << "\n"; + ofs.close(); +} + +int main (int argc, char * const argv[]) +{ + if (argc != 3) + { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file nbL \n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]); + clock_t start, end; + + // Construct the Simplex Tree + Simplex_tree<> simplex_tree; + + // Read the point file + Point_Vector point_vector; + read_points_cust(file_name, point_vector); + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + std::vector<std::vector< int > > knn; + Landmark_choice_by_random_point(point_vector, nbL, knn); + end = clock(); + std::cout << "Landmark choice for " << nbL << " landmarks took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + + // Compute witness complex + start = clock(); + WitnessComplex(knn, simplex_tree, nbL, point_vector[0].size()); + end = clock(); + std::cout << "Witness complex took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + +} diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp new file mode 100644 index 00000000..1fa0fc42 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -0,0 +1,126 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ +#define BOOST_PARAMETER_MAX_ARITY 12 + + +#include <iostream> +#include <fstream> +#include <ctime> +#include <utility> + +#include <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/reader_utils.h> + +#include "generators.h" + +using namespace Gudhi; + +typedef std::vector< Vertex_handle > typeVectorVertex; +//typedef std::vector< std::vector <double> > Point_Vector; + +typedef Witness_complex< Simplex_tree<> > WitnessComplex; + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , std::vector< std::vector< double > > & 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); } + if (point.size() != 1) + points.push_back(point); + } + in_file.close(); +} + +/** Write a gnuplot readable file. + * Data range is a random access range of pairs (arg, value) + */ +template < typename Data_range > +void write_data( Data_range & data, std::string filename ) +{ + std::ofstream ofs(filename, std::ofstream::out); + for (auto entry: data) + ofs << entry.first << ", " << entry.second << "\n"; + ofs.close(); +} + +int main (int argc, char * const argv[]) +{ + if (argc != 3) + { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file nbL \n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]); + clock_t start, end; + + // Construct the Simplex Tree + Simplex_tree<> simplex_tree; + + std::vector< std::pair<int, double> > l_time; + + // Read the point file + for (int nbP = 500; nbP < 10000; nbP += 500) + { + Point_Vector point_vector; + generate_points_sphere(point_vector, nbP, 4); + std::cout << "Successfully generated " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + std::vector<std::vector< int > > knn; + Landmark_choice_by_random_point(point_vector, nbL, knn); + + // Compute witness complex + WitnessComplex(knn, simplex_tree, nbL, point_vector[0].size()); + end = clock(); + double time = (double)(end-start)/CLOCKS_PER_SEC; + std::cout << "Witness complex for " << nbL << " landmarks took " + << time << " s. \n"; + l_time.push_back(std::make_pair(nbP,time)); + } + write_data(l_time, "w_time.dat"); +} diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h new file mode 100644 index 00000000..6ac59ae9 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -0,0 +1,99 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ +#define GUDHI_LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ + +/** + * \class Landmark_choice_by_furthest_point + * \brief The class `Landmark_choice_by_furthest_point` allows to construct the matrix + * of closest landmarks per witness by iteratively choosing the furthest witness + * from the set of already chosen landmarks as the new landmark. + * \ingroup witness_complex + */ + +class Landmark_choice_by_furthest_point { + +private: + typedef std::vector<int> typeVectorVertex; + +public: + +/** + * \brief Landmark choice strategy by iteratively adding the furthest witness from the + * current landmark set as the new landmark. + * \details It chooses nbL landmarks from a random access range `points` and + * writes {witness}*{closest landmarks} matrix in `knn`. + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + Landmark_choice_by_furthest_point(Point_random_access_range &points, + int nbL, + KNearestNeighbours &knn) + { + int nb_points = points.end() - points.begin(); + std::vector<std::vector<double>> wit_land_dist(nb_points, std::vector<double>()); // distance matrix witness x landmarks + typeVectorVertex chosen_landmarks; // landmark list + + knn = KNearestNeighbours(nb_points, std::vector<int>()); + int current_number_of_landmarks=0; // counter for landmarks + double curr_max_dist = 0; // used for defining the furhest point from L + const double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nb_points,infty); // vector of current distances to L from points + //int dim = points.begin()->size(); + + int rand_int = rand() % nb_points; + int curr_max_w = rand_int; //For testing purposes a pseudo-random number is used here + + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + { + //curr_max_w at this point is the next landmark + chosen_landmarks.push_back(curr_max_w); + unsigned i = 0; + for (auto& p: points) + { + double curr_dist = euclidean_distance(p, *(points.begin() + chosen_landmarks[current_number_of_landmarks])); + wit_land_dist[i].push_back(curr_dist); + knn[i].push_back(current_number_of_landmarks); + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + ++i; + } + curr_max_dist = 0; + for (i = 0; i < dist_to_L.size(); i++) + if (dist_to_L[i] > curr_max_dist) + { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + for (unsigned i = 0; i < points.size(); ++i) + std::sort(knn[i].begin(), + knn[i].end(), + [&wit_land_dist, i](int a, int b) + { return wit_land_dist[i][a] < wit_land_dist[i][b]; }); + } + +}; + +#endif diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h new file mode 100644 index 00000000..fa822591 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h @@ -0,0 +1,84 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_LANDMARK_CHOICE_BY_RANDOM_POINT_H_ +#define GUDHI_LANDMARK_CHOICE_BY_RANDOM_POINT_H_ + +/** + * \class Landmark_choice_by_random_point + * \brief The class `Landmark_choice_by_random_point` allows to construct the matrix + * of closest landmarks per witness by iteratively choosing a random non-chosen witness + * as a new landmark. + * \ingroup witness_complex + */ + +class Landmark_choice_by_random_point { + +public: + + /** \brief Landmark choice strategy by taking random vertices for landmarks. + * \details It chooses nbL distinct landmarks from a random access range `points` + * and outputs a matrix {witness}*{closest landmarks} in knn. + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + Landmark_choice_by_random_point(Point_random_access_range &points, int nbL, KNearestNeighbours &knn) + { + int nbP = points.end() - points.begin(); + std::set<int> landmarks; + int current_number_of_landmarks=0; // counter for landmarks + + int chosen_landmark = rand()%nbP; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + { + while (landmarks.find(chosen_landmark) != landmarks.end()) + chosen_landmark = rand()% nbP; + landmarks.insert(chosen_landmark); + } + + int dim = points.begin()->size(); + typedef std::pair<double,int> dist_i; + typedef bool (*comp)(dist_i,dist_i); + knn = KNearestNeighbours(nbP); + for (int points_i = 0; points_i < nbP; points_i++) + { + std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([&](dist_i j1, dist_i j2){return j1.first > j2.first;}); + std::set<int>::iterator landmarks_it; + int landmarks_i = 0; + for (landmarks_it = landmarks.begin(), landmarks_i=0; landmarks_it != landmarks.end(); landmarks_it++, landmarks_i++) + { + dist_i dist = std::make_pair(euclidean_distance(points[points_i],points[*landmarks_it]), landmarks_i); + l_heap.push(dist); + } + for (int i = 0; i < dim+1; i++) + { + dist_i dist = l_heap.top(); + knn[points_i].push_back(dist.second); + l_heap.pop(); + } + } + } + +}; + +#endif diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h new file mode 100644 index 00000000..915e445c --- /dev/null +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -0,0 +1,269 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_WITNESS_COMPLEX_H_ +#define GUDHI_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include <gudhi/distance_functions.h> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +namespace Gudhi { + + + /** + \class Witness_complex + \brief Constructs the witness complex for the given set of witnesses and landmarks. + \ingroup witness_complex + */ + template< class Simplicial_complex> + class Witness_complex { + + private: + + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) + {} + }; + + + private: + + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL; // Number of landmarks + double density; // Desired density + Simplicial_complex& sc; // Simplicial complex + + public: + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * Parameter dim serves as the limit for the number of closest landmarks to consider. + * Landmarks are supposed to be in [0,nbL_-1] + */ + template< typename KNearestNeighbours > + Witness_complex(KNearestNeighbours & knn, + Simplicial_complex & sc_, + int nbL_, + int dim ): nbL(nbL_), sc(sc_) + { + //Construction of the active witness list + int nbW = knn.size(); + typeVectorVertex vv; + typeSimplex simplex; + typePairSimplexBool returnValue; + int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w;// = new ActiveWitnessList(); + for (int i=0; i != nbL; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + counter++; + vv = {i}; + returnValue = sc.insert_simplex(vv, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + int k=1; /* current dimension in iterative construction */ + for (int i=0; i != nbW; ++i) + active_w.push_back(i); + //std::cout << "Successfully added edges" << std::endl; + while (!active_w.empty() && k < dim ) + { + //std::cout << "Started the step k=" << k << std::endl; + typename ActiveWitnessList::iterator it = active_w.begin(); + while (it != active_w.end()) + { + typeVectorVertex simplex_vector; + /* THE INSERTION: Checking if all the subfaces are in the simplex tree*/ + bool ok = all_faces_in(knn, *it, k); + if (ok) + { + for (int i = 0; i != k+1; ++i) + simplex_vector.push_back(knn[*it][i]); + returnValue = sc.insert_simplex(simplex_vector,0.0); + it++; + } + else + active_w.erase(it++); //First increase the iterator and then erase the previous element + } + k++; + } + } + + //@} + + private: + + /** \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + template <typename KNearestNeighbours> + bool all_faces_in(KNearestNeighbours &knn, int witness_id, int k) + { + //std::cout << "All face in with the landmark " << inserted_vertex << std::endl; + std::vector< Vertex_handle > facet; + //Vertex_handle curr_vh = curr_sh->first; + // CHECK ALL THE FACETS + for (int i = 0; i != k+1; ++i) + { + facet = {}; + for (int j = 0; j != k+1; ++j) + { + if (j != i) + { + facet.push_back(knn[witness_id][j]); + } + }//endfor + if (sc.find(facet) == sc.null_simplex()) + return false; + //std::cout << "++++ finished loop safely\n"; + } //endfor + return true; + } + + template <typename T> + void print_vector(std::vector<T> v) + { + 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 << "]"; + } + + public: + + /** + * \brief Verification if every simplex in the complex is witnessed by witnesses in knn. + * \param print_output =true will print the witnesses for each simplex + * \remark Added for debugging purposes. + */ + template< class KNearestNeighbors > + bool is_witness_complex(KNearestNeighbors & knn, bool print_output) + { + //bool final_result = true; + for (Simplex_handle sh: sc.complex_simplex_range()) + { + bool is_witnessed = false; + typeVectorVertex simplex; + int nbV = 0; //number of verticed in the simplex + for (int v: sc.simplex_vertex_range(sh)) + simplex.push_back(v); + nbV = simplex.size(); + for (typeVectorVertex w: knn) + { + bool has_vertices = true; + for (int v: simplex) + if (std::find(w.begin(), w.begin()+nbV, v) == w.begin()+nbV) + { + has_vertices = false; + //break; + } + if (has_vertices) + { + is_witnessed = true; + if (print_output) + { + std::cout << "The simplex "; + print_vector(simplex); + std::cout << " is witnessed by the witness "; + print_vector(w); + std::cout << std::endl; + } + break; + } + } + if (!is_witnessed) + { + if (print_output) + { + std::cout << "The following simplex is not witnessed "; + print_vector(simplex); + std::cout << std::endl; + } + assert(is_witnessed); + return false; + } + } + return true; // Arrive here if the not_witnessed check failed all the time + } + + +}; //class Witness_complex + + + +} // namespace Guhdi + +#endif diff --git a/src/Witness_complex/include/gudhi/Witness_complex_doc.h b/src/Witness_complex/include/gudhi/Witness_complex_doc.h new file mode 100644 index 00000000..71c8776b --- /dev/null +++ b/src/Witness_complex/include/gudhi/Witness_complex_doc.h @@ -0,0 +1,38 @@ +#ifndef WITNESS_COMPLEX_DOC_ +#define WITNESS_COMPLEX_DOC_ + +/** + \defgroup witness_complex Witness complex + + \author Siargey Kachanovich + + \section Definitions + + Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$: + + \li \f$W\f$ set of **witnesses** and + \li \f$L \subseteq W\f$ set of **landmarks**. + + The simplices are based on landmarks + and a simplex belongs to the witness complex if and only if it is witnessed, that is: + + \f$ \sigma \subset L \f$ is witnessed if there exists a point \f$w \in W\f$ such that + w is closer to the vertices of \f$ \sigma \f$ than other points in \f$ L \f$ and all of its faces are witnessed as well. + + \section Implementation + + The principal class of this module is Gudhi::Witness_complex. + + In both cases, the constructor for this class takes a {witness}x{closest_landmarks} table, where each row represents a witness and consists of landmarks sorted by distance to this witness. + This table can be constructed by two additional classes Landmark_choice_by_furthest_point and Landmark_choice_by_random_point also included in the module. + + *\image html "bench_Cy8.png" "Running time as function on number of landmarks" width=10cm + *\image html "bench_sphere.png" "Running time as function on number of witnesses for |L|=300" width=10cm + + + \copyright GNU General Public License v3. + + + */ + +#endif diff --git a/src/Witness_complex/test/CMakeLists.txt b/src/Witness_complex/test/CMakeLists.txt new file mode 100644 index 00000000..a13fd94a --- /dev/null +++ b/src/Witness_complex/test/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHITestWitnessComplex) + +#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") +#endif() + +add_executable ( simple_witness_complex simple_witness_complex.cpp ) +add_test(test ${CMAKE_CURRENT_BINARY_DIR}/simple_witness_complex) + +add_executable ( witness_complex_points witness_complex_points.cpp ) +add_test(test ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_points) + +#cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Simplex_tree/include/gudhi") diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp new file mode 100644 index 00000000..ea38b5c0 --- /dev/null +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -0,0 +1,54 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <ctime> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> + +using namespace Gudhi; + +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Witness_complex<Simplex_tree<>> WitnessComplex; +//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; +//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + +int main (int argc, char * const argv[]) +{ + Simplex_tree<> complex; + std::vector< typeVectorVertex > knn; + typeVectorVertex witness0 = {1,0,5,2,6,3,4}; knn.push_back(witness0 ); + typeVectorVertex witness1 = {2,6,4,5,0,1,3}; knn.push_back(witness1 ); + typeVectorVertex witness2 = {3,4,2,1,5,6,0}; knn.push_back(witness2 ); + typeVectorVertex witness3 = {4,2,1,3,5,6,0}; knn.push_back(witness3 ); + typeVectorVertex witness4 = {5,1,6,0,2,3,4}; knn.push_back(witness4 ); + typeVectorVertex witness5 = {6,0,5,2,1,3,4}; knn.push_back(witness5 ); + typeVectorVertex witness6 = {0,5,6,1,2,3,4}; knn.push_back(witness6 ); + typeVectorVertex witness7 = {2,6,4,5,3,1,0}; knn.push_back(witness7 ); + typeVectorVertex witness8 = {1,2,5,4,3,6,0}; knn.push_back(witness8 ); + typeVectorVertex witness9 = {3,4,0,6,5,1,2}; knn.push_back(witness9 ); + typeVectorVertex witness10 = {5,0,1,3,6,2,4}; knn.push_back(witness10); + typeVectorVertex witness11 = {5,6,1,0,2,3,4}; knn.push_back(witness11); + typeVectorVertex witness12 = {1,6,0,5,2,3,4}; knn.push_back(witness12); + WitnessComplex witnessComplex(knn, complex, 7, 7); + assert(witnessComplex.is_witness_complex(knn, true)); +} diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp new file mode 100644 index 00000000..0a50101d --- /dev/null +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -0,0 +1,63 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * 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 + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <ctime> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/Landmark_choice_by_furthest_point.h> + + +using namespace Gudhi; + +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Witness_complex<Simplex_tree<>> WitnessComplex; +typedef std::vector<double> Point; +//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; +//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + +int main (int argc, char * const argv[]) +{ + std::vector< typeVectorVertex > knn; + std::vector< Point > points; + // Add grid points as witnesses + for (double i = 0; i < 10; i += 1.0) + for (double j = 0; j < 10; j += 1.0) + for (double k = 0; k < 10; k += 1.0) + points.push_back(Point({i,j,k})); + + bool b_print_output = false; + // First test: random choice + Simplex_tree<> complex1; + Landmark_choice_by_random_point lcrp(points, 100, knn); + assert(!knn.empty()); + WitnessComplex witnessComplex1(knn, complex1, 100, 3); + assert(witnessComplex1.is_witness_complex(knn, b_print_output)); + + // Second test: furthest choice + knn.clear(); + Simplex_tree<> complex2; + Landmark_choice_by_furthest_point lcfp(points, 100, knn); + WitnessComplex witnessComplex2(knn, complex2, 100, 3); + assert(witnessComplex2.is_witness_complex(knn, b_print_output)); +} |
