diff options
Diffstat (limited to 'src/Witness_complex')
| -rw-r--r-- | src/Witness_complex/example/CMakeLists.txt | 11 | ||||
| -rw-r--r-- | src/Witness_complex/example/relaxed_witness_complex_sphere.cpp | 461 | ||||
| -rw-r--r-- | src/Witness_complex/example/simple_witness_complex.cpp | 52 | ||||
| -rw-r--r-- | src/Witness_complex/example/witness_complex_from_file.cpp | 156 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Relaxed_witness_complex.h | 886 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Witness_complex.h | 1111 |
6 files changed, 2677 insertions, 0 deletions
diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt new file mode 100644 index 00000000..83f9c71c --- /dev/null +++ b/src/Witness_complex/example/CMakeLists.txt @@ -0,0 +1,11 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIWitnessComplex) + +# A simple example + add_executable ( simple_witness_complex simple_witness_complex.cpp ) + add_test(simple_witness_complex ${CMAKE_CURRENT_BINARY_DIR}/simple_witness_complex) + + add_executable( witness_complex_from_file witness_complex_from_file.cpp ) + #target_link_libraries(witness_complex_from_file ${EIGEN3_LIBRARIES} ${CGAL_LIBRARY}) + add_test( witness_complex_from_bunny &{CMAKE_CURRENT_BINARY_DIR}/witness_complex_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + 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/simple_witness_complex.cpp b/src/Witness_complex/example/simple_witness_complex.cpp new file mode 100644 index 00000000..e95f67a8 --- /dev/null +++ b/src/Witness_complex/example/simple_witness_complex.cpp @@ -0,0 +1,52 @@ +/* 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/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" + +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[]) +{ + Witness_complex<> witnessComplex = Witness_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.witness_complex(KNN); +} 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..70c81528 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -0,0 +1,156 @@ +/* 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 <stdlib.h> + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +//#include <boost/filesystem.hpp> + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< std::vector <double> > Point_Vector; +//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; +//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + + +/** + * \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(); +} + +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) + { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file nbL \n"; + return 0; + } + /* + boost::filesystem::path p; + + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + std::string file_name = argv[1]; + int nbL = atoi(argv[2]); + + 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); + //std::cout << "Successfully read the points\n"; + witnessComplex.setNbL(nbL); + // witnessComplex.witness_complex_from_points(point_vector); + std::vector<std::vector< int > > WL; + std::set<int> L; + start = clock(); + //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); + witnessComplex.landmark_choice_by_random_points(point_vector, point_vector.size(), L); + witnessComplex.nearest_landmarks(point_vector,L,WL); + end = clock(); + std::cout << "Landmark choice for " << nbL << " landmarks took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + // Write the WL matrix in a file + mkdir("output", S_IRWXU); + const size_t last_slash_idx = file_name.find_last_of("/"); + if (std::string::npos != last_slash_idx) + { + file_name.erase(0, last_slash_idx + 1); + } + std::string out_file = "output/"+file_name+"_"+argv[2]+".wl"; + write_wl(out_file,WL); + start = clock(); + witnessComplex.witness_complex(WL); + // + end = clock(); + std::cout << "Howdy world! The process took " + << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; + /* + char buffer[100]; + int i = sprintf(buffer,"%s_%s_result.txt",argv[1],argv[2]); + 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(); + } + */ + + 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/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h new file mode 100644 index 00000000..c869628f --- /dev/null +++ b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h @@ -0,0 +1,886 @@ +/* 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_RELAXED_WITNESS_COMPLEX_H_ +#define GUDHI_RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#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_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ + template<typename FiltrationValue = double, + typename SimplexKey = int, + typename VertexHandle = int> + class Witness_complex: public Simplex_tree<> { + + private: + + struct Active_witness { + int witness_id; + int landmark_id; + Simplex_handle simplex_handle; + + Active_witness(int witness_id_, int landmark_id_, Simplex_handle simplex_handle_) + : witness_id(witness_id_), + landmark_id(landmark_id_), + simplex_handle(simplex_handle_) + {} + }; + + + + + public: + + + /** \brief Type for the vertex handle. + * + * Must be a signed integer type. It admits a total order <. */ + typedef VertexHandle Vertex_handle; + + /* Type of node in the simplex tree. */ + typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node; + /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ + typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary; + typedef typename Dictionary::iterator Simplex_handle; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + /** Number of landmarks + */ + int nbL; + /** Desired density + */ + double density; + + public: + + /** \brief Set number of landmarks to nbL_ + */ + void setNbL(int nbL_) + { + nbL = nbL_; + } + + /** \brief Set density to density_ + */ + void setDensity(double density_) + { + density = density_; + } + + /** + * /brief Iterative construction of the relaxed witness complex basing on a matrix of k nearest neighbours of the form {witnesses}x{landmarks} and (1+epsilon)-limit table {witnesses}*{landmarks} consisting of iterators of k nearest neighbor matrix. + * The line lengths can differ, however both matrices have the same corresponding line lengths. + */ + + template< typename KNearestNeighbours, typename OPELimits > + void relaxed_witness_complex(KNearestNeighbours & knn, OPELimits & rl) + //void witness_complex(std::vector< std::vector< Vertex_handle > > & knn) + { + std::cout << "**Start the procedure witness_complex" << std::endl; + //Construction of the active witness list + int nbW = knn.size(); + //int nbL = knn.at(0).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}; + 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 */ + //std::cout << "Successfully added landmarks" << std::endl; + // PRINT2 + //print_sc(root()); std::cout << std::endl; + for (int i=0; i != nbW; ++i) + active_w.push_back(i); + /* + int u,v; // two extremities of an edge + if (nbL > 1) // if the supposed dimension of the complex is >0 + { + for (int i=0; i != nbW; ++i) + { + // initial fill of active witnesses list + u = knn[i][0]; + v = knn[i][1]; + vv = {u,v}; + this->insert_simplex(vv,Filtration_value(0.0)); + //print_sc(root()); std::cout << std::endl; + //std::cout << "Added edges" << std::endl; + } + //print_sc(root()); + + } + */ + std::cout << "k=0, active witnesses: " << active_w.size() << std::endl; + //std::cout << "Successfully added edges" << std::endl; + //count_good = {0,0}; + //count_bad = {0,0}; + while (!active_w.empty() && k < nbL ) + { + //count_good.push_back(0); + //count_bad.push_back(0); + //std::cout << "Started the step k=" << k << std::endl; + typename ActiveWitnessList::iterator aw_it = active_w.begin(); + while (aw_it != active_w.end()) + { + std::vector<int> simplex; + bool ok = add_all_faces_of_dimension(k, knn[*aw_it].begin(), rl[*aw_it].begin(), simplex, knn[*aw_it].end(), knn[*aw_it].end()); + if (!ok) + active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + else + aw_it++; + } + std::cout << "k=" << k << ", active witnesses: " << active_w.size() << std::endl; + k++; + } + //print_sc(root()); std::cout << std::endl; + } + + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness + * Iterator is needed to know until how far we can take landmarks to form simplexes + * simplex is the prefix of the simplexes to insert + * The output value indicates if the witness rests active or not + */ + bool add_all_faces_of_dimension(int dim, std::vector<int>::iterator curr_l, typename std::vector< std::vector<int>::iterator >::iterator curr_until, std::vector<int>& simplex, std::vector<int>::iterator until, std::vector<int>::iterator end) + { + /* + std::ofstream ofs ("stree_result.txt", std::ofstream::out); + st_to_file(ofs); + ofs.close(); + */ + //print_sc(root()); + bool will_be_active = false; + if (dim > 0) + for (std::vector<int>::iterator it = curr_l; it != until && it != end; ++it, ++curr_until) + { + simplex.push_back(*it); + if (find(simplex) != null_simplex()) + will_be_active = will_be_active || add_all_faces_of_dimension(dim-1, it+1, curr_until+1, simplex, until, end); + simplex.pop_back(); + if (until == end) + until = *curr_until; + } + else if (dim == 0) + for (std::vector<int>::iterator it = curr_l; it != until && it != end; ++it, ++curr_until) + { + simplex.push_back(*it); + if (all_faces_in(simplex)) + { + will_be_active = true; + insert_simplex(simplex, 0.0); + } + simplex.pop_back(); + if (until == end) + until = *curr_until; + } + return will_be_active; + } + + /** \brief Construction of witness complex from points given explicitly + * nbL must be set to the right value of landmarks for strategies + * FURTHEST_POINT_STRATEGY and RANDOM_POINT_STRATEGY and + * density must be set to the right value for DENSITY_STRATEGY + */ + // void witness_complex_from_points(Point_Vector point_vector) + // { + // std::vector<std::vector< int > > WL; + // landmark_choice_by_random_points(point_vector, point_vector.size(), WL); + // witness_complex(WL); + // } + +private: + + /** \brief Print functions + */ + void print_sc(Siblings * sibl) + { + if (sibl == NULL) + std::cout << "&"; + else + print_children(sibl->members_); + } + + void print_children(Dictionary map) + { + std::cout << "("; + if (!map.empty()) + { + std::cout << map.begin()->first; + if (has_children(map.begin())) + print_sc(map.begin()->second.children()); + typename Dictionary::iterator it; + for (it = map.begin()+1; it != map.end(); ++it) + { + std::cout << "," << it->first; + if (has_children(it)) + print_sc(it->second.children()); + } + } + std::cout << ")"; + } + + public: + /** \brief Print functions + */ + + void st_to_file(std::ofstream& out_file) + { + sc_to_file(out_file, root()); + } + + private: + void sc_to_file(std::ofstream& out_file, Siblings * sibl) + { + assert(sibl); + children_to_file(out_file, sibl->members_); + } + + void children_to_file(std::ofstream& out_file, Dictionary& map) + { + out_file << "(" << std::flush; + if (!map.empty()) + { + out_file << map.begin()->first << std::flush; + if (has_children(map.begin())) + sc_to_file(out_file, map.begin()->second.children()); + typename Dictionary::iterator it; + for (it = map.begin()+1; it != map.end(); ++it) + { + out_file << "," << it->first << std::flush; + if (has_children(it)) + sc_to_file(out_file, it->second.children()); + } + } + out_file << ")" << std::flush; + } + + + /** \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 + */ + bool all_faces_in(std::vector<int>& simplex) + { + //std::cout << "All face in with the landmark " << inserted_vertex << std::endl; + std::vector< VertexHandle > facet; + //VertexHandle curr_vh = curr_sh->first; + // CHECK ALL THE FACETS + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + //facet = {}; + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + if (find(facet) == null_simplex()) + return false; + } //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 << "]"; + } + + template <typename T> + void print_vvector(std::vector< std::vector <T> > vv) + { + std::cout << "["; + if (!vv.empty()) + { + print_vector(*(vv.begin())); + for (auto it = vv.begin()+1; it != vv.end(); ++it) + { + std::cout << ","; + print_vector(*it); + } + } + std::cout << "]\n"; + } + + public: +/** + * \brief Landmark choice strategy by iteratively adding the landmark the furthest from the + * current landmark set + * \arg W is the vector of points which will be the witnesses + * \arg nbP is the number of witnesses + * \arg nbL is the number of landmarks + * \arg WL is the matrix of the nearest landmarks with respect to witnesses (output) + */ + + template <typename KNearestNeighbours> + void landmark_choice_by_furthest_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) + { + //std::cout << "Enter landmark_choice_by_furthest_points "<< std::endl; + //std::cout << "W="; print_vvector(W); + //double density = 5.; + Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks + typeVectorVertex chosen_landmarks; // landmark list + + WL = KNearestNeighbours(nbP,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 + double curr_dist; // used to stock the distance from the current point to L + double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nbP,infty); // vector of current distances to L from points + // double mindist = infty; + int curr_max_w=0; // the point currently furthest from L + int j; + int temp_swap_int; + double temp_swap_double; + + //CHOICE OF THE FIRST LANDMARK + std::cout << "Enter the first landmark stage\n"; + srand(354698); + int rand_int = rand()% nbP; + 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); + //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; + //std::cout << "WL="; print_vvector(WL); + //std::cout << "WLD="; print_vvector(wit_land_dist); + //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; + for (auto v: WL) + v.push_back(current_number_of_landmarks); + for (int i = 0; i < nbP; ++i) + { + // iteration on points in W. update of distance vectors + + //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; + //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; + curr_dist = euclidean_distance(W[i],W[chosen_landmarks[current_number_of_landmarks]]); + //std::cout << "The problem is not in distance function\n"; + wit_land_dist[i].push_back(curr_dist); + WL[i].push_back(current_number_of_landmarks); + //std::cout << "Push't back\n"; + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + j = current_number_of_landmarks; + //std::cout << "First half complete\n"; + while (j > 0 && wit_land_dist[i][j-1] > wit_land_dist[i][j]) + { + // sort the closest landmark vector for every witness + temp_swap_int = WL[i][j]; + WL[i][j] = WL[i][j-1]; + WL[i][j-1] = temp_swap_int; + temp_swap_double = wit_land_dist[i][j]; + wit_land_dist[i][j] = wit_land_dist[i][j-1]; + wit_land_dist[i][j-1] = temp_swap_double; + --j; + } + //std::cout << "result WL="; print_vvector(WL); + //std::cout << "result WLD="; print_vvector(wit_land_dist); + //std::cout << "result distL="; print_vector(dist_to_L); std::cout << std::endl; + //std::cout << "End loop\n"; + } + //std::cout << "Distance to landmarks="; print_vector(dist_to_L); std::cout << std::endl; + curr_max_dist = 0; + for (int i = 0; i < nbP; ++i) { + if (dist_to_L[i] > curr_max_dist) + { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + //std::cout << "Chose " << curr_max_w << " as new landmark\n"; + } + //std::cout << endl; + } + + /** \brief Landmark choice strategy by taking random vertices for landmarks. + * + */ + + // template <typename KNearestNeighbours> + // void landmark_choice_by_random_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) + // { + // std::cout << "Enter landmark_choice_by_random_points "<< std::endl; + // //std::cout << "W="; print_vvector(W); + // std::unordered_set< int > chosen_landmarks; // landmark set + + // Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks + + // WL = KNearestNeighbours(nbP,std::vector<int>()); + // int current_number_of_landmarks=0; // counter for landmarks + + // srand(24660); + // int chosen_landmark = rand()%nbP; + // double curr_dist; + + // //int j; + // //int temp_swap_int; + // //double temp_swap_double; + + + // for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + // { + // while (chosen_landmarks.find(chosen_landmark) != chosen_landmarks.end()) + // { + // srand((int)clock()); + // chosen_landmark = rand()% nbP; + // //std::cout << chosen_landmark << "\n"; + // } + // chosen_landmarks.insert(chosen_landmark); + // //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; + // //std::cout << "WL="; print_vvector(WL); + // //std::cout << "WLD="; print_vvector(wit_land_dist); + // //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; + // for (auto v: WL) + // v.push_back(current_number_of_landmarks); + // for (int i = 0; i < nbP; ++i) + // { + // // iteration on points in W. update of distance vectors + + // //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; + // //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; + // curr_dist = euclidean_distance(W[i],W[chosen_landmark]); + // //std::cout << "The problem is not in distance function\n"; + // wit_land_dist[i].push_back(curr_dist); + // WL[i].push_back(current_number_of_landmarks); + // //std::cout << "Push't back\n"; + // //j = current_number_of_landmarks; + // //std::cout << "First half complete\n"; + // //std::cout << "result WL="; print_vvector(WL); + // //std::cout << "result WLD="; print_vvector(wit_land_dist); + // //std::cout << "End loop\n"; + // } + // } + // for (int i = 0; i < nbP; i++) + // { + // sort(WL[i].begin(), WL[i].end(), [&](int j1, int j2){return wit_land_dist[i][j1] < wit_land_dist[i][j2];}); + // } + // //std::cout << endl; + // } + + /** \brief Landmark choice strategy by taking random vertices for landmarks. + * + */ + + // template <typename KNearestNeighbours> + void landmark_choice_by_random_points(Point_Vector &W, int nbP, std::set<int> &L) + { + std::cout << "Enter landmark_choice_by_random_points "<< std::endl; + //std::cout << "W="; print_vvector(W); + //std::unordered_set< int > chosen_landmarks; // landmark set + + //Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks + + //WL = KNearestNeighbours(nbP,std::vector<int>()); + int current_number_of_landmarks=0; // counter for landmarks + + srand(24660); + int chosen_landmark = rand()%nbP; + //double curr_dist; + //int j; + //int temp_swap_int; + //double temp_swap_double; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + { + while (L.find(chosen_landmark) != L.end()) + { + srand((int)clock()); + chosen_landmark = rand()% nbP; + //std::cout << chosen_landmark << "\n"; + } + L.insert(chosen_landmark); + //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; + //std::cout << "WL="; print_vvector(WL); + //std::cout << "WLD="; print_vvector(wit_land_dist); + //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; + // for (auto v: WL) + // v.push_back(current_number_of_landmarks); + // for (int i = 0; i < nbP; ++i) + // { + // // iteration on points in W. update of distance vectors + + // //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; + // //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; + // curr_dist = euclidean_distance(W[i],W[chosen_landmark]); + // //std::cout << "The problem is not in distance function\n"; + // wit_land_dist[i].push_back(curr_dist); + // WL[i].push_back(current_number_of_landmarks); + // //std::cout << "Push't back\n"; + // //j = current_number_of_landmarks; + // //std::cout << "First half complete\n"; + // //std::cout << "result WL="; print_vvector(WL); + // //std::cout << "result WLD="; print_vvector(wit_land_dist); + // //std::cout << "End loop\n"; + // } + } + // for (int i = 0; i < nbP; i++) + // { + // sort(WL[i].begin(), WL[i].end(), [&](int j1, int j2){return wit_land_dist[i][j1] < wit_land_dist[i][j2];}); + // } + //std::cout << endl; + } + + + /** \brief Construct the matrix |W|x(D+1) of D+1 closest landmarks + * where W is the set of witnesses and D is the ambient dimension + */ + template <typename KNearestNeighbours> + void nearest_landmarks(Point_Vector &W, std::set<int> &L, KNearestNeighbours &WL) + { + int D = W[0].size(); + int nbP = W.size(); + WL = KNearestNeighbours(nbP,std::vector<int>()); + typedef std::pair<double,int> dist_i; + typedef bool (*comp)(dist_i,dist_i); + for (int W_i = 0; W_i < nbP; W_i++) + { + //std::cout << "<<<<<<<<<<<<<<" << W_i <<"\n"; + 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 L_it; + int L_i; + for (L_it = L.begin(), L_i=0; L_it != L.end(); L_it++, L_i++) + { + dist_i dist = std::make_pair(euclidean_distance(W[W_i],W[*L_it]), L_i); + l_heap.push(dist); + } + for (int i = 0; i < D+1; i++) + { + dist_i dist = l_heap.top(); + WL[W_i].push_back(dist.second); + //WL[W_i].insert(WL[W_i].begin(),dist.second); + //std::cout << dist.first << " " << dist.second << std::endl; + l_heap.pop(); + } + } + } + + /** \brief Search and output links around vertices that are not pseudomanifolds + * + */ + void write_bad_links(std::ofstream& out_file) + { + out_file << "Bad links list\n"; + std::cout << "Entered write_bad_links\n"; + //typeVectorVertex testv = {9,15,17}; + //int count = 0; + for (auto v: complex_vertex_range()) + { + //std::cout << "Vertex " << v << ":\n"; + std::vector< Vertex_handle > link_vertices; + // Fill link_vertices + for (auto u: complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u != v && find(edge) != null_simplex()) + link_vertices.push_back(u); + } + /* + print_vector(link_vertices); + std::cout << "\n"; + */ + // Find the dimension + typeVectorVertex empty_simplex = {}; + int d = link_dim(link_vertices, link_vertices.begin(),-1, empty_simplex); + //std::cout << " dim " << d << "\n"; + //Siblings* curr_sibl = root(); + if (link_is_pseudomanifold(link_vertices,d)) + count_good[d]++; + //out_file << "Bad link at " << v << "\n"; + } + //out_file << "Number of bad links: " << count << "/" << root()->size(); + //std::cout << "Number of bad links: " << count << "/" << root()->size() << std::endl; + nc = nbL; + for (unsigned int i = 0; i != count_good.size(); i++) + { + out_file << "count_good[" << i << "] = " << count_good[i] << std::endl; + nc -= count_good[i]; + if (count_good[i] != 0) + std::cout << "count_good[" << i << "] = " << count_good[i] << std::endl; + } + for (unsigned int i = 0; i != count_bad.size(); i++) + { + out_file << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + nc -= count_bad[i]; + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + } + std::cout << "not_connected = " << nc << std::endl; + } + + private: + + std::vector<int> count_good; + std::vector<int> count_bad; + int nc; + + int link_dim(std::vector< Vertex_handle >& link_vertices, + typename std::vector< Vertex_handle >::iterator curr_v, + int curr_d, + typeVectorVertex& curr_simplex) + { + //std::cout << "Entered link_dim for " << *(curr_v-1) << "\n"; + Simplex_handle sh; + int final_d = curr_d; + typename std::vector< Vertex_handle >::iterator it; + for (it = curr_v; it != link_vertices.end(); ++it) + { + curr_simplex.push_back(*it); + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " final_dim " << final_d; + */ + sh = find(curr_simplex); + if (sh != null_simplex()) + { + //std::cout << " -> " << *it << "\n"; + int d = link_dim(link_vertices, it+1, curr_d+1, curr_simplex); + if (d > final_d) + final_d = d; + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); + } + return final_d; + } + + // color is false is a (d-1)-dim face, true is a d-dim face + //typedef bool Color; + // graph is an adjacency list + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph; + // map that gives to a certain simplex its node in graph and its dimension + //typedef std::pair<boost::vecS,Color> Reference; + typedef boost::graph_traits<Adj_graph>::vertex_descriptor Vertex_t; + typedef boost::graph_traits<Adj_graph>::edge_descriptor Edge_t; + + typedef boost::container::flat_map<Simplex_handle, Vertex_t> Graph_map; + + /* \brief Verifies if the simplices formed by vertices given by link_vertices + * form a pseudomanifold. + * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional + * faces and edges represent adjacency between them. + */ + bool link_is_pseudomanifold(std::vector< Vertex_handle >& link_vertices, + int dimension) + { + Adj_graph adj_graph; + Graph_map d_map, f_map; // d_map = map for d-dimensional simplices + // f_map = map for its facets + typeVectorVertex empty_vector = {}; + add_vertices(link_vertices, + link_vertices.begin(), + adj_graph, + d_map, + f_map, + empty_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges(adj_graph, d_map, f_map); + for (auto f_map_it : f_map) + { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2) + { + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + std::vector<int> components(boost::num_vertices(adj_graph)); + return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + void add_vertices(typeVectorVertex& link_vertices, + typename typeVectorVertex::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + typeVectorVertex& curr_simplex, + int curr_d, + int dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename typeVectorVertex::iterator it; + std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != link_vertices.end(); ++it) + { + curr_simplex.push_back(*it); + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + sh = find(curr_simplex); + if (sh != null_simplex()) + { + //std::cout << " added\n"; + if (curr_d == dimension) + { + vert = boost::add_vertex(adj_graph); + resPair = d_map.emplace(sh,vert); + } + else + { + if (curr_d == dimension-1) + { + vert = boost::add_vertex(adj_graph); + resPair = f_map.emplace(sh,vert); + } + add_vertices(link_vertices, + it+1, + adj_graph, + d_map, + f_map, + curr_simplex, + curr_d+1, dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); + } + } + + void add_edges(Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map) + { + Simplex_handle sh; + // Add edges + //std::cout << "Entered add edges:\n"; + typename Graph_map::iterator map_it; + for (auto d_map_pair : d_map) + { + //std::cout << "*"; + sh = d_map_pair.first; + Vertex_t d_vert = d_map_pair.second; + for (auto facet_sh : boundary_simplex_range(sh)) + //for (auto f_map_it : f_map) + { + //std::cout << "'"; + map_it = f_map.find(facet_sh); + //We must have all the facets in the graph at this point + assert(map_it != f_map.end()); + Vertex_t f_vert = map_it->second; + //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; + boost::add_edge(d_vert,f_vert,adj_graph); + } + } + } + + ////////////////////////////////////////////////////////////////////////////////////////////////// + //***********COLLAPSES**************************************************************************// + ////////////////////////////////////////////////////////////////////////////////////////////////// + + + + + + + +}; //class Witness_complex + + + +} // namespace Guhdi + +#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..201d6525 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -0,0 +1,1111 @@ +/* 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/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +//#include <CGAL/Delaunay_triangulation.h> +//#include <CGAL/Epick_d.h> +//#include <CGAL/K_neighbor_search.h> +//#include <CGAL/Search_traits_d.h> + +// 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 { + + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ + template<typename FiltrationValue = double, + typename SimplexKey = int, + typename VertexHandle = int> + class Witness_complex: public Simplex_tree<> { + + private: + + struct Active_witness { + int witness_id; + int landmark_id; + Simplex_handle simplex_handle; + + Active_witness(int witness_id_, int landmark_id_, Simplex_handle simplex_handle_) + : witness_id(witness_id_), + landmark_id(landmark_id_), + simplex_handle(simplex_handle_) + {} + }; + + + + + public: + + + /** \brief Type for the vertex handle. + * + * Must be a signed integer type. It admits a total order <. */ + typedef VertexHandle Vertex_handle; + + /* Type of node in the simplex tree. */ + typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node; + /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ + typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary; + typedef typename Dictionary::iterator Simplex_handle; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + /** Number of landmarks + */ + int nbL; + /** Desired density + */ + double density; + + public: + + /** \brief Set number of landmarks to nbL_ + */ + void setNbL(int nbL_) + { + nbL = nbL_; + } + + /** \brief Set density to density_ + */ + void setDensity(double density_) + { + density = density_; + } + + /** + * /brief Iterative construction of the witness complex basing on a matrix of k nearest neighbours of the form {witnesses}x{landmarks}. + * Landmarks are supposed to be in [0,nbL-1] + */ + + template< typename KNearestNeighbours > + void witness_complex(KNearestNeighbours & knn) + //void witness_complex(std::vector< std::vector< Vertex_handle > > & knn) + { + std::cout << "**Start the procedure witness_complex" << std::endl; + //Construction of the active witness list + int nbW = knn.size(); + //int nbL = knn.at(0).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 = 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 */ + //std::cout << "Successfully added landmarks" << std::endl; + // PRINT2 + //print_sc(root()); std::cout << std::endl; + /* + int u,v; // two extremities of an edge + int count = 0; + if (nbL > 1) // if the supposed dimension of the complex is >0 + { + for (int i=0; i != nbW; ++i) + { + // initial fill of active witnesses list + u = knn[i][0]; + v = knn[i][1]; + vv = {u,v}; + returnValue = this->insert_simplex(vv,Filtration_value(0.0)); + if (returnValue.second) + count++; + //print_sc(root()); std::cout << std::endl; + //std::cout << "Added edges" << std::endl; + } + std::cout << "The number of edges = " << count << std::endl; + count = 0; + //print_sc(root()); + for (int i=0; i != nbW; ++i) + { + // initial fill of active witnesses list + u = knn[i][0]; + v = knn[i][1]; + if ( u > v) + { + u = v; + v = knn[i][0]; + knn[i][0] = knn[i][1]; + knn[i][1] = v; + } + Simplex_handle sh; + vv = {u,v}; + //if (u==v) std::cout << "Bazzinga!\n"; + sh = (root()->find(u))->second.children()->find(v); + active_w.push_back(i); + } + } + */ + for (int i=0; i != nbW; ++i) + active_w.push_back(i); + std::cout << "k=0, active witnesses: " << active_w.size() << std::endl; + //std::cout << "Successfully added edges" << std::endl; + count_good = {0}; + count_bad = {0}; + int D = knn[0].size(); + while (!active_w.empty() && k < D ) + { + count_good.push_back(0); + count_bad.push_back(0); + //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 = insert_simplex(simplex_vector,0.0); + it++; + } + else + active_w.erase(it++); //First increase the iterator and then erase the previous element + } + std::cout << "k=" << k << ", active witnesses: " << active_w.size() << std::endl; + //std::cout << "** k=" << k << ", num_simplices: " <<count << std::endl; + k++; + } + //print_sc(root()); std::cout << std::endl; + } + + /** \brief Construction of witness complex from points given explicitly + * nbL must be set to the right value of landmarks for strategies + * FURTHEST_POINT_STRATEGY and RANDOM_POINT_STRATEGY and + * density must be set to the right value for DENSITY_STRATEGY + */ + // void witness_complex_from_points(Point_Vector point_vector) + // { + // std::vector<std::vector< int > > WL; + // landmark_choice_by_random_points(point_vector, point_vector.size(), WL); + // witness_complex(WL); + // } + +private: + + /** \brief Print functions + */ + void print_sc(Siblings * sibl) + { + if (sibl == NULL) + std::cout << "&"; + else + print_children(sibl->members_); + } + + void print_children(Dictionary map) + { + std::cout << "("; + if (!map.empty()) + { + std::cout << map.begin()->first; + if (has_children(map.begin())) + print_sc(map.begin()->second.children()); + typename Dictionary::iterator it; + for (it = map.begin()+1; it != map.end(); ++it) + { + std::cout << "," << it->first; + if (has_children(it)) + print_sc(it->second.children()); + } + } + std::cout << ")"; + } + + public: + /** \brief Print functions + */ + + void st_to_file(std::ofstream& out_file) + { + sc_to_file(out_file, root()); + } + + private: + void sc_to_file(std::ofstream& out_file, Siblings * sibl) + { + assert(sibl); + children_to_file(out_file, sibl->members_); + } + + void children_to_file(std::ofstream& out_file, Dictionary& map) + { + out_file << "(" << std::flush; + if (!map.empty()) + { + out_file << map.begin()->first << std::flush; + if (has_children(map.begin())) + sc_to_file(out_file, map.begin()->second.children()); + typename Dictionary::iterator it; + for (it = map.begin()+1; it != map.end(); ++it) + { + out_file << "," << it->first << std::flush; + if (has_children(it)) + sc_to_file(out_file, it->second.children()); + } + } + out_file << ")" << std::flush; + } + + + /** \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< VertexHandle > facet; + //VertexHandle 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 (find(facet) == 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 << "]"; + } + + template <typename T> + void print_vvector(std::vector< std::vector <T> > vv) + { + std::cout << "["; + if (!vv.empty()) + { + print_vector(*(vv.begin())); + for (auto it = vv.begin()+1; it != vv.end(); ++it) + { + std::cout << ","; + print_vector(*it); + } + } + std::cout << "]\n"; + } + + public: +/** + * \brief Landmark choice strategy by iteratively adding the landmark the furthest from the + * current landmark set + * \arg W is the vector of points which will be the witnesses + * \arg nbP is the number of witnesses + * \arg nbL is the number of landmarks + * \arg WL is the matrix of the nearest landmarks with respect to witnesses (output) + */ + + template <typename KNearestNeighbours> + void landmark_choice_by_furthest_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) + { + //std::cout << "Enter landmark_choice_by_furthest_points "<< std::endl; + //std::cout << "W="; print_vvector(W); + //double density = 5.; + Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks + typeVectorVertex chosen_landmarks; // landmark list + + WL = KNearestNeighbours(nbP,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 + double curr_dist; // used to stock the distance from the current point to L + double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nbP,infty); // vector of current distances to L from points + // double mindist = infty; + int curr_max_w=0; // the point currently furthest from L + int j; + int temp_swap_int; + double temp_swap_double; + + //CHOICE OF THE FIRST LANDMARK + std::cout << "Enter the first landmark stage\n"; + srand(354698); + int rand_int = rand()% nbP; + 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); + //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; + //std::cout << "WL="; print_vvector(WL); + //std::cout << "WLD="; print_vvector(wit_land_dist); + //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; + for (auto v: WL) + v.push_back(current_number_of_landmarks); + for (int i = 0; i < nbP; ++i) + { + // iteration on points in W. update of distance vectors + + //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; + //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; + curr_dist = euclidean_distance(W[i],W[chosen_landmarks[current_number_of_landmarks]]); + //std::cout << "The problem is not in distance function\n"; + wit_land_dist[i].push_back(curr_dist); + WL[i].push_back(current_number_of_landmarks); + //std::cout << "Push't back\n"; + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + j = current_number_of_landmarks; + //std::cout << "First half complete\n"; + while (j > 0 && wit_land_dist[i][j-1] > wit_land_dist[i][j]) + { + // sort the closest landmark vector for every witness + temp_swap_int = WL[i][j]; + WL[i][j] = WL[i][j-1]; + WL[i][j-1] = temp_swap_int; + temp_swap_double = wit_land_dist[i][j]; + wit_land_dist[i][j] = wit_land_dist[i][j-1]; + wit_land_dist[i][j-1] = temp_swap_double; + --j; + } + //std::cout << "result WL="; print_vvector(WL); + //std::cout << "result WLD="; print_vvector(wit_land_dist); + //std::cout << "result distL="; print_vector(dist_to_L); std::cout << std::endl; + //std::cout << "End loop\n"; + } + //std::cout << "Distance to landmarks="; print_vector(dist_to_L); std::cout << std::endl; + curr_max_dist = 0; + for (int i = 0; i < nbP; ++i) { + if (dist_to_L[i] > curr_max_dist) + { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + //std::cout << "Chose " << curr_max_w << " as new landmark\n"; + } + //std::cout << endl; + } + + /** \brief Landmark choice strategy by taking random vertices for landmarks. + * + */ + + // template <typename KNearestNeighbours> + // void landmark_choice_by_random_points(Point_Vector &W, int nbP, KNearestNeighbours &WL) + // { + // std::cout << "Enter landmark_choice_by_random_points "<< std::endl; + // //std::cout << "W="; print_vvector(W); + // std::unordered_set< int > chosen_landmarks; // landmark set + + // Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks + + // WL = KNearestNeighbours(nbP,std::vector<int>()); + // int current_number_of_landmarks=0; // counter for landmarks + + // srand(24660); + // int chosen_landmark = rand()%nbP; + // double curr_dist; + + // //int j; + // //int temp_swap_int; + // //double temp_swap_double; + + + // for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + // { + // while (chosen_landmarks.find(chosen_landmark) != chosen_landmarks.end()) + // { + // srand((int)clock()); + // chosen_landmark = rand()% nbP; + // //std::cout << chosen_landmark << "\n"; + // } + // chosen_landmarks.insert(chosen_landmark); + // //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; + // //std::cout << "WL="; print_vvector(WL); + // //std::cout << "WLD="; print_vvector(wit_land_dist); + // //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; + // for (auto v: WL) + // v.push_back(current_number_of_landmarks); + // for (int i = 0; i < nbP; ++i) + // { + // // iteration on points in W. update of distance vectors + + // //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; + // //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; + // curr_dist = euclidean_distance(W[i],W[chosen_landmark]); + // //std::cout << "The problem is not in distance function\n"; + // wit_land_dist[i].push_back(curr_dist); + // WL[i].push_back(current_number_of_landmarks); + // //std::cout << "Push't back\n"; + // //j = current_number_of_landmarks; + // //std::cout << "First half complete\n"; + // //std::cout << "result WL="; print_vvector(WL); + // //std::cout << "result WLD="; print_vvector(wit_land_dist); + // //std::cout << "End loop\n"; + // } + // } + // for (int i = 0; i < nbP; i++) + // { + // sort(WL[i].begin(), WL[i].end(), [&](int j1, int j2){return wit_land_dist[i][j1] < wit_land_dist[i][j2];}); + // } + // //std::cout << endl; + // } + + /** \brief Landmark choice strategy by taking random vertices for landmarks. + * + */ + + // template <typename KNearestNeighbours> + void landmark_choice_by_random_points(Point_Vector &W, int nbP, std::set<int> &L) + { + std::cout << "Enter landmark_choice_by_random_points "<< std::endl; + //std::cout << "W="; print_vvector(W); + //std::unordered_set< int > chosen_landmarks; // landmark set + + //Point_Vector wit_land_dist(nbP,std::vector<double>()); // distance matrix witness x landmarks + + //WL = KNearestNeighbours(nbP,std::vector<int>()); + int current_number_of_landmarks=0; // counter for landmarks + + srand(24660); + int chosen_landmark = rand()%nbP; + //double curr_dist; + //int j; + //int temp_swap_int; + //double temp_swap_double; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) + { + while (L.find(chosen_landmark) != L.end()) + { + srand((int)clock()); + chosen_landmark = rand()% nbP; + //std::cout << chosen_landmark << "\n"; + } + L.insert(chosen_landmark); + //std::cout << "**********Entered loop with current number of landmarks = " << current_number_of_landmarks << std::endl; + //std::cout << "WL="; print_vvector(WL); + //std::cout << "WLD="; print_vvector(wit_land_dist); + //std::cout << "landmarks="; print_vector(chosen_landmarks); std::cout << std::endl; + // for (auto v: WL) + // v.push_back(current_number_of_landmarks); + // for (int i = 0; i < nbP; ++i) + // { + // // iteration on points in W. update of distance vectors + + // //std::cout << "In the loop with i=" << i << " and landmark=" << chosen_landmarks[current_number_of_landmarks] << std::endl; + // //std::cout << "W[i]="; print_vector(W[i]); std::cout << " W[landmark]="; print_vector(W[chosen_landmarks[current_number_of_landmarks]]); std::cout << std::endl; + // curr_dist = euclidean_distance(W[i],W[chosen_landmark]); + // //std::cout << "The problem is not in distance function\n"; + // wit_land_dist[i].push_back(curr_dist); + // WL[i].push_back(current_number_of_landmarks); + // //std::cout << "Push't back\n"; + // //j = current_number_of_landmarks; + // //std::cout << "First half complete\n"; + // //std::cout << "result WL="; print_vvector(WL); + // //std::cout << "result WLD="; print_vvector(wit_land_dist); + // //std::cout << "End loop\n"; + // } + } + // for (int i = 0; i < nbP; i++) + // { + // sort(WL[i].begin(), WL[i].end(), [&](int j1, int j2){return wit_land_dist[i][j1] < wit_land_dist[i][j2];}); + // } + //std::cout << endl; + } + + + /** \brief Construct the matrix |W|x(D+1) of D+1 closest landmarks + * where W is the set of witnesses and D is the ambient dimension + */ + template <typename KNearestNeighbours> + void nearest_landmarks(Point_Vector &W, std::set<int> &L, KNearestNeighbours &WL) + { + int D = W[0].size(); + int nbP = W.size(); + WL = KNearestNeighbours(nbP,std::vector<int>()); + typedef std::pair<double,int> dist_i; + typedef bool (*comp)(dist_i,dist_i); + for (int W_i = 0; W_i < nbP; W_i++) + { + //std::cout << "<<<<<<<<<<<<<<" << W_i <<"\n"; + 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 L_it; + int L_i; + for (L_it = L.begin(), L_i=0; L_it != L.end(); L_it++, L_i++) + { + dist_i dist = std::make_pair(euclidean_distance(W[W_i],W[*L_it]), L_i); + l_heap.push(dist); + } + for (int i = 0; i < D+1; i++) + { + dist_i dist = l_heap.top(); + WL[W_i].push_back(dist.second); + //WL[W_i].insert(WL[W_i].begin(),dist.second); + //std::cout << dist.first << " " << dist.second << std::endl; + l_heap.pop(); + } + } + } + + /** \brief Returns true if the link is good + */ + bool has_good_link(Vertex_handle v, std::vector< int >& bad_count, std::vector< int >& good_count) + { + std::vector< Vertex_handle > star_vertices; + // Fill star_vertices + star_vertices.push_back(v); + for (auto u: complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u != v && find(edge) != null_simplex()) + star_vertices.push_back(u); + } + // Find the dimension + typeVectorVertex init_simplex = {star_vertices[0]}; + bool is_pure = true; + std::vector<int> dim_coface(star_vertices.size(), 1); + int d = star_dim(star_vertices, star_vertices.begin()+1, 0, init_simplex, dim_coface.begin()+1) - 1; //link_dim = star_dim - 1 + assert(init_simplex.size() == 1); + if (!is_pure) + std::cout << "Found an impure star around " << v << "\n"; + for (int dc: dim_coface) + is_pure = (dc == dim_coface[0]); + /* + if (d == count_good.size()) + { + std::cout << "Found a star of dimension " << (d+1) << " around " << v << "\nThe star is "; + print_vector(star_vertices); std::cout << std::endl; + } + */ + //if (d == -1) bad_count[0]++; + bool b= (is_pure && link_is_pseudomanifold(star_vertices,d)); + if (d != -1) {if (b) good_count[d]++; else bad_count[d]++;} + if (!is_pure) bad_count[0]++; + return (d != -1 && b && is_pure); + + } + + /** \brief Search and output links around vertices that are not pseudomanifolds + * + */ + /* + void write_bad_links(std::ofstream& out_file) + { + out_file << "Bad links list\n"; + std::cout << "Entered write_bad_links\n"; + for (auto v: complex_vertex_range()) + { + std::cout << "Vertex " << v << ": "; + std::vector< Vertex_handle > link_vertices; + // Fill link_vertices + for (auto u: complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u != v && find(edge) != null_simplex()) + link_vertices.push_back(u); + } + + print_vector(link_vertices); + std::cout << "\n"; + + // Find the dimension + typeVectorVertex empty_simplex = {}; + int d = link_dim(link_vertices, link_vertices.begin(),-1, empty_simplex); + if (link_is_pseudomanifold(link_vertices,d)) + count_good[d]++; + } + nc = nbL; + for (unsigned int i = 0; i != count_good.size(); i++) + { + out_file << "count_good[" << i << "] = " << count_good[i] << std::endl; + nc -= count_good[i]; + if (count_good[i] != 0) + std::cout << "count_good[" << i << "] = " << count_good[i] << std::endl; + } + for (unsigned int i = 0; i != count_bad.size(); i++) + { + out_file << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + nc -= count_bad[i]; + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + } + std::cout << "not_connected = " << nc << std::endl; + } + */ + private: + + std::vector<int> count_good; + std::vector<int> count_bad; + int nc; + + int star_dim(std::vector< Vertex_handle >& star_vertices, + typename std::vector< Vertex_handle >::iterator curr_v, + int curr_d, + typeVectorVertex& curr_simplex, + typename std::vector< int >::iterator curr_dc) + { + //std::cout << "Entered star_dim for " << *(curr_v-1) << "\n"; + Simplex_handle sh; + int final_d = curr_d; + typename std::vector< Vertex_handle >::iterator it; + typename std::vector< Vertex_handle >::iterator dc_it; + //std::cout << "Current vertex is " << + for (it = curr_v, dc_it = curr_dc; it != star_vertices.end(); ++it, ++dc_it) + { + curr_simplex.push_back(*it); + typeVectorVertex curr_simplex_copy(curr_simplex); + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " final_dim " << final_d; + */ + sh = find(curr_simplex_copy); //Need a copy because find sorts the vector and I want star center to be the first + if (sh != null_simplex()) + { + //std::cout << " -> " << *it << "\n"; + int d = star_dim(star_vertices, it+1, curr_d+1, curr_simplex, dc_it); + if (d >= final_d) + { + final_d = d; + //std::cout << d << " "; + //print_vector(curr_simplex); + //std::cout << std::endl; + } + if (d >= *dc_it) + *dc_it = d; + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); + } + return final_d; + } + + // color is false is a (d-1)-dim face, true is a d-dim face + //typedef bool Color; + // graph is an adjacency list + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph; + // map that gives to a certain simplex its node in graph and its dimension + //typedef std::pair<boost::vecS,Color> Reference; + typedef boost::graph_traits<Adj_graph>::vertex_descriptor Vertex_t; + typedef boost::graph_traits<Adj_graph>::edge_descriptor Edge_t; + typedef boost::graph_traits<Adj_graph>::adjacency_iterator Adj_it; + typedef std::pair<Adj_it, Adj_it> Out_edge_it; + + typedef boost::container::flat_map<Simplex_handle, Vertex_t> Graph_map; + typedef boost::container::flat_map<Vertex_t, Simplex_handle> Inv_graph_map; + + /* \brief Verifies if the simplices formed by vertices given by link_vertices + * form a pseudomanifold. + * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional + * faces and edges represent adjacency between them. + */ + bool link_is_pseudomanifold(std::vector< Vertex_handle >& star_vertices, + int dimension) + { + Adj_graph adj_graph; + Graph_map d_map, f_map; // d_map = map for d-dimensional simplices + // f_map = map for its facets + typeVectorVertex init_vector = {}; + add_vertices_to_link_graph(star_vertices, + star_vertices.begin()+1, + adj_graph, + d_map, + f_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, d_map, f_map); + for (auto f_map_it : f_map) + { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2) + { + /* + if (boost::out_degree(f_map_it.second, adj_graph) >= 3) + { + std::cout << "This simplex has 3+ cofaces: "; + for(auto v : simplex_vertex_range(f_map_it.first)) + std::cout << v << " "; + std::cout << std::endl; + Adj_it ai, ai_end; + for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) + { + + } + } + */ + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector<int> components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + public: +bool complex_is_pseudomanifold(int dimension) + { + Adj_graph adj_graph; + Graph_map d_map, f_map; // d_map = map for d-dimensional simplices + // f_map = map for its facets + Inv_graph_map inv_d_map; + typeVectorVertex init_vector = {}; + std::vector<int> star_vertices; + for (int v: complex_vertex_range()) + star_vertices.push_back(v); + add_max_simplices_to_graph(star_vertices, + star_vertices.begin(), + adj_graph, + d_map, + f_map, + inv_d_map, + init_vector, + 0, dimension); + std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, d_map, f_map); + for (auto f_map_it : f_map) + { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2) + { + if (boost::out_degree(f_map_it.second, adj_graph) >= 3) + { + std::cout << "This simplex has 3+ cofaces: "; + for(auto v : simplex_vertex_range(f_map_it.first)) + std::cout << v << " "; + std::cout << std::endl; + Adj_it ai, ai_end; + for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) + { + auto it = inv_d_map.find(*ai); + assert (it != inv_d_map.end()); + Simplex_handle sh = it->second; + for(auto v : simplex_vertex_range(sh)) + std::cout << v << " "; + std::cout << std::endl; + } + } + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector<int> components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + private: + void add_vertices_to_link_graph(typeVectorVertex& star_vertices, + typename typeVectorVertex::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + typeVectorVertex& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename typeVectorVertex::iterator it; + //std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) + { + curr_simplex.push_back(*it); //push next vertex in question + curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + typeVectorVertex curr_simplex_copy(curr_simplex); + sh = find(curr_simplex_copy); //a simplex of the star + curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = typeVectorVertex(curr_simplex); + if (sh != null_simplex()) + { + //std::cout << " added\n"; + if (curr_d == link_dimension) + { + sh = find(curr_simplex_copy); //a simplex of the link + assert(sh != null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + } + else + { + + if (curr_d == link_dimension-1) + { + sh = find(curr_simplex_copy); //a simplex of the link + assert(sh != null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_vertices_to_link_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + curr_simplex, + curr_d+1, link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + + void add_edges_to_link_graph(Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map) + { + Simplex_handle sh; + // Add edges + //std::cout << "Entered add edges:\n"; + typename Graph_map::iterator map_it; + for (auto d_map_pair : d_map) + { + //std::cout << "*"; + sh = d_map_pair.first; + Vertex_t d_vert = d_map_pair.second; + for (auto facet_sh : boundary_simplex_range(sh)) + //for (auto f_map_it : f_map) + { + //std::cout << "'"; + map_it = f_map.find(facet_sh); + //We must have all the facets in the graph at this point + assert(map_it != f_map.end()); + Vertex_t f_vert = map_it->second; + //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; + boost::add_edge(d_vert,f_vert,adj_graph); + } + } + } + + void add_max_simplices_to_graph(typeVectorVertex& star_vertices, + typename typeVectorVertex::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Inv_graph_map& inv_d_map, + typeVectorVertex& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename typeVectorVertex::iterator it; + //std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) + { + curr_simplex.push_back(*it); //push next vertex in question + //curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + typeVectorVertex curr_simplex_copy(curr_simplex); + sh = find(curr_simplex_copy); //a simplex of the star + //curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = typeVectorVertex(curr_simplex); + if (sh != null_simplex()) + { + //std::cout << " added\n"; + if (curr_d == link_dimension) + { + sh = find(curr_simplex_copy); //a simplex of the link + assert(sh != null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + inv_d_map.emplace(vert,sh); + } + else + { + + if (curr_d == link_dimension-1) + { + sh = find(curr_simplex_copy); //a simplex of the link + assert(sh != null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_max_simplices_to_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + inv_d_map, + curr_simplex, + curr_d+1, link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + + public: + /** \brief Verification if every simplex in the complex is witnessed + */ + template< class KNearestNeighbors > + bool is_witness_complex(KNearestNeighbors WL) + { + //bool final_result = true; + for (Simplex_handle sh: complex_simplex_range()) + { + bool is_witnessed = false; + typeVectorVertex simplex; + int nbV = 0; //number of verticed in the simplex + for (int v: simplex_vertex_range(sh)) + simplex.push_back(v); + nbV = simplex.size(); + for (typeVectorVertex w: WL) + { + 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; + 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) + { + 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 |