From 5cae0e771c7d4603d25ccf06c782fe383f91e74a Mon Sep 17 00:00:00 2001 From: skachano Date: Fri, 29 May 2015 09:22:36 +0000 Subject: Added the test with spheres git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@598 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f35dd2d082e318b912f581ee86e1407fd39109eb --- .../example/witness_complex_sphere.cpp | 748 +++++++++++++++++++++ 1 file changed, 748 insertions(+) create mode 100644 src/Witness_complex/example/witness_complex_sphere.cpp (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') 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..c9a9119a --- /dev/null +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -0,0 +1,748 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +//#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +//#include + +//#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d 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 Euclidean_distance; + +/** + * \brief Class of distance in a flat torus in dimaension D + * + */ +//class Torus_distance : public Euclidean_distance { + class Torus_distance { + +public: + typedef K::FT FT; + typedef K::Point_d Point_d; + typedef Point_d Query_item; + typedef typename CGAL::Dynamic_dimension_tag D; + + double box_length = 2; + + FT transformed_distance(Query_item q, Point_d p) const + { + FT distance = FT(0); + FT coord = FT(0); + //std::cout << "Hello skitty!\n"; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1), pit = construct_it(p); + for(; qit != qe; qit++, pit++) + { + coord = sqrt(((*qit)-(*pit))*((*qit)-(*pit))); + if (coord*coord <= (box_length-coord)*(box_length-coord)) + distance += coord*coord; + else + distance += (box_length-coord)*(box_length-coord); + } + return distance; + } + + FT min_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r) const { + FT distance = FT(0); + FT dist1, dist2; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if((*qit) < r.min_coord(i)) + { + dist1 = (r.min_coord(i)-(*qit)); + dist2 = (box_length - r.max_coord(i)+(*qit)); + if (dist1 < dist2) + distance += dist1*dist1; + else + distance += dist2*dist2; + } + else if ((*qit) > r.max_coord(i)) + { + dist1 = (box_length - (*qit)+r.min_coord(i)); + dist2 = ((*qit) - r.max_coord(i)); + if (dist1 < dist2) + distance += dist1*dist1; + else + distance += dist2*dist2; + } + } + return distance; + } + + FT min_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r, + std::vector& dists) const { + FT distance = FT(0); + FT dist1, dist2; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + //std::cout << r.max_coord(0) << std::endl; + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if((*qit) < r.min_coord(i)) + { + dist1 = (r.min_coord(i)-(*qit)); + dist2 = (box_length - r.max_coord(i)+(*qit)); + if (dist1 < dist2) + { + dists[i] = dist1; + distance += dist1*dist1; + } + else + { + dists[i] = dist2; + distance += dist2*dist2; + //std::cout << "Good stuff1\n"; + } + } + else if ((*qit) > r.max_coord(i)) + { + dist1 = (box_length - (*qit)+r.min_coord(i)); + dist2 = ((*qit) - r.max_coord(i)); + if (dist1 < dist2) + { + dists[i] = dist1; + distance += dist1*dist1; + //std::cout << "Good stuff2\n"; + } + else + { + dists[i] = dist2; + distance += dist2*dist2; + } + } + }; + return distance; + } + + FT max_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r) const { + FT distance=FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if (box_length <= (r.min_coord(i)+r.max_coord(i))) + if ((r.max_coord(i)+r.min_coord(i)-box_length)/FT(2.0) <= (*qit) && + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + else + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + else + if ((box_length-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + else + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + return distance; + } + + + FT max_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r, + std::vector& dists) const { + FT distance=FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if (box_length <= (r.min_coord(i)+r.max_coord(i))) + if ((r.max_coord(i)+r.min_coord(i)-box_length)/FT(2.0) <= (*qit) && + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + { + dists[i] = r.max_coord(i)-(*qit); + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + } + else + { + dists[i] = sqrt(((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i))); + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + else + if ((box_length-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + { + dists[i] = sqrt((r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit))); + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + + } + else + { + dists[i] = (*qit)-r.min_coord(i); + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + } + return distance; + } + + inline FT new_distance(FT dist, FT old_off, FT new_off, + int /* cutting_dimension */) const { + + FT new_dist = dist + (new_off*new_off - old_off*old_off); + return new_dist; + } + + inline FT transformed_distance(FT d) const { + return d*d; + } + + inline FT inverse_of_transformed_distance(FT d) const { + return sqrt(d); + } + +}; + + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair 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 Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef boost::container::flat_map Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +//typedef K::Equal_d Equal_d; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; + +bool toric=false; + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , Point_Vector & points) +{ + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector< double > point; + std::istringstream iss( line ); + while(iss >> x) { point.push_back(x); } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +void generate_points_grid(Point_Vector& W, int width, int D) +{ + +} + +void generate_points_random_box(Point_Vector& W, int nbP, int dim) +{ + Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + +/* NOT TORUS RELATED + */ +void generate_points_sphere(Point_Vector& W, int nbP, int dim) +{ + CGAL::Random_points_on_sphere_d rp(dim,1); + for (int i = 0; i < nbP; i++) + W.push_back(*rp++); +} +/* +void read_points_to_tree (std::string file_name, Tree& tree) +{ + //I assume here that tree is empty + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector coords; + std::istringstream iss( line ); + while(iss >> x) { coords.push_back(x); } + if (coords.size() != 1) + { + Point_d point(coords.begin(), coords.end()); + tree.insert(point); + } + } + in_file.close(); +} +*/ + +void write_wl( std::string file_name, std::vector< std::vector > & WL) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : WL) + { + for (auto l: w) + ofs << l << " "; + ofs << "\n"; + } + ofs.close(); +} + + +std::vector 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 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 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& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + //std::vector landmarks; + int chosen_landmark; + //std::pair res = std::make_pair(L_i.begin(),false); + Point_d* p; + CGAL::Random rand; + for (int i = 0; i < nbL; i++) + { + // while (!res.second) + // { + while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0) + chosen_landmark = rand.get_int(0,nbP); + //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; + } + } + + +int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + //********************Preface: origin point + int D = W[0].size(); + std::vector orig_vector; + for (int i=0; i({0.8,0.8})), p2(std::vector({0.1,0.1})); + FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); + //std::cout << "Lambda=" << lambda << std::endl; + //FT lambda = 0.1;//Euclidean_distance(); + STraits traits(&(landmarks[0])); + std::vector< std::vector > WL(nbP); + Tree L(boost::counting_iterator(0), + boost::counting_iterator(nbL), + typename Tree::Splitter(), + traits); + /*Tree2 L2(boost::counting_iterator(0), + boost::counting_iterator(nbL), + typename Tree::Splitter(), + STraits(&(landmarks[0]))); + */ + 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::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+1, FT(0), true, + //CGAL::Distance_adapter(&(landmarks[0])) ); + CGAL::Distance_adapter(&(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; + } + if (i == landmarks_ind[WL[i][0]]) + { + //std::cout << "'"; + FT dist = ed.transformed_distance(W[i], landmarks[WL[i][1]]); + if (dist < lambda) + lambda = dist; + } + } + //std::cout << "\n"; + + std::string out_file = "wl_result"; + write_wl(out_file,WL); + + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.witness_complex(WL); + //******************** Making a set of bad link landmarks + std::cout << "Entered bad links\n"; + std::set< int > perturbL; + int count_badlinks = 0; + //std::cout << "Bad links around "; + for (auto u: witnessComplex.complex_vertex_range()) + if (!witnessComplex.has_good_link(u)) + { + //std::cout << "Landmark " << u << " start!" << std::endl; + //perturbL.insert(u); + count_badlinks++; + //std::cout << u << " "; + Point_d& l = landmarks[u]; + Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + std::vector curr_perturb; + L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); + for (int i: curr_perturb) + perturbL.insert(i%nbL); + //L.search(std::inserter(perturbL,perturbL.begin()),fs); + //L.search(std::ostream_iterator(std::cout,"\n"),fs); + //std::cout << "PerturbL size is " << perturbL.size() << std::endl; + } + std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; + //std::cout << "landmark[0][0] before" << landmarks[0][0] << std::endl; + //*********************** Perturb bad link landmarks + + for (auto u: perturbL) + { + Random_point_iterator rp(D,sqrt(lambda)/8); + //std::cout << landmarks[u] << std::endl; + + std::vector point; + for (int i = 0; i < D; i++) + { + while (K().squared_distance_d_object()(*rp,origin) < lambda/256) + rp++; + //FT coord = W[landmarks_ind[u]][i] + (*rp)[i]; + FT coord = landmarks[u][i] + (*rp)[i]; + if (coord > 1) + point.push_back(coord-1); + else if (coord < -1) + point.push_back(coord+1); + else + point.push_back(coord); + } + landmarks[u] = Point_d(point); + //std::cout << landmarks[u] << std::endl; + } + + //std::cout << "landmark[0][0] after" << landmarks[0][0] << std::endl; + std::cout << "lambda=" << lambda << std::endl; + + //std::cout << "WL size" << WL.size() << std::endl; + /* + std::cout << "L:" << std::endl; + for (int i = 0; i < landmarks.size(); i++) + std::cout << landmarks[i] << std::endl; + */ + + 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({0.1,0.1})), Point_d(std::vector({1.9,1.9}))) << std::endl; + return count_badlinks; +} + + +int main (int argc, char * const argv[]) +{ + + if (argc != 4) + { + std::cerr << "Usage: " << argv[0] + << " nbP nbL dim\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]); + //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); + // witnessComplex.witness_complex_from_points(point_vector); + //int nbP = point_vector.size(); + //std::vector > WL(nbP); + //std::set L; + Point_Vector L; + std::vector chosen_landmarks; + //Point_etiquette_map L_i; + //start = clock(); + //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); + bool ok=false; + while (!ok) + { + ok = true; + L = {}; + chosen_landmarks = {}; + landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); + for (auto i: chosen_landmarks) + { + ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); + if (!ok) break; + } + } + int bl = nbL, curr_min = bl; + write_points("landmarks/initial_pointset",point_vector); + write_points("landmarks/initial_landmarks",L); + + for (int i = 0; bl > 0; i++) + { + std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + bl=landmark_perturbation(point_vector, L, chosen_landmarks); + if (bl < curr_min) + curr_min=bl; + write_points("landmarks/landmarks0",L); + } + //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(); + */ +} -- cgit v1.2.3 From 3efbb8c26aecc9f88dae4ec69d3149a4f5f1bd40 Mon Sep 17 00:00:00 2001 From: skachano Date: Sat, 30 May 2015 15:36:35 +0000 Subject: witness: Testing continues git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@600 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 57d8d66a90e9839887289e4efd7fbddabe009d94 --- .../example/witness_complex_perturbations.cpp | 2 +- .../example/witness_complex_sphere.cpp | 15 ++++++++-- .../include/gudhi/Witness_complex.h | 32 ++++++++++++++++------ 3 files changed, 37 insertions(+), 12 deletions(-) (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp index 07b1eb8d..9d39906b 100644 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -328,7 +328,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< for (auto u: perturbL) { - Random_point_iterator rp(D,sqrt(lambda)/2); + Random_point_iterator rp(D,sqrt(lambda)/4); //std::cout << landmarks[u] << std::endl; std::vector point; diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index c9a9119a..7ba1dd61 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -489,8 +489,9 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, { // while (!res.second) // { - while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0) + 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]; @@ -546,7 +547,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< if (w[0]>0.95) std::cout << i << std::endl; */ - K_neighbor_search search(L, w, D+1, FT(0), true, + K_neighbor_search search(L, w, D, FT(0), true, //CGAL::Distance_adapter(&(landmarks[0])) ); CGAL::Distance_adapter(&(landmarks[0])) ); //std::cout << "Safely found nearest landmarks\n"; @@ -583,8 +584,10 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::set< int > perturbL; int count_badlinks = 0; //std::cout << "Bad links around "; + std::vector< int > count_bad(D+3); + std::vector< int > count_good(D+3); for (auto u: witnessComplex.complex_vertex_range()) - if (!witnessComplex.has_good_link(u)) + if (!witnessComplex.has_good_link(u, count_bad, count_good)) { //std::cout << "Landmark " << u << " start!" << std::endl; //perturbL.insert(u); @@ -600,6 +603,12 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //L.search(std::ostream_iterator(std::cout,"\n"),fs); //std::cout << "PerturbL size is " << perturbL.size() << std::endl; } + for (unsigned int i = 0; i != count_good.size(); 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++) + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; //std::cout << "landmark[0][0] before" << landmarks[0][0] << std::endl; //*********************** Perturb bad link landmarks diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index c8eaa3ef..c45d144d 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -628,7 +628,7 @@ private: /** \brief Returns true if the link is good */ - bool has_good_link(Vertex_handle v) + bool has_good_link(Vertex_handle v, std::vector< int >& bad_count, std::vector< int >& good_count) { std::vector< Vertex_handle > link_vertices; // Fill link_vertices @@ -642,14 +642,28 @@ private: print_vector(link_vertices); std::cout << "\n"; */ - print_vector(link_vertices); std::cout << "\n"; + //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 (d >= good_count.size()) + { + std::cout << "d=" << d << std::endl; + std::cout << "gc.size=" << good_count.size() << std::endl; + print_vector(link_vertices); + std::cout << std::endl; + } + */ + //assert(d < good_count.size()); + if (d == -1) bad_count[0]++; //std::cout << " dim " << d << "\n"; //Siblings* curr_sibl = root(); //std::cout << "Currently at vertex " - return (link_is_pseudomanifold(link_vertices,d)); + bool b= (link_is_pseudomanifold(link_vertices,d)); + if (d != -1) {if (b) good_count[d]++; else bad_count[d]++;} + return b; + } /** \brief Search and output links around vertices that are not pseudomanifolds @@ -737,9 +751,9 @@ private: if (d >= final_d) { final_d = d; - std::cout << d << " "; - print_vector(curr_simplex); - std::cout << std::endl; + //std::cout << d << " "; + //print_vector(curr_simplex); + // std::cout << std::endl; } } /* @@ -789,13 +803,15 @@ private: //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; - } + + }*/ count_bad[dimension]++; return false; } -- cgit v1.2.3 From 9e99b02fa2d91aced59c5fded3f1a2e07afe30aa Mon Sep 17 00:00:00 2001 From: skachano Date: Mon, 1 Jun 2015 09:00:41 +0000 Subject: Corrected bad-good links git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@601 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 004b949f2d74b803483535b126e76b84b27f4a7e --- .../example/witness_complex_flat_torus.cpp | 14 ++- .../example/witness_complex_sphere.cpp | 9 +- .../include/gudhi/Witness_complex.h | 110 ++++++++++----------- 3 files changed, 71 insertions(+), 62 deletions(-) (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index ff995a4f..14d6426d 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -601,10 +601,12 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::set< int > perturbL; int count_badlinks = 0; //std::cout << "Bad links around "; + std::vector< int > count_bad(D); + std::vector< int > count_good(D); for (auto u: witnessComplex.complex_vertex_range()) { - std::cout << "Vertex " << u << " "; - if (!witnessComplex.has_good_link(u)) + //std::cout << "Vertex " << u << " "; + if (!witnessComplex.has_good_link(u, count_bad, count_good)) { //std::cout << "Landmark " << u << " start!" << std::endl; //perturbL.insert(u); @@ -621,6 +623,12 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //std::cout << "PerturbL size is " << perturbL.size() << std::endl; } } + for (unsigned int i = 0; i != count_good.size(); 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++) + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; //std::cout << "landmark[0][0] before" << landmarks[0][0] << std::endl; //*********************** Perturb bad link landmarks @@ -744,7 +752,7 @@ int main (int argc, char * const argv[]) write_points("landmarks/initial_pointset",point_vector); //write_points("landmarks/initial_landmarks",L); - for (int i = 0; i < 1; i++) + for (int i = 0; bl > 0; i++) { std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; bl=landmark_perturbation(point_vector, L, chosen_landmarks); diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index 7ba1dd61..f6345411 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -584,8 +584,8 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::set< int > perturbL; int count_badlinks = 0; //std::cout << "Bad links around "; - std::vector< int > count_bad(D+3); - std::vector< int > count_good(D+3); + std::vector< int > count_bad(D); + std::vector< int > count_good(D); for (auto u: witnessComplex.complex_vertex_range()) if (!witnessComplex.has_good_link(u, count_bad, count_good)) { @@ -623,8 +623,8 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< { while (K().squared_distance_d_object()(*rp,origin) < lambda/256) rp++; - //FT coord = W[landmarks_ind[u]][i] + (*rp)[i]; - FT coord = landmarks[u][i] + (*rp)[i]; + FT coord = W[landmarks_ind[u]][i] + (*rp)[i]; + //FT coord = landmarks[u][i] + (*rp)[i]; if (coord > 1) point.push_back(coord-1); else if (coord < -1) @@ -723,6 +723,7 @@ int main (int argc, char * const argv[]) write_points("landmarks/initial_landmarks",L); for (int i = 0; bl > 0; i++) + //for (int i = 0; i < 1; i++) { std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; bl=landmark_perturbation(point_vector, L, chosen_landmarks); diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index c45d144d..0888c086 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -630,51 +630,41 @@ private: */ bool has_good_link(Vertex_handle v, std::vector< int >& bad_count, std::vector< int >& good_count) { - std::vector< Vertex_handle > link_vertices; - // Fill link_vertices + 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()) - link_vertices.push_back(u); + star_vertices.push_back(u); } - /* - print_vector(link_vertices); - std::cout << "\n"; - */ - //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); + typeVectorVertex init_simplex = {star_vertices[0]}; + int d = star_dim(star_vertices, star_vertices.begin()+1, 0, init_simplex) - 1; //link_dim = star_dim - 1 + assert(init_simplex.size() == 1); /* - if (d >= good_count.size()) - { - std::cout << "d=" << d << std::endl; - std::cout << "gc.size=" << good_count.size() << std::endl; - print_vector(link_vertices); - std::cout << std::endl; - } - */ - //assert(d < good_count.size()); + 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]++; - //std::cout << " dim " << d << "\n"; - //Siblings* curr_sibl = root(); - //std::cout << "Currently at vertex " - bool b= (link_is_pseudomanifold(link_vertices,d)); + bool b= (link_is_pseudomanifold(star_vertices,d)); if (d != -1) {if (b) good_count[d]++; else bad_count[d]++;} - return b; + return (d != -1 && b); } /** \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 << ": "; @@ -693,14 +683,9 @@ private: // 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++) { @@ -718,42 +703,43 @@ private: } std::cout << "not_connected = " << nc << std::endl; } - + */ private: std::vector count_good; std::vector count_bad; int nc; - int link_dim(std::vector< Vertex_handle >& link_vertices, + int star_dim(std::vector< Vertex_handle >& star_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"; + //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; //std::cout << "Current vertex is " << - for (it = curr_v; it != link_vertices.end(); ++it) + for (it = curr_v; it != star_vertices.end(); ++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); + 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 = link_dim(link_vertices, it+1, curr_d+1, curr_simplex); + int d = star_dim(star_vertices, it+1, curr_d+1, curr_simplex); if (d >= final_d) { final_d = d; //std::cout << d << " "; //print_vector(curr_simplex); - // std::cout << std::endl; + //std::cout << std::endl; } } /* @@ -781,19 +767,19 @@ private: * 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, + 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 empty_vector = {}; - add_vertices(link_vertices, - link_vertices.begin(), + typeVectorVertex init_vector = {}; + add_vertices(star_vertices, + star_vertices.begin()+1, adj_graph, d_map, f_map, - empty_vector, + init_vector, 0, dimension); //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; @@ -823,60 +809,74 @@ private: //return (boost::connected_components(adj_graph, &components[0]) == 1); } - void add_vertices(typeVectorVertex& link_vertices, + void add_vertices(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 dimension) + int link_dimension) { Simplex_handle sh; Vertex_t vert; typename typeVectorVertex::iterator it; - std::pair resPair; + //std::pair resPair; //typename Graph_map::iterator resPair; //Add vertices //std::cout << "Entered add vertices\n"; - for (it = curr_v; it != link_vertices.end(); ++it) + for (it = curr_v; it != star_vertices.end(); ++it) { - curr_simplex.push_back(*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 << ""; */ - sh = find(curr_simplex); + 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 == dimension) + 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); - resPair = d_map.emplace(sh,vert); + d_map.emplace(sh,vert); + for (Simplex_handle sh_b: boundary_simplex_range(sh)) + { + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh_b,vert); + } } else { - if (curr_d == dimension-1) + /* + if (curr_d == link_dimension-1) { vert = boost::add_vertex(adj_graph); resPair = f_map.emplace(sh,vert); } - add_vertices(link_vertices, + */ + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_vertices(star_vertices, it+1, adj_graph, d_map, f_map, curr_simplex, - curr_d+1, dimension); + curr_d+1, link_dimension); } } /* else std::cout << "\n"; */ - curr_simplex.pop_back(); + curr_simplex.pop_back(); //pop the vertex in question } } -- cgit v1.2.3 From 81d222a2222bb7a4e6304ed1ce708598a8e29dc6 Mon Sep 17 00:00:00 2001 From: skachano Date: Fri, 12 Jun 2015 09:44:31 +0000 Subject: Added grid witnesses for flat torus git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@610 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 34d0bfe869ee04bd0dfa832e0e36b245f33b062f --- src/Witness_complex/example/CMakeLists.txt | 3 + .../example/relaxed_witness_complex_sphere.cpp | 460 +++++++++++++++++++++ .../example/witness_complex_flat_torus.cpp | 23 +- .../example/witness_complex_perturbations.cpp | 18 +- .../example/witness_complex_sphere.cpp | 3 +- 5 files changed, 493 insertions(+), 14 deletions(-) create mode 100644 src/Witness_complex/example/relaxed_witness_complex_sphere.cpp (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 09167cbe..77f95c79 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -86,6 +86,9 @@ if(CGAL_FOUND) add_executable ( witness_complex_sphere witness_complex_sphere.cpp ) target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(witness_complex_sphere ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( relaxed_witness_complex_sphere relaxed_witness_complex_sphere.cpp ) + target_link_libraries(relaxed_witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(relaxed_witness_complex_sphere ${CMAKE_CURRENT_BINARY_DIR}/relaxed_witness_complex_sphere ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") 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..53380124 --- /dev/null +++ b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp @@ -0,0 +1,460 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +//#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Relaxed_witness_complex.h" +#include "gudhi/reader_utils.h" +//#include + +//#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d 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 Euclidean_distance; + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair 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 Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_incremental_neighbor_search> 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 Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +//typedef K::Equal_d Equal_d; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_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 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 > & 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 ::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 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 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 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& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + //std::vector landmarks; + int chosen_landmark; + //std::pair 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& landmarks_ind, FT alpha) +{ + //********************Preface: origin point + unsigned D = W[0].size(); + std::vector 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 > WL(nbP); + std::vector< std::vector< typename std::vector::iterator > > ope_limits(nbP); + Tree L(boost::counting_iterator(0), + boost::counting_iterator(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::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark + Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter(&(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::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(&(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({0.1,0.1})), Point_d(std::vector({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 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_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index 14d6426d..42bf5e7e 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -323,7 +323,20 @@ read_points_cust ( std::string file_name , Point_Vector & points) void generate_points_grid(Point_Vector& W, int width, int D) { - + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(0.01*(cell_i%width)); + cell_i /= width; + } + W.push_back(point); + } } void generate_points_random_box(Point_Vector& W, int nbP, int dim) @@ -490,8 +503,8 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, { // while (!res.second) // { - while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0) - chosen_landmark = rand.get_int(0,nbP); + do chosen_landmark = rand.get_int(0,nbP); + while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); //rand++; //std::cout << "Chose " << chosen_landmark << std::endl; p = &W[chosen_landmark]; @@ -717,7 +730,9 @@ int main (int argc, char * const argv[]) std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; //read_points_cust(file_name, point_vector); - generate_points_random_box(point_vector, nbP, dim); + //generate_points_random_box(point_vector, nbP, dim); + generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim); + nbP = (int)(pow((int)pow(nbP, 1.0/dim), dim)); /* for (auto &p: point_vector) { diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp index ec7a5abe..88a7510a 100644 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -194,7 +194,7 @@ void write_edges_gnuplot(std::string file_name, Witness_complex<>& witness_compl /** 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& landmarks_ind) { std::cout << "Enter landmark choice to kd tree\n"; @@ -216,12 +216,13 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, //std::cout << "Added landmark " << chosen_landmark << std::endl; } } -/* +*/ + void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) { std::cout << "Enter landmark choice to kd tree\n"; //std::vector landmarks; - int chosen_landmark; + int chosen_landmark = 0; //std::pair res = std::make_pair(L_i.begin(),false); Point_d* p; CGAL::Random rand; @@ -229,9 +230,9 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, { // while (!res.second) // { - while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0) - chosen_landmark = rand.get_int(0,nbP); - //rand++; + do chosen_landmark = rand.uniform_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); @@ -241,7 +242,7 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, //std::cout << "Added landmark " << chosen_landmark << std::endl; } } -*/ + int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind) { @@ -417,7 +418,8 @@ int main (int argc, char * const argv[]) } //write_points("landmarks/initial_pointset",point_vector); write_points("landmarks/initial_landmarks",L); - for (int i = 0; bl != 0; i++) + //for (int i = 0; bl != 0; i++) + for (int i = 0; i < 1; i++) { std::cout << "========== Start iteration " << i << " ========\n"; bl = landmark_perturbation(point_vector, L, chosen_landmarks); diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index f6345411..f2bb9819 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -489,8 +489,7 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, { // while (!res.second) // { - do - chosen_landmark = rand.get_int(0,nbP); + 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; -- cgit v1.2.3 From c88e3b93de22be92cc7027f4c14ea4294f8c366f Mon Sep 17 00:00:00 2001 From: skachano Date: Tue, 16 Jun 2015 09:36:00 +0000 Subject: Fixed the changing WL matrix bug in witness complex. Added is_witness_complex test git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@615 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 60917e75766b53d08240b510c0508e7d781b56d6 --- .../example/relaxed_witness_complex_sphere.cpp | 1 + .../example/witness_complex_sphere.cpp | 10 +- .../include/gudhi/Relaxed_witness_complex.h | 886 +++++++++++++++++++++ .../include/gudhi/Witness_complex.h | 93 ++- 4 files changed, 962 insertions(+), 28 deletions(-) create mode 100644 src/Witness_complex/include/gudhi/Relaxed_witness_complex.h (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp index 53380124..067321ce 100644 --- a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp +++ b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp @@ -36,6 +36,7 @@ //#include "gudhi/graph_simplicial_complex.h" #include "gudhi/Relaxed_witness_complex.h" #include "gudhi/reader_utils.h" +#include "gudhi/Collapse/Collapse.h" //#include //#include diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index f2bb9819..d08c763f 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -578,6 +578,10 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< Witness_complex<> witnessComplex; witnessComplex.setNbL(nbL); witnessComplex.witness_complex(WL); + if (witnessComplex.is_witness_complex(WL)) + std::cout << "!!YES. IT IS A WITNESS COMPLEX!!\n"; + else + std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n"; //******************** Making a set of bad link landmarks std::cout << "Entered bad links\n"; std::set< int > perturbL; @@ -614,7 +618,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< for (auto u: perturbL) { - Random_point_iterator rp(D,sqrt(lambda)/8); + Random_point_iterator rp(D,sqrt(lambda)/8*nbL/count_badlinks); //std::cout << landmarks[u] << std::endl; std::vector point; @@ -721,8 +725,8 @@ int main (int argc, char * const argv[]) write_points("landmarks/initial_pointset",point_vector); write_points("landmarks/initial_landmarks",L); - for (int i = 0; bl > 0; i++) - //for (int i = 0; i < 1; i++) + //for (int i = 0; bl > 0; i++) + for (int i = 0; i < 1; i++) { std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; bl=landmark_perturbation(point_vector, L, chosen_landmarks); 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 . + */ + +#ifndef GUDHI_RELAXED_WITNESS_COMPLEX_H_ +#define GUDHI_RELAXED_WITNESS_COMPLEX_H_ + +#include +#include +#include +#include +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include +#include +#include +#include +#include +#include +#include +#include + +// Needed for nearest neighbours +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +// Needed for the adjacency graph in bad link search +#include +#include +#include + +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 + 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 Node; + /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */ + typedef typename boost::container::flat_map 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 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::iterator curr_l, typename std::vector< std::vector::iterator >::iterator curr_until, std::vector& simplex, std::vector::iterator until, std::vector::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::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::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 > 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& 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::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + //facet = {}; + for (std::vector::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 + void print_vector(std::vector 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 + void print_vvector(std::vector< std::vector > 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 + 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()); // distance matrix witness x landmarks + typeVectorVertex chosen_landmarks; // landmark list + + WL = KNearestNeighbours(nbP,std::vector()); + 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::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 + // 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()); // distance matrix witness x landmarks + + // WL = KNearestNeighbours(nbP,std::vector()); + // 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 + void landmark_choice_by_random_points(Point_Vector &W, int nbP, std::set &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()); // distance matrix witness x landmarks + + //WL = KNearestNeighbours(nbP,std::vector()); + 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 + void nearest_landmarks(Point_Vector &W, std::set &L, KNearestNeighbours &WL) + { + int D = W[0].size(); + int nbP = W.size(); + WL = KNearestNeighbours(nbP,std::vector()); + typedef std::pair 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, comp> l_heap([&](dist_i j1, dist_i j2){return j1.first > j2.first;}); + std::set::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 count_good; + std::vector 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 Adj_graph; + // map that gives to a certain simplex its node in graph and its dimension + //typedef std::pair Reference; + typedef boost::graph_traits::vertex_descriptor Vertex_t; + typedef boost::graph_traits::edge_descriptor Edge_t; + + typedef boost::container::flat_map 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 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 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 index 04fcc98f..201d6525 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -140,7 +140,6 @@ namespace Gudhi { //void witness_complex(std::vector< std::vector< Vertex_handle > > & knn) { std::cout << "**Start the procedure witness_complex" << std::endl; - int k=2; /* current dimension in iterative construction */ //Construction of the active witness list int nbW = knn.size(); //int nbL = knn.at(0).size(); @@ -157,13 +156,14 @@ namespace Gudhi { // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore counter++; vv = {i}; - /* TODO Filtration */ 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 @@ -202,39 +202,37 @@ namespace Gudhi { active_w.push_back(i); } } - std::cout << "k=1, active witnesses: " << active_w.size() << std::endl; + */ + 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,0}; - count_bad = {0,0}; + 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); - count++; //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*/ - // First sort the first k landmarks - VertexHandle inserted_vertex = knn[*it][k]; - bool ok = all_faces_in(knn, *it, k, inserted_vertex); + 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); - if (returnValue.second) - count++; 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: " < - bool all_faces_in(KNearestNeighbours &knn, int witness_id, int k, VertexHandle inserted_vertex) + 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; @@ -332,20 +330,17 @@ private: // CHECK ALL THE FACETS for (int i = 0; i != k+1; ++i) { - if (knn[witness_id][i] != inserted_vertex) + facet = {}; + for (int j = 0; j != k+1; ++j) { - facet = {}; - for (int j = 0; j != k+1; ++j) + if (j != i) { - if (j != i) - { - facet.push_back(knn[witness_id][j]); - } - }//endfor - if (find(facet) == null_simplex()) - return false; - //std::cout << "++++ finished loop safely\n"; - }//endif + facet.push_back(knn[witness_id][j]); + } + }//endfor + if (find(facet) == null_simplex()) + return false; + //std::cout << "++++ finished loop safely\n"; } //endfor return true; } @@ -1059,6 +1054,54 @@ bool complex_is_pseudomanifold(int dimension) } } + 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 -- cgit v1.2.3 From d4c823dce7436f561c8cd545f5a4f80d2d301ca1 Mon Sep 17 00:00:00 2001 From: skachano Date: Sat, 4 Jul 2015 09:49:32 +0000 Subject: Added files for testing on didgeridoo git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@676 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 6d8aeda398803d49ee9218b99939095098ceb92e --- src/Witness_complex/example/CMakeLists.txt | 15 +- src/Witness_complex/example/Torus_distance.h | 209 ++++++++ .../example/witness_complex_cube.cpp | 541 +++++++++++++++++++ .../example/witness_complex_cubic_systems.cpp | 547 +++++++++++++++++++ .../example/witness_complex_epsilon.cpp | 566 ++++++++++++++++++++ .../example/witness_complex_flat_torus.cpp | 73 ++- .../example/witness_complex_perturbations.cpp | 2 +- .../example/witness_complex_protected_delaunay.cpp | 594 +++++++++++++++++++++ .../example/witness_complex_sphere.cpp | 248 ++------- 9 files changed, 2582 insertions(+), 213 deletions(-) create mode 100644 src/Witness_complex/example/Torus_distance.h create mode 100644 src/Witness_complex/example/witness_complex_cube.cpp create mode 100644 src/Witness_complex/example/witness_complex_cubic_systems.cpp create mode 100644 src/Witness_complex/example/witness_complex_epsilon.cpp create mode 100644 src/Witness_complex/example/witness_complex_protected_delaunay.cpp (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 77f95c79..23919b4a 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -87,8 +87,19 @@ if(CGAL_FOUND) target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(witness_complex_sphere ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) add_executable ( relaxed_witness_complex_sphere relaxed_witness_complex_sphere.cpp ) - target_link_libraries(relaxed_witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - add_test(relaxed_witness_complex_sphere ${CMAKE_CURRENT_BINARY_DIR}/relaxed_witness_complex_sphere ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_test(witness_complex_sphere ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( witness_complex_protected_delaunay witness_complex_protected_delaunay.cpp ) + target_link_libraries(witness_complex_protected_delaunay ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(witness_complex_protected_delaunay ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_protected_delaunay ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( witness_complex_cubic_systems witness_complex_cubic_systems.cpp ) + target_link_libraries(witness_complex_cubic_systems ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(witness_complex_cubic_systems ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_cubic_systems ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( witness_complex_cube witness_complex_cube.cpp ) + target_link_libraries(witness_complex_cube ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(witness_complex_cube ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_cube ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + add_executable ( witness_complex_epsilon witness_complex_epsilon.cpp ) + target_link_libraries(witness_complex_epsilon ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(witness_complex_epsilon ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_epsilon ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") endif() diff --git a/src/Witness_complex/example/Torus_distance.h b/src/Witness_complex/example/Torus_distance.h new file mode 100644 index 00000000..5ae127df --- /dev/null +++ b/src/Witness_complex/example/Torus_distance.h @@ -0,0 +1,209 @@ +#ifndef GUDHI_TORUS_DISTANCE_H_ +#define GUDHI_TORUS_DISTANCE_H_ + +#include + +#include +#include +#include + +typedef CGAL::Epick_d K; +typedef K::Point_d Point_d; +typedef K::FT FT; +typedef CGAL::Search_traits< + FT, Point_d, + typename K::Cartesian_const_iterator_d, + typename K::Construct_cartesian_const_iterator_d> Traits_base; + +/** + * \brief Class of distance in a flat torus in dimension D + * + */ +class Torus_distance { + +public: + typedef K::FT FT; + typedef K::Point_d Point_d; + typedef Point_d Query_item; + typedef typename CGAL::Dynamic_dimension_tag D; + + double box_length = 2; + + FT transformed_distance(Query_item q, Point_d p) const + { + FT distance = FT(0); + FT coord = FT(0); + //std::cout << "Hello skitty!\n"; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1), pit = construct_it(p); + for(; qit != qe; qit++, pit++) + { + coord = sqrt(((*qit)-(*pit))*((*qit)-(*pit))); + if (coord*coord <= (box_length-coord)*(box_length-coord)) + distance += coord*coord; + else + distance += (box_length-coord)*(box_length-coord); + } + return distance; + } + + FT min_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r) const { + FT distance = FT(0); + FT dist1, dist2; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if((*qit) < r.min_coord(i)) + { + dist1 = (r.min_coord(i)-(*qit)); + dist2 = (box_length - r.max_coord(i)+(*qit)); + if (dist1 < dist2) + distance += dist1*dist1; + else + distance += dist2*dist2; + } + else if ((*qit) > r.max_coord(i)) + { + dist1 = (box_length - (*qit)+r.min_coord(i)); + dist2 = ((*qit) - r.max_coord(i)); + if (dist1 < dist2) + distance += dist1*dist1; + else + distance += dist2*dist2; + } + } + return distance; + } + + FT min_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r, + std::vector& dists) const { + FT distance = FT(0); + FT dist1, dist2; + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + //std::cout << r.max_coord(0) << std::endl; + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if((*qit) < r.min_coord(i)) + { + dist1 = (r.min_coord(i)-(*qit)); + dist2 = (box_length - r.max_coord(i)+(*qit)); + if (dist1 < dist2) + { + dists[i] = dist1; + distance += dist1*dist1; + } + else + { + dists[i] = dist2; + distance += dist2*dist2; + //std::cout << "Good stuff1\n"; + } + } + else if ((*qit) > r.max_coord(i)) + { + dist1 = (box_length - (*qit)+r.min_coord(i)); + dist2 = ((*qit) - r.max_coord(i)); + if (dist1 < dist2) + { + dists[i] = dist1; + distance += dist1*dist1; + //std::cout << "Good stuff2\n"; + } + else + { + dists[i] = dist2; + distance += dist2*dist2; + } + } + }; + return distance; + } + + FT max_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r) const { + FT distance=FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if (box_length <= (r.min_coord(i)+r.max_coord(i))) + if ((r.max_coord(i)+r.min_coord(i)-box_length)/FT(2.0) <= (*qit) && + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + else + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + else + if ((box_length-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + else + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + return distance; + } + + + FT max_distance_to_rectangle(const Query_item& q, + const CGAL::Kd_tree_rectangle& r, + std::vector& dists) const { + FT distance=FT(0); + typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); + typename K::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++) + { + if (box_length <= (r.min_coord(i)+r.max_coord(i))) + if ((r.max_coord(i)+r.min_coord(i)-box_length)/FT(2.0) <= (*qit) && + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + { + dists[i] = r.max_coord(i)-(*qit); + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + } + else + { + dists[i] = sqrt(((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i))); + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + else + if ((box_length-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) + { + dists[i] = sqrt((r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit))); + distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); + + } + else + { + dists[i] = (*qit)-r.min_coord(i); + distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); + } + } + return distance; + } + + inline FT new_distance(FT dist, FT old_off, FT new_off, + int ) const { + + FT new_dist = dist + (new_off*new_off - old_off*old_off); + return new_dist; + } + + inline FT transformed_distance(FT d) const { + return d*d; + } + + inline FT inverse_of_transformed_distance(FT d) const { + return sqrt(d); + } + +}; + +#endif diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp new file mode 100644 index 00000000..7545f156 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_cube.cpp @@ -0,0 +1,541 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +//#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +#include "Torus_distance.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d K; +typedef K::Point_d Point_d; +//typedef CGAL::Cartesian_d K; +//typedef CGAL::Point_d Point_d; +typedef K::FT FT; +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 Euclidean_distance; + + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair 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 Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef boost::container::flat_map Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +//typedef K::Equal_d Equal_d; +//typedef CGAL::Random_points_in_cube_d > > Random_cube_iterator; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; + +typedef CGAL::Delaunay_triangulation Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef CGAL::Sphere_d Sphere_d; + +bool toric=false; + + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , Point_Vector & points) +{ + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector< double > point; + std::istringstream iss( line ); + while(iss >> x) { point.push_back(x); } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +void generate_points_grid(Point_Vector& W, int width, int D) +{ + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(0.01*(cell_i%width)); + cell_i /= width; + } + W.push_back(point); + } +} + +void generate_points_random_box(Point_Vector& W, int nbP, int dim) +{ + /* + Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + std::vector point; + for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it) + point.push_back(*it); + W.push_back(Point_d(point)); + rp++; + } + */ + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + + +void write_wl( std::string file_name, std::vector< std::vector > & 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_points( std::string file_name, std::vector< Point_d > & points) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : points) + { + for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + } + ofs.close(); +} + +void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto u: witness_complex.complex_vertex_range()) + for (auto v: witness_complex.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) + { + for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n\n\n"; + } + } + ofs.close(); +} + + +void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count) +{ + delaunay.insert(W[chosen_landmark]); + landmarks_ind.push_back(chosen_landmark); + landmark_count++; +} + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + c = t.locate(p, lt, f, ft, v); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + if (!t.is_infinite(fc_it)) + { + std::vector vertices; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT dist2 = ed.transformed_distance(center_cs, p); + //if the new point is inside the protection ball of a non conflicting simplex + if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) + return true; + } + return false; +} + +bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) +{ + Euclidean_distance ed; + int D = t.current_dimension(); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + if (!t.is_infinite(fc_it)) + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + { + //check if vertex belongs to the face + bool belongs = false; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + if (v_it == *fc_v_it) + { + belongs = true; + break; + } + if (!belongs) + { + std::vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT dist2 = ed.transformed_distance(center_cs, v_it->point()); + //if the new point is inside the protection ball of a non conflicting simplex + if (dist2 <= (r+delta)*(r+delta)) + return false; + } + } + return true; +} + +void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + for (int j = 0; j < landmarks_ind.size(); ++j) + landmarks.push_back(W[landmarks_ind[j]][l]); +} + +void landmark_choice_by_delaunay(Point_Vector& W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind, FT delta) +{ + int D = W[0].size(); + Delaunay_triangulation t(D); + CGAL::Random rand; + int chosen_landmark; + int landmark_count = 0; + for (int i = 0; i <= D+1; ++i) + { + do chosen_landmark = rand.get_int(0,nbP); + while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + insert_delaunay_landmark_with_copies(W, chosen_landmark, landmarks_ind, t, landmark_count); + } + while (landmark_count < nbL) + { + do chosen_landmark = rand.get_int(0,nbP); + while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + // If no conflicts then insert in every copy of T^3 + if (!is_violating_protection(W[chosen_landmark], t, D, delta)) + insert_delaunay_landmark_with_copies(W, chosen_landmark, landmarks_ind, t, landmark_count); + } +} + + +void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector& landmarks_ind, FT delta) +{ + int D = W[0].size(); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + CGAL::Random rand; + int landmark_count = 0; + std::list index_list; + // shuffle the list of indexes (via a vector) + { + std::vector temp_vector; + for (int i = 0; i < nbP; ++i) + temp_vector.push_back(i); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + for (std::vector::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + index_list.push_front(*it); + } + // add the first D+1 vertices to form one non-empty cell + for (int i = 0; i <= D+1; ++i) + { + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); + index_list.pop_front(); + } + // add other vertices if they don't violate protection + std::list::iterator list_it = index_list.begin(); + while (list_it != index_list.end()) + if (!is_violating_protection(W[*list_it], t, D, delta)) + { + // If no conflicts then insert in every copy of T^3 + insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + list_it++; + fill_landmark_copies(W, landmarks, landmarks_ind); +} + + +int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + //******************** Preface: origin point + int D = W[0].size(); + std::vector orig_vector; + for (int i=0; i landmarks_ext; + int nb_cells = 1; + for (int i = 0; i < D; ++i) + nb_cells *= 3; + for (int i = 0; i < nb_cells; ++i) + for (int k = 0; k < nbL; ++k) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(landmarks[k][l] + 2.0*((cell_i%3)-1.0)); + cell_i /= 3; + } + landmarks_ext.push_back(point); + } + write_points("landmarks/initial_landmarks",landmarks_ext); + STraits traits(&(landmarks_ext[0])); + std::vector< std::vector > WL(nbP); + + //********************** Neighbor search in a Kd tree + Tree L(boost::counting_iterator(0), + boost::counting_iterator(nb_cells*nbL), + typename Tree::Splitter(), + traits); + std::cout << "Enter (D+1) nearest landmarks\n"; + for (int i = 0; i < nbP; i++) + { + Point_d& w = W[i]; + ////Search D+1 nearest neighbours from the tree of landmarks L + K_neighbor_search search(L, w, D+1, FT(0), true, + CGAL::Distance_adapter(&(landmarks_ext[0])) ); + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) + { + if (std::find(WL[i].begin(), WL[i].end(), (it->first)%nbL) == WL[i].end()) + WL[i].push_back((it->first)%nbL); + } + if (i == landmarks_ind[WL[i][0]]) + { + FT dist = ed.transformed_distance(W[i], landmarks[WL[i][1]]); + if (dist < lambda) + lambda = dist; + } + } + std::string out_file = "wl_result"; + write_wl(out_file,WL); + + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.witness_complex(WL); + + //******************** Making a set of bad link landmarks + std::cout << "Entered bad links\n"; + std::set< int > perturbL; + int count_badlinks = 0; + //std::cout << "Bad links around "; + std::vector< int > count_bad(D); + std::vector< int > count_good(D); + for (auto u: witnessComplex.complex_vertex_range()) + { + if (!witnessComplex.has_good_link(u, count_bad, count_good)) + { + count_badlinks++; + Point_d& l = landmarks[u]; + Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + std::vector curr_perturb; + L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); + for (int i: curr_perturb) + perturbL.insert(i%nbL); + } + } + for (unsigned int i = 0; i != count_good.size(); 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++) + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; + + //*********************** Perturb bad link landmarks + for (auto u: perturbL) + { + Random_point_iterator rp(D,sqrt(lambda)/8); + std::vector point; + for (int i = 0; i < D; i++) + { + while (K().squared_distance_d_object()(*rp,origin) < lambda/256) + rp++; + FT coord = landmarks[u][i] + (*rp)[i]; + if (coord > 1) + point.push_back(coord-1); + else if (coord < -1) + point.push_back(coord+1); + else + point.push_back(coord); + } + landmarks[u] = Point_d(point); + } + std::cout << "lambda=" << lambda << std::endl; + 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); + return count_badlinks; +} + + +int main (int argc, char * const argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] + << " nbP nbL dim delta\n"; + return 0; + } + int nbP = atoi(argv[1]); + int nbL = atoi(argv[2]); + int dim = atoi(argv[3]); + FT delta = atof(argv[4]); + + std::cout << "Let the carnage begin!\n"; + Point_Vector point_vector; + generate_points_random_box(point_vector, nbP, dim); + Point_Vector L; + std::vector chosen_landmarks; + bool ok=false; + while (!ok) + { + ok = true; + L = {}; + chosen_landmarks = {}; + //landmark_choice_by_delaunay(point_vector, nbP, nbL, L, chosen_landmarks, delta); + landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta); + nbL = chosen_landmarks.size(); + std::cout << "Number of landmarks is " << nbL << std::endl; + //int width = (int)pow(nbL, 1.0/dim); landmark_choice_bcc(point_vector, nbP, width, L, chosen_landmarks); + for (auto i: chosen_landmarks) + { + ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); + if (!ok) break; + } + + } + int bl = nbL, curr_min = bl; + write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_landmarks",L); + //for (int i = 0; i < 1; i++) + for (int i = 0; bl > 0; i++) + { + std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + bl=landmark_perturbation(point_vector, nbL, L, chosen_landmarks); + if (bl < curr_min) + curr_min=bl; + write_points("landmarks/landmarks0",L); + } + +} diff --git a/src/Witness_complex/example/witness_complex_cubic_systems.cpp b/src/Witness_complex/example/witness_complex_cubic_systems.cpp new file mode 100644 index 00000000..2f4ee1cb --- /dev/null +++ b/src/Witness_complex/example/witness_complex_cubic_systems.cpp @@ -0,0 +1,547 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +#include "Torus_distance.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d K; +typedef K::Point_d Point_d; +//typedef CGAL::Cartesian_d K; +//typedef CGAL::Point_d Point_d; +typedef K::FT FT; +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 Euclidean_distance; + + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair 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 Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef boost::container::flat_map Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +//typedef K::Equal_d Equal_d; +//typedef CGAL::Random_points_in_cube_d > > Random_cube_iterator; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; + +typedef CGAL::Delaunay_triangulation Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef CGAL::Sphere_d Sphere_d; + +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_random_box(Point_Vector& W, int nbP, int dim) +{ + /* + Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + std::vector point; + for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it) + point.push_back(*it); + W.push_back(Point_d(point)); + rp++; + } + */ + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + + +void write_wl( std::string file_name, std::vector< std::vector > & 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_points( std::string file_name, std::vector< Point_d > & points) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : points) + { + for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + } + ofs.close(); +} + +void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto u: witness_complex.complex_vertex_range()) + for (auto v: witness_complex.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) + { + for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n\n\n"; + } + } + ofs.close(); +} + + +/** 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& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + int chosen_landmark; + 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::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + //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 aux_fill_grid(Point_Vector& W, int& width, Point_Vector& landmarks, std::vector& landmarks_ind, std::vector & curr_pattern) +{ + int D = W[0].size(); + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + if (curr_pattern[l]) + point.push_back(-1.0+(2.0/width)*(cell_i%width)+(1.0/width)); + else + point.push_back(-1.0+(2.0/width)*(cell_i%width)); + cell_i /= width; + } + landmarks.push_back(Point_d(point)); + landmarks_ind.push_back(0);//landmarks_ind.push_back(W.size()); + //std::cout << "Added point " << W.size() << std::endl;; + //W.push_back(Point_d(point)); + } +} + +void aux_put_halves(Point_Vector& W, int& width, Point_Vector& landmarks, std::vector& landmarks_ind, std::vector& curr_pattern, std::vector::iterator curr_pattern_it, std::vector::iterator bool_it, std::vector::iterator bool_end) +{ + if (curr_pattern_it != curr_pattern.end()) + { + if (bool_it != bool_end) + { + *curr_pattern_it = false; + aux_put_halves(W, width, landmarks, landmarks_ind, curr_pattern, curr_pattern_it+1, bool_it, bool_end); + *curr_pattern_it = true; + aux_put_halves(W, width, landmarks, landmarks_ind, curr_pattern, curr_pattern_it+1, bool_it+1, bool_end); + } + } + else + if (*bool_it) + { + std::cout << "Filling the pattern "; + for (bool b: curr_pattern) + if (b) std::cout << '1'; + else std::cout << '0'; + std::cout << "\n"; + aux_fill_grid(W, width, landmarks, landmarks_ind, curr_pattern); + } +} + +void landmark_choice_cs(Point_Vector& W, int width, Point_Vector& landmarks, std::vector& landmarks_ind, std::vector& face_centers) +{ + std::cout << "Enter landmark choice to kd tree\n"; + //int chosen_landmark; + CGAL::Random rand; + //To speed things up check the last true in the code and put it as the finishing condition + unsigned last_true = face_centers.size()-1; + while (!face_centers[last_true] && last_true != 0) + last_true--; + //Recursive procedure to understand where we put +1/2 in centers' coordinates + std::vector curr_pattern(W[0].size(), false); + aux_put_halves(W, width, landmarks, landmarks_ind, curr_pattern, curr_pattern.begin(), face_centers.begin(), face_centers.begin()+(last_true+1)); + std::cout << "The number of landmarks is: " << landmarks.size() << std::endl; + + } + +int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + //******************** Preface: origin point + int D = W[0].size(); + std::vector orig_vector; + for (int i=0; i landmarks_ext; + int nb_cells = 1; + for (int i = 0; i < D; ++i) + nb_cells *= 3; + for (int i = 0; i < nb_cells; ++i) + for (int k = 0; k < nbL; ++k) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(landmarks[k][l] + 2.0*((cell_i%3)-1.0)); + cell_i /= 3; + } + landmarks_ext.push_back(point); + } + write_points("landmarks/initial_landmarks",landmarks_ext); + STraits traits(&(landmarks_ext[0])); + std::vector< std::vector > WL(nbP); + + //********************** Neighbor search in a Kd tree + Tree L(boost::counting_iterator(0), + boost::counting_iterator(nb_cells*nbL), + typename Tree::Splitter(), + traits); + std::cout << "Enter (D+1) nearest landmarks\n"; + for (int i = 0; i < nbP; i++) + { + Point_d& w = W[i]; + ////Search D+1 nearest neighbours from the tree of landmarks L + K_neighbor_search search(L, w, D+1, FT(0), true, + CGAL::Distance_adapter(&(landmarks_ext[0])) ); + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) + { + if (std::find(WL[i].begin(), WL[i].end(), (it->first)%nbL) == WL[i].end()) + WL[i].push_back((it->first)%nbL); + } + if (i == landmarks_ind[WL[i][0]]) + { + FT dist = ed.transformed_distance(W[i], landmarks[WL[i][1]]); + if (dist < lambda) + lambda = dist; + } + } + std::string out_file = "wl_result"; + write_wl(out_file,WL); + + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.witness_complex(WL); + + //******************** Making a set of bad link landmarks + std::cout << "Entered bad links\n"; + std::set< int > perturbL; + int count_badlinks = 0; + //std::cout << "Bad links around "; + std::vector< int > count_bad(D); + std::vector< int > count_good(D); + for (auto u: witnessComplex.complex_vertex_range()) + { + if (!witnessComplex.has_good_link(u, count_bad, count_good, D)) + { + count_badlinks++; + Point_d& l = landmarks[u]; + Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + std::vector curr_perturb; + L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); + for (int i: curr_perturb) + perturbL.insert(i%nbL); + } + } + for (unsigned int i = 0; i != count_good.size(); 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++) + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; + + //*********************** Perturb bad link landmarks + for (auto u: perturbL) + { + Random_point_iterator rp(D,sqrt(lambda)/8); + std::vector point; + for (int i = 0; i < D; i++) + { + while (K().squared_distance_d_object()(*rp,origin) < lambda/256) + rp++; + FT coord = landmarks[u][i] + (*rp)[i]; + if (coord > 1) + point.push_back(coord-1); + else if (coord < -1) + point.push_back(coord+1); + else + point.push_back(coord); + } + landmarks[u] = Point_d(point); + } + std::cout << "lambda=" << lambda << std::endl; + 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); + return count_badlinks; +} + +void exaustive_search(Point_Vector& W, int width) +{ + int D = W[0].size()+1; + int nb_points = pow(2,D); + std::vector face_centers(D, false); + int bl = 0; //Bad links + std::vector> good_patterns; + for (int i = 0; i < nb_points; ++i) + { + int cell_i = i; + for (int l = 0; l < D; ++l) + { + if (cell_i%2 == 0) + face_centers[l] = false; + else + face_centers[l] = true; + cell_i /= 2; + } + std::cout << "**Current pattern "; + for (bool b: face_centers) + if (b) std::cout << '1'; + else std::cout << '0'; + std::cout << "\n"; + Point_Vector landmarks; + std::vector landmarks_ind; + Point_Vector W_copy(W); + landmark_choice_cs(W_copy, width, landmarks, landmarks_ind, face_centers); + if (landmarks.size() != 0) + { + bl = landmark_perturbation(W_copy, landmarks, landmarks_ind); + if ((1.0*bl)/landmarks.size() < 0.5) + good_patterns.push_back(face_centers); + } + } + std::cout << "The following patterns worked: "; + for (std::vector pattern : good_patterns) + { + std::cout << "["; + for (bool b: pattern) + if (b) std::cout << '1'; + else std::cout << '0'; + std::cout << "] "; + } + std::cout << "\n"; +} + +int main (int argc, char * const argv[]) +{ + unsigned nbP = atoi(argv[1]); + unsigned width = atoi(argv[2]); + unsigned dim = atoi(argv[3]); + std::string code = (std::string) argv[4]; + bool e_option = false; + int c; + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] + << "witness_complex_cubic_systems nbP width dim code || witness_complex_systems -e nbP width dim\n" + << "where nbP stands for the number of witnesses, width for the width of the grid, dim for dimension " + << "and code is a sequence of (dim+1) symbols 0 and 1 representing if we take the centers of k-dimensional faces of the cubic system depending if it is 0 or 1." + << "-e stands for the 'exaustive' option"; + return 0; + } + while ((c = getopt (argc, argv, "e::")) != -1) + switch(c) + { + case 'e' : + e_option = true; + nbP = atoi(argv[2]); + width = atoi(argv[3]); + dim = atoi(argv[4]); + break; + default : + nbP = atoi(argv[1]); + width = atoi(argv[2]); + dim = atoi(argv[3]); + code = (std::string) argv[4]; + } + Point_Vector point_vector; + generate_points_random_box(point_vector, nbP, dim); + + // Exaustive search + if (e_option) + { + std::cout << "Start exaustive search!\n"; + exaustive_search(point_vector, width); + return 0; + } + // Search with a specific cubic system + std::vector face_centers; + if (code.size() != dim+1) + { + std::cerr << "The code should contain (dim+1) symbols"; + return 1; + } + for (char c: code) + if (c == '0') + face_centers.push_back(false); + else + face_centers.push_back(true); + std::cout << "Let the carnage begin!\n"; + Point_Vector L; + std::vector chosen_landmarks; + + landmark_choice_cs(point_vector, width, L, chosen_landmarks, face_centers); + + int nbL = width; //!!!!!!!!!!!!! + int bl = nbL, curr_min = bl; + write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_landmarks",L); + //for (int i = 0; i < 1; i++) + for (int i = 0; bl > 0; i++) + { + std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + bl=landmark_perturbation(point_vector, L, chosen_landmarks); + if (bl < curr_min) + curr_min=bl; + write_points("landmarks/landmarks0",L); + } + +} diff --git a/src/Witness_complex/example/witness_complex_epsilon.cpp b/src/Witness_complex/example/witness_complex_epsilon.cpp new file mode 100644 index 00000000..d091bdb7 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_epsilon.cpp @@ -0,0 +1,566 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +//#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +#include "Torus_distance.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d K; +typedef K::Point_d Point_d; +//typedef CGAL::Cartesian_d K; +//typedef CGAL::Point_d Point_d; +typedef K::FT FT; +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 Euclidean_distance; + + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair 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 Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef boost::container::flat_map Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +//typedef K::Equal_d Equal_d; +//typedef CGAL::Random_points_in_cube_d > > Random_cube_iterator; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; + +typedef CGAL::Delaunay_triangulation Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef CGAL::Sphere_d Sphere_d; + +bool toric=false; + + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , Point_Vector & points) +{ + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector< double > point; + std::istringstream iss( line ); + while(iss >> x) { point.push_back(x); } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +void generate_points_grid(Point_Vector& W, int width, int D) +{ + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(0.01*(cell_i%width)); + cell_i /= width; + } + W.push_back(point); + } +} + +void generate_points_random_box(Point_Vector& W, int nbP, int dim) +{ + /* + Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + std::vector point; + for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it) + point.push_back(*it); + W.push_back(Point_d(point)); + rp++; + } + */ + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + + +void write_wl( std::string file_name, std::vector< std::vector > & 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_points( std::string file_name, std::vector< Point_d > & points) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : points) + { + for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + } + ofs.close(); +} + +void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto u: witness_complex.complex_vertex_range()) + for (auto v: witness_complex.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) + { + for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n\n\n"; + } + } + ofs.close(); +} + + +/** 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& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + int chosen_landmark; + 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::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + //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 insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count) +{ + delaunay.insert(W[chosen_landmark]); + landmarks_ind.push_back(chosen_landmark); + landmark_count++; +} + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + c = t.locate(p, lt, f, ft, v); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + if (!t.is_infinite(fc_it)) + { + std::vector vertices; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT dist2 = ed.transformed_distance(center_cs, p); + //if the new point is inside the protection ball of a non conflicting simplex + if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) + return true; + } + return false; +} + +bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) +{ + Euclidean_distance ed; + int D = t.current_dimension(); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + if (!t.is_infinite(fc_it)) + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + { + //check if vertex belongs to the face + bool belongs = false; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + if (v_it == *fc_v_it) + { + belongs = true; + break; + } + if (!belongs) + { + std::vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT dist2 = ed.transformed_distance(center_cs, v_it->point()); + //if the new point is inside the protection ball of a non conflicting simplex + if (dist2 <= (r+delta)*(r+delta)) + return false; + } + } + return true; +} + +FT sampling_radius(Delaunay_triangulation& t) +{ + int D = t.current_dimension(); + FT epsilon2 = 4.0; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin())); + if (epsilon2 > r2) + epsilon2 = r2; + } + return sqrt(epsilon2); +} + +FT point_sampling_radius_by_delaunay(Point_Vector& points) +{ + Delaunay_triangulation t(points[0].size()); + t.insert(points.begin(), points.end()); + return sampling_radius(t); +} + +void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector& landmarks_ind, FT delta) +{ + int D = W[0].size(); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + CGAL::Random rand; + int landmark_count = 0; + std::list index_list; + // shuffle the list of indexes (via a vector) + { + std::vector temp_vector; + for (int i = 0; i < nbP; ++i) + temp_vector.push_back(i); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + for (std::vector::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + index_list.push_front(*it); + } + // add the first D+1 vertices to form one non-empty cell + for (int i = 0; i <= D+1; ++i) + { + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); + index_list.pop_front(); + } + // add other vertices if they don't violate protection + std::list::iterator list_it = index_list.begin(); + while (list_it != index_list.end()) + if (!is_violating_protection(W[*list_it], t, D, delta)) + { + // If no conflicts then insert in every copy of T^3 + insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + list_it++; + for (std::vector::iterator it = landmarks_ind.begin(); it != landmarks_ind.end(); ++it) + landmarks.push_back(W[*it]); +} + + +int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + //******************** Preface: origin point + int D = W[0].size(); + std::vector orig_vector; + for (int i=0; i landmarks_ext; + int nb_cells = 1; + for (int i = 0; i < D; ++i) + nb_cells *= 3; + for (int i = 0; i < nb_cells; ++i) + for (int k = 0; k < nbL; ++k) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(landmarks[k][l] + 2.0*((cell_i%3)-1.0)); + cell_i /= 3; + } + landmarks_ext.push_back(point); + } + write_points("landmarks/initial_landmarks",landmarks_ext); + STraits traits(&(landmarks_ext[0])); + std::vector< std::vector > WL(nbP); + + //********************** Neighbor search in a Kd tree + Tree L(boost::counting_iterator(0), + boost::counting_iterator(nb_cells*nbL), + typename Tree::Splitter(), + traits); + std::cout << "Enter (D+1) nearest landmarks\n"; + for (int i = 0; i < nbP; i++) + { + Point_d& w = W[i]; + ////Search D+1 nearest neighbours from the tree of landmarks L + K_neighbor_search search(L, w, D+1, FT(0), true, + CGAL::Distance_adapter(&(landmarks_ext[0])) ); + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) + { + if (std::find(WL[i].begin(), WL[i].end(), (it->first)%nbL) == WL[i].end()) + WL[i].push_back((it->first)%nbL); + } + if (i == landmarks_ind[WL[i][0]]) + { + FT dist = ed.transformed_distance(W[i], landmarks[WL[i][1]]); + if (dist < lambda) + lambda = dist; + } + } + std::string out_file = "wl_result"; + write_wl(out_file,WL); + + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.witness_complex(WL); + + //******************** Making a set of bad link landmarks + std::cout << "Entered bad links\n"; + std::set< int > perturbL; + int count_badlinks = 0; + //std::cout << "Bad links around "; + std::vector< int > count_bad(D); + std::vector< int > count_good(D); + for (auto u: witnessComplex.complex_vertex_range()) + { + if (!witnessComplex.has_good_link(u, count_bad, count_good)) + { + count_badlinks++; + Point_d& l = landmarks[u]; + Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + std::vector curr_perturb; + L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); + for (int i: curr_perturb) + perturbL.insert(i%nbL); + } + } + for (unsigned int i = 0; i != count_good.size(); 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++) + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; + + //*********************** Perturb bad link landmarks + for (auto u: perturbL) + { + Random_point_iterator rp(D,sqrt(lambda)/8); + std::vector point; + for (int i = 0; i < D; i++) + { + while (K().squared_distance_d_object()(*rp,origin) < lambda/256) + rp++; + FT coord = landmarks[u][i] + (*rp)[i]; + if (coord > 1) + point.push_back(coord-1); + else if (coord < -1) + point.push_back(coord+1); + else + point.push_back(coord); + } + landmarks[u] = Point_d(point); + } + std::cout << "lambda=" << lambda << std::endl; + 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); + return count_badlinks; +} + + +int main (int argc, char * const argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] + << " nbP nbL dim delta\n"; + return 0; + } + int nbP = atoi(argv[1]); + int nbL = atoi(argv[2]); + int dim = atoi(argv[3]); + FT delta = atof(argv[4]); + + std::cout << "Let the carnage begin!\n"; + Point_Vector point_vector; + generate_points_random_box(point_vector, nbP, dim); + FT epsilon = point_sampling_radius_by_delaunay(point_vector); + std::cout << "Initial epsilon = " << epsilon << std::endl; + Point_Vector L; + std::vector chosen_landmarks; + bool ok=false; + while (!ok) + { + ok = true; + L = {}; + chosen_landmarks = {}; + //landmark_choice_by_delaunay(point_vector, nbP, nbL, L, chosen_landmarks, delta); + landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta); + nbL = chosen_landmarks.size(); + std::cout << "Number of landmarks is " << nbL << std::endl; + //int width = (int)pow(nbL, 1.0/dim); landmark_choice_bcc(point_vector, nbP, width, L, chosen_landmarks); + for (auto i: chosen_landmarks) + { + ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); + if (!ok) break; + } + + } + FT epsilon2 = point_sampling_radius_by_delaunay(L); + std::cout << "Final epsilon = " << epsilon2 << ". Ratio = " << epsilon/epsilon2 << std::endl; + int bl = nbL, curr_min = bl; + write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_landmarks",L); + //for (int i = 0; i < 1; i++) + for (int i = 0; bl > 0; i++) + { + std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + bl=landmark_perturbation(point_vector, nbL, L, chosen_landmarks); + if (bl < curr_min) + curr_min=bl; + write_points("landmarks/landmarks0",L); + } + +} diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index 42bf5e7e..06bf5a9f 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -63,8 +63,10 @@ using namespace Gudhi; //using namespace boost::filesystem; typedef CGAL::Epick_d K; -typedef K::FT FT; typedef K::Point_d Point_d; +//typedef CGAL::Cartesian_d K; +//typedef CGAL::Point_d Point_d; +typedef K::FT FT; typedef CGAL::Search_traits< FT, Point_d, typename K::Cartesian_const_iterator_d, @@ -72,7 +74,7 @@ typedef CGAL::Search_traits< typedef CGAL::Euclidean_distance Euclidean_distance; /** - * \brief Class of distance in a flat torus in dimaension D + * \brief Class of distance in a flat torus in dimension D * */ //class Torus_distance : public Euclidean_distance { @@ -288,10 +290,11 @@ typedef CGAL::Fuzzy_sphere Fuzzy_sphere; typedef std::vector Point_Vector; //typedef K::Equal_d Equal_d; +//typedef CGAL::Random_points_in_cube_d > > Random_cube_iterator; typedef CGAL::Random_points_in_cube_d Random_cube_iterator; typedef CGAL::Random_points_in_ball_d Random_point_iterator; -bool toric=true; +bool toric=false; /** * \brief Customized version of read_points @@ -341,7 +344,18 @@ void generate_points_grid(Point_Vector& W, int width, int D) void generate_points_random_box(Point_Vector& W, int nbP, int dim) { + /* Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + std::vector point; + for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it) + point.push_back(*it); + W.push_back(Point_d(point)); + rp++; + } + */ + Random_cube_iterator rp(dim, 1.0); for (int i = 0; i < nbP; i++) { W.push_back(*rp++); @@ -494,9 +508,7 @@ void write_edges(std::string file_name, Witness_complex<>& witness_complex, Poin void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) { std::cout << "Enter landmark choice to kd tree\n"; - //std::vector landmarks; int chosen_landmark; - //std::pair res = std::make_pair(L_i.begin(),false); Point_d* p; CGAL::Random rand; for (int i = 0; i < nbL; i++) @@ -516,6 +528,33 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, } } +/** \brief Choose landmarks on a body-central cubic system + */ +void landmark_choice_bcc(Point_Vector &W, int nbP, int width, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + int D = W[0].size(); + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) + { + std::vector point; + std::vector cpoint; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(-1.0+(2.0/width)*(cell_i%width)); + cpoint.push_back(-1.0+(2.0/width)*(cell_i%width)+(1.0/width)); + cell_i /= width; + } + landmarks.push_back(point); + landmarks.push_back(cpoint); + landmarks_ind.push_back(2*i); + landmarks_ind.push_back(2*i+1); + } + std::cout << "The number of landmarks is: " << landmarks.size() << std::endl; +} + int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind) { @@ -548,7 +587,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< int cell_i = i; for (int l = 0; l < D; ++l) { - point.push_back(landmarks[k][l] + 2.0*((cell_i%3)-1)); + point.push_back(landmarks[k][l] + 2.0*((cell_i%3)-1.0)); cell_i /= 3; } landmarks_ext.push_back(point); @@ -587,7 +626,8 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //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)%nbL); + if (std::find(WL[i].begin(), WL[i].end(), (it->first)%nbL) == WL[i].end()) + WL[i].push_back((it->first)%nbL); //std::cout << "ITFIRST " << it->first << std::endl; //std::cout << i << " " << it->first << ": " << it->second << std::endl; } @@ -609,6 +649,12 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< Witness_complex<> witnessComplex; witnessComplex.setNbL(nbL); witnessComplex.witness_complex(WL); + /* + if (witnessComplex.is_witness_complex(WL)) + std::cout << "!!YES. IT IS A WITNESS COMPLEX!!\n"; + else + std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n"; + */ //******************** Making a set of bad link landmarks std::cout << "Entered bad links\n"; std::set< int > perturbL; @@ -730,9 +776,9 @@ int main (int argc, char * const argv[]) std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; //read_points_cust(file_name, point_vector); - //generate_points_random_box(point_vector, nbP, dim); - generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim); - nbP = (int)(pow((int)pow(nbP, 1.0/dim), dim)); + generate_points_random_box(point_vector, nbP, dim); + //generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim); + //nbP = (int)(pow((int)pow(nbP, 1.0/dim), dim)); /* for (auto &p: point_vector) { @@ -757,17 +803,20 @@ int main (int argc, char * const argv[]) L = {}; chosen_landmarks = {}; landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); + + //int width = (int)pow(nbL, 1.0/dim); landmark_choice_bcc(point_vector, nbP, width, L, chosen_landmarks); for (auto i: chosen_landmarks) { ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); if (!ok) break; } + } int bl = nbL, curr_min = bl; write_points("landmarks/initial_pointset",point_vector); //write_points("landmarks/initial_landmarks",L); - - for (int i = 0; bl > 0; i++) + for (int i = 0; i < 1; i++) + //for (int i = 0; bl > 0; i++) { std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; bl=landmark_perturbation(point_vector, L, chosen_landmarks); diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp index 88a7510a..b3b84b1f 100644 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -416,7 +416,7 @@ int main (int argc, char * const argv[]) { file_name.erase(0, last_slash_idx + 1); } - //write_points("landmarks/initial_pointset",point_vector); + write_points("landmarks/initial_pointset",point_vector); write_points("landmarks/initial_landmarks",L); //for (int i = 0; bl != 0; i++) for (int i = 0; i < 1; i++) diff --git a/src/Witness_complex/example/witness_complex_protected_delaunay.cpp b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp new file mode 100644 index 00000000..2f795a5f --- /dev/null +++ b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp @@ -0,0 +1,594 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +//#include + +//#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Witness_complex.h" +#include "gudhi/reader_utils.h" +#include "Torus_distance.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +using namespace Gudhi; +//using namespace boost::filesystem; + +typedef CGAL::Epick_d K; +typedef K::Point_d Point_d; +//typedef CGAL::Cartesian_d K; +//typedef CGAL::Point_d Point_d; +typedef K::FT FT; +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 Euclidean_distance; + + +typedef std::vector< Vertex_handle > typeVectorVertex; + +//typedef std::pair 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 Euclidean_adapter; +//typedef CGAL::Kd_tree Kd_tree; +typedef CGAL::Orthogonal_k_neighbor_search> K_neighbor_search; +typedef K_neighbor_search::Tree Tree; +typedef K_neighbor_search::Distance Distance; +typedef K_neighbor_search::iterator KNS_iterator; +typedef K_neighbor_search::iterator KNS_range; +typedef boost::container::flat_map Point_etiquette_map; +typedef CGAL::Kd_tree Tree2; + +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; + +typedef std::vector Point_Vector; + +//typedef K::Equal_d Equal_d; +//typedef CGAL::Random_points_in_cube_d > > Random_cube_iterator; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; + +typedef CGAL::Delaunay_triangulation Delaunay_triangulation; +typedef Delaunay_triangulation::Facet Facet; +typedef CGAL::Sphere_d Sphere_d; + +bool toric=false; + + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust ( std::string file_name , Point_Vector & points) +{ + std::ifstream in_file (file_name.c_str(),std::ios::in); + if(!in_file.is_open()) + { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while( getline ( in_file , line ) ) + { + std::vector< double > point; + std::istringstream iss( line ); + while(iss >> x) { point.push_back(x); } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +void generate_points_grid(Point_Vector& W, int width, int D) +{ + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(0.01*(cell_i%width)); + cell_i /= width; + } + W.push_back(point); + } +} + +void generate_points_random_box(Point_Vector& W, int nbP, int dim) +{ + /* + Random_cube_iterator rp(dim, 1); + for (int i = 0; i < nbP; i++) + { + std::vector point; + for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it) + point.push_back(*it); + W.push_back(Point_d(point)); + rp++; + } + */ + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) + { + W.push_back(*rp++); + } +} + + +void write_wl( std::string file_name, std::vector< std::vector > & 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_points( std::string file_name, std::vector< Point_d > & points) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : points) + { + for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + } + ofs.close(); +} + +void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto u: witness_complex.complex_vertex_range()) + for (auto v: witness_complex.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) + { + for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n\n\n"; + } + } + ofs.close(); +} + + +/** 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& landmarks_ind) +{ + std::cout << "Enter landmark choice to kd tree\n"; + int chosen_landmark; + 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::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + //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 insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count) +{ + int D = W[0].size(); + int nb_cells = pow(3, D); + for (int i = 0; i < nb_cells; ++i) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1)); + cell_i /= 3; + } + delaunay.insert(point); + } + landmarks_ind.push_back(chosen_landmark); + landmark_count++; +} + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) +{ + Euclidean_distance ed; + Delaunay_triangulation::Vertex_handle v; + Delaunay_triangulation::Face f(t.current_dimension()); + Delaunay_triangulation::Facet ft; + Delaunay_triangulation::Full_cell_handle c; + Delaunay_triangulation::Locate_type lt; + c = t.locate(p, lt, f, ft, v); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + if (!t.is_infinite(fc_it)) + { + std::vector vertices; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT dist2 = ed.transformed_distance(center_cs, p); + //if the new point is inside the protection ball of a non conflicting simplex + if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) + return true; + } + return false; +} + +bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) +{ + Euclidean_distance ed; + int D = t.current_dimension(); + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + if (!t.is_infinite(fc_it)) + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + { + //check if vertex belongs to the face + bool belongs = false; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + if (v_it == *fc_v_it) + { + belongs = true; + break; + } + if (!belongs) + { + std::vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs(D, vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT dist2 = ed.transformed_distance(center_cs, v_it->point()); + //if the new point is inside the protection ball of a non conflicting simplex + if (dist2 <= (r+delta)*(r+delta)) + return false; + } + } + return true; +} + +void fill_landmark_copies(Point_Vector& W, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + int D = W[0].size(); + int nb_cells = pow(3, D); + int nbL = landmarks_ind.size(); + // Fill landmarks + for (int i = 0; i < nb_cells-1; ++i) + for (int j = 0; j < nbL; ++j) + { + int cell_i = i; + Point_d point; + for (int l = 0; l < D; ++l) + { + point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1)); + cell_i /= 3; + } + landmarks.push_back(point); + } +} + +void landmark_choice_by_delaunay(Point_Vector& W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind, FT delta) +{ + int D = W[0].size(); + Delaunay_triangulation t(D); + CGAL::Random rand; + int chosen_landmark; + int landmark_count = 0; + for (int i = 0; i <= D+1; ++i) + { + do chosen_landmark = rand.get_int(0,nbP); + while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + insert_delaunay_landmark_with_copies(W, chosen_landmark, landmarks_ind, t, landmark_count); + } + while (landmark_count < nbL) + { + do chosen_landmark = rand.get_int(0,nbP); + while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + // If no conflicts then insert in every copy of T^3 + if (!is_violating_protection(W[chosen_landmark], t, D, delta)) + insert_delaunay_landmark_with_copies(W, chosen_landmark, landmarks_ind, t, landmark_count); + } + fill_landmark_copies(W, landmarks, landmarks_ind); +} + + +void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector& landmarks_ind, FT delta) +{ + int D = W[0].size(); + Torus_distance td; + Euclidean_distance ed; + Delaunay_triangulation t(D); + CGAL::Random rand; + int landmark_count = 0; + std::list index_list; + // shuffle the list of indexes (via a vector) + { + std::vector temp_vector; + for (int i = 0; i < nbP; ++i) + temp_vector.push_back(i); + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + for (std::vector::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) + index_list.push_front(*it); + } + // add the first D+1 vertices to form one non-empty cell + for (int i = 0; i <= D+1; ++i) + { + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); + index_list.pop_front(); + } + // add other vertices if they don't violate protection + std::list::iterator list_it = index_list.begin(); + while (list_it != index_list.end()) + if (!is_violating_protection(W[*list_it], t, D, delta)) + { + // If no conflicts then insert in every copy of T^3 + insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); + index_list.erase(list_it); + list_it = index_list.begin(); + } + else + list_it++; + fill_landmark_copies(W, landmarks, landmarks_ind); +} + + +int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + //******************** Preface: origin point + int D = W[0].size(); + std::vector orig_vector; + for (int i=0; i landmarks_ext; + int nb_cells = 1; + for (int i = 0; i < D; ++i) + nb_cells *= 3; + for (int i = 0; i < nb_cells; ++i) + for (int k = 0; k < nbL; ++k) + { + std::vector point; + int cell_i = i; + for (int l = 0; l < D; ++l) + { + point.push_back(landmarks[k][l] + 2.0*((cell_i%3)-1.0)); + cell_i /= 3; + } + landmarks_ext.push_back(point); + } + write_points("landmarks/initial_landmarks",landmarks_ext); + STraits traits(&(landmarks_ext[0])); + std::vector< std::vector > WL(nbP); + + //********************** Neighbor search in a Kd tree + Tree L(boost::counting_iterator(0), + boost::counting_iterator(nb_cells*nbL), + typename Tree::Splitter(), + traits); + std::cout << "Enter (D+1) nearest landmarks\n"; + for (int i = 0; i < nbP; i++) + { + Point_d& w = W[i]; + ////Search D+1 nearest neighbours from the tree of landmarks L + K_neighbor_search search(L, w, D+1, FT(0), true, + CGAL::Distance_adapter(&(landmarks_ext[0])) ); + for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) + { + if (std::find(WL[i].begin(), WL[i].end(), (it->first)%nbL) == WL[i].end()) + WL[i].push_back((it->first)%nbL); + } + if (i == landmarks_ind[WL[i][0]]) + { + FT dist = ed.transformed_distance(W[i], landmarks[WL[i][1]]); + if (dist < lambda) + lambda = dist; + } + } + std::string out_file = "wl_result"; + write_wl(out_file,WL); + + //******************** Constructng a witness complex + std::cout << "Entered witness complex construction\n"; + Witness_complex<> witnessComplex; + witnessComplex.setNbL(nbL); + witnessComplex.witness_complex(WL); + + //******************** Making a set of bad link landmarks + std::cout << "Entered bad links\n"; + std::set< int > perturbL; + int count_badlinks = 0; + //std::cout << "Bad links around "; + std::vector< int > count_bad(D); + std::vector< int > count_good(D); + for (auto u: witnessComplex.complex_vertex_range()) + { + if (!witnessComplex.has_good_link(u, count_bad, count_good)) + { + count_badlinks++; + Point_d& l = landmarks[u]; + Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + std::vector curr_perturb; + L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); + for (int i: curr_perturb) + perturbL.insert(i%nbL); + } + } + for (unsigned int i = 0; i != count_good.size(); 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++) + if (count_bad[i] != 0) + std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; + std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; + + //*********************** Perturb bad link landmarks + for (auto u: perturbL) + { + Random_point_iterator rp(D,sqrt(lambda)/8); + std::vector point; + for (int i = 0; i < D; i++) + { + while (K().squared_distance_d_object()(*rp,origin) < lambda/256) + rp++; + FT coord = landmarks[u][i] + (*rp)[i]; + if (coord > 1) + point.push_back(coord-1); + else if (coord < -1) + point.push_back(coord+1); + else + point.push_back(coord); + } + landmarks[u] = Point_d(point); + } + std::cout << "lambda=" << lambda << std::endl; + 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); + return count_badlinks; +} + + +int main (int argc, char * const argv[]) +{ + if (argc != 5) + { + std::cerr << "Usage: " << argv[0] + << " nbP nbL dim delta\n"; + return 0; + } + int nbP = atoi(argv[1]); + int nbL = atoi(argv[2]); + int dim = atoi(argv[3]); + FT delta = atof(argv[4]); + + std::cout << "Let the carnage begin!\n"; + Point_Vector point_vector; + generate_points_random_box(point_vector, nbP, dim); + Point_Vector L; + std::vector chosen_landmarks; + bool ok=false; + while (!ok) + { + ok = true; + L = {}; + chosen_landmarks = {}; + //landmark_choice_by_delaunay(point_vector, nbP, nbL, L, chosen_landmarks, delta); + landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta); + nbL = chosen_landmarks.size(); + std::cout << "Number of landmarks is " << nbL << std::endl; + //int width = (int)pow(nbL, 1.0/dim); landmark_choice_bcc(point_vector, nbP, width, L, chosen_landmarks); + for (auto i: chosen_landmarks) + { + ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); + if (!ok) break; + } + + } + int bl = nbL, curr_min = bl; + write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_landmarks",L); + //for (int i = 0; i < 1; i++) + for (int i = 0; bl > 0; i++) + { + std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + bl=landmark_perturbation(point_vector, nbL, L, chosen_landmarks); + if (bl < curr_min) + curr_min=bl; + write_points("landmarks/landmarks0",L); + } + +} diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index d08c763f..74aae875 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -71,199 +71,6 @@ typedef CGAL::Search_traits< typename K::Construct_cartesian_const_iterator_d> Traits_base; typedef CGAL::Euclidean_distance Euclidean_distance; -/** - * \brief Class of distance in a flat torus in dimaension D - * - */ -//class Torus_distance : public Euclidean_distance { - class Torus_distance { - -public: - typedef K::FT FT; - typedef K::Point_d Point_d; - typedef Point_d Query_item; - typedef typename CGAL::Dynamic_dimension_tag D; - - double box_length = 2; - - FT transformed_distance(Query_item q, Point_d p) const - { - FT distance = FT(0); - FT coord = FT(0); - //std::cout << "Hello skitty!\n"; - typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); - typename K::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1), pit = construct_it(p); - for(; qit != qe; qit++, pit++) - { - coord = sqrt(((*qit)-(*pit))*((*qit)-(*pit))); - if (coord*coord <= (box_length-coord)*(box_length-coord)) - distance += coord*coord; - else - distance += (box_length-coord)*(box_length-coord); - } - return distance; - } - - FT min_distance_to_rectangle(const Query_item& q, - const CGAL::Kd_tree_rectangle& r) const { - FT distance = FT(0); - FT dist1, dist2; - typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); - typename K::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1); - for(unsigned int i = 0;qit != qe; i++, qit++) - { - if((*qit) < r.min_coord(i)) - { - dist1 = (r.min_coord(i)-(*qit)); - dist2 = (box_length - r.max_coord(i)+(*qit)); - if (dist1 < dist2) - distance += dist1*dist1; - else - distance += dist2*dist2; - } - else if ((*qit) > r.max_coord(i)) - { - dist1 = (box_length - (*qit)+r.min_coord(i)); - dist2 = ((*qit) - r.max_coord(i)); - if (dist1 < dist2) - distance += dist1*dist1; - else - distance += dist2*dist2; - } - } - return distance; - } - - FT min_distance_to_rectangle(const Query_item& q, - const CGAL::Kd_tree_rectangle& r, - std::vector& dists) const { - FT distance = FT(0); - FT dist1, dist2; - typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); - typename K::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1); - //std::cout << r.max_coord(0) << std::endl; - for(unsigned int i = 0;qit != qe; i++, qit++) - { - if((*qit) < r.min_coord(i)) - { - dist1 = (r.min_coord(i)-(*qit)); - dist2 = (box_length - r.max_coord(i)+(*qit)); - if (dist1 < dist2) - { - dists[i] = dist1; - distance += dist1*dist1; - } - else - { - dists[i] = dist2; - distance += dist2*dist2; - //std::cout << "Good stuff1\n"; - } - } - else if ((*qit) > r.max_coord(i)) - { - dist1 = (box_length - (*qit)+r.min_coord(i)); - dist2 = ((*qit) - r.max_coord(i)); - if (dist1 < dist2) - { - dists[i] = dist1; - distance += dist1*dist1; - //std::cout << "Good stuff2\n"; - } - else - { - dists[i] = dist2; - distance += dist2*dist2; - } - } - }; - return distance; - } - - FT max_distance_to_rectangle(const Query_item& q, - const CGAL::Kd_tree_rectangle& r) const { - FT distance=FT(0); - typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); - typename K::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1); - for(unsigned int i = 0;qit != qe; i++, qit++) - { - if (box_length <= (r.min_coord(i)+r.max_coord(i))) - if ((r.max_coord(i)+r.min_coord(i)-box_length)/FT(2.0) <= (*qit) && - (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) - distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); - else - distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); - else - if ((box_length-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || - (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) - distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); - else - distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); - } - return distance; - } - - - FT max_distance_to_rectangle(const Query_item& q, - const CGAL::Kd_tree_rectangle& r, - std::vector& dists) const { - FT distance=FT(0); - typename K::Construct_cartesian_const_iterator_d construct_it=Traits_base().construct_cartesian_const_iterator_d_object(); - typename K::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1); - for(unsigned int i = 0;qit != qe; i++, qit++) - { - if (box_length <= (r.min_coord(i)+r.max_coord(i))) - if ((r.max_coord(i)+r.min_coord(i)-box_length)/FT(2.0) <= (*qit) && - (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) - { - dists[i] = r.max_coord(i)-(*qit); - distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); - } - else - { - dists[i] = sqrt(((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i))); - distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); - } - else - if ((box_length-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || - (*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) - { - dists[i] = sqrt((r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit))); - distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit)); - - } - else - { - dists[i] = (*qit)-r.min_coord(i); - distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); - } - } - return distance; - } - - inline FT new_distance(FT dist, FT old_off, FT new_off, - int /* cutting_dimension */) const { - - FT new_dist = dist + (new_off*new_off - old_off*old_off); - return new_dist; - } - - inline FT transformed_distance(FT d) const { - return d*d; - } - - inline FT inverse_of_transformed_distance(FT d) const { - return sqrt(d); - } - -}; - - typedef std::vector< Vertex_handle > typeVectorVertex; //typedef std::pair typeSimplex; @@ -502,6 +309,49 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, } } +/** \brief A test with 600cell, the generalisation of icosaedre in 4d + */ +void landmark_choice_600cell(Point_Vector&W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind) +{ + assert(W[0].size() == 4); //4-dimensionality required + FT phi = (1+sqrt(5))/2; + FT phi_1 = FT(1)/phi; + std::vector p; + // 16 vertices + for (FT a = -0.5; a < 1; a += 1) + for (FT b = -0.5; b < 1; b += 1) + for (FT c = -0.5; c < 1; c += 1) + for (FT d = -0.5; d < 1; d += 1) + landmarks.push_back(Point_d(std::vector({a,b,c,d}))); + // 8 vertices + for (FT a = -0.5; a < 1; a += 1) + { + landmarks.push_back(Point_d(std::vector({a,0,0,0}))); + landmarks.push_back(Point_d(std::vector({0,a,0,0}))); + landmarks.push_back(Point_d(std::vector({0,0,a,0}))); + landmarks.push_back(Point_d(std::vector({0,0,0,a}))); + } + // 96 vertices + for (FT a = -phi/2; a < phi; a += phi) + for (FT b = -0.5; b < 1; b += 1) + for (FT c = -phi_1/2; c < phi_1; c += phi_1) + { + landmarks.push_back(Point_d(std::vector({a,b,c,0}))); + landmarks.push_back(Point_d(std::vector({b,a,0,c}))); + landmarks.push_back(Point_d(std::vector({c,0,a,b}))); + landmarks.push_back(Point_d(std::vector({0,c,b,a}))); + landmarks.push_back(Point_d(std::vector({a,c,0,b}))); + landmarks.push_back(Point_d(std::vector({a,0,b,c}))); + landmarks.push_back(Point_d(std::vector({c,b,0,a}))); + landmarks.push_back(Point_d(std::vector({0,b,a,c}))); + landmarks.push_back(Point_d(std::vector({b,0,c,a}))); + landmarks.push_back(Point_d(std::vector({0,a,c,b}))); + landmarks.push_back(Point_d(std::vector({b,c,a,0}))); + landmarks.push_back(Point_d(std::vector({c,a,b,0}))); + } + for (int i = 0; i < 120; ++i) + landmarks_ind.push_back(i); +} int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector& landmarks_ind) { @@ -516,10 +366,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //******************** Constructing a WL matrix int nbP = W.size(); int nbL = landmarks.size(); - //Point_Vector landmarks_ = landmarks; - Torus_distance ed; - //Equal_d ed; - //Point_d p1(std::vector({0.8,0.8})), p2(std::vector({0.1,0.1})); + Euclidean_distance ed; FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); //std::cout << "Lambda=" << lambda << std::endl; //FT lambda = 0.1;//Euclidean_distance(); @@ -578,10 +425,12 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< Witness_complex<> witnessComplex; witnessComplex.setNbL(nbL); witnessComplex.witness_complex(WL); + /* if (witnessComplex.is_witness_complex(WL)) std::cout << "!!YES. IT IS A WITNESS COMPLEX!!\n"; else - std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n"; + std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n"; + */ //******************** Making a set of bad link landmarks std::cout << "Entered bad links\n"; std::set< int > perturbL; @@ -715,11 +564,14 @@ int main (int argc, char * const argv[]) L = {}; chosen_landmarks = {}; landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); + //landmark_choice_600cell(point_vector, nbP, nbL, L, chosen_landmarks); + /* for (auto i: chosen_landmarks) { ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); if (!ok) break; } + */ } int bl = nbL, curr_min = bl; write_points("landmarks/initial_pointset",point_vector); -- cgit v1.2.3 From 1bcbd5e259d3c13211a9b71c66cacadaad35404c Mon Sep 17 00:00:00 2001 From: skachano Date: Mon, 6 Jul 2015 07:04:53 +0000 Subject: Corrected errors in witness_complex_cube git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@680 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: e0a33544c992889344cdaaed947479a5bc458ce3 --- .../example/witness_complex_cube.cpp | 59 +++++++++------------- .../example/witness_complex_flat_torus.cpp | 6 +-- .../example/witness_complex_perturbations.cpp | 20 +++++--- .../example/witness_complex_sphere.cpp | 7 +-- 4 files changed, 44 insertions(+), 48 deletions(-) (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp index 7545f156..b6051a5f 100644 --- a/src/Witness_complex/example/witness_complex_cube.cpp +++ b/src/Witness_complex/example/witness_complex_cube.cpp @@ -291,8 +291,8 @@ bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector& landmarks_ind) { - for (int j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]][l]); + for (unsigned j = 0; j < landmarks_ind.size(); ++j) + landmarks.push_back(W[landmarks_ind[j]]); } void landmark_choice_by_delaunay(Point_Vector& W, int nbP, int nbL, Point_Vector& landmarks, std::vector& landmarks_ind, FT delta) @@ -356,7 +356,7 @@ void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& } else list_it++; - fill_landmark_copies(W, landmarks, landmarks_ind); + fill_landmarks(W, landmarks, landmarks_ind); } @@ -427,6 +427,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std witnessComplex.witness_complex(WL); //******************** Making a set of bad link landmarks + /* std::cout << "Entered bad links\n"; std::set< int > perturbL; int count_badlinks = 0; @@ -453,8 +454,9 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std if (count_bad[i] != 0) std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; - + */ //*********************** Perturb bad link landmarks + /* for (auto u: perturbL) { Random_point_iterator rp(D,sqrt(lambda)/8); @@ -486,56 +488,41 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std } write_edges("landmarks/edges", witnessComplex, landmarks); return count_badlinks; + */ + return 0; } int main (int argc, char * const argv[]) { - if (argc != 5) + if (argc != 3) { std::cerr << "Usage: " << argv[0] - << " nbP nbL dim delta\n"; + << " nbP dim\n"; return 0; } int nbP = atoi(argv[1]); - int nbL = atoi(argv[2]); - int dim = atoi(argv[3]); - FT delta = atof(argv[4]); + int dim = atoi(argv[2]); std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; generate_points_random_box(point_vector, nbP, dim); Point_Vector L; std::vector chosen_landmarks; - bool ok=false; - while (!ok) - { - ok = true; - L = {}; - chosen_landmarks = {}; - //landmark_choice_by_delaunay(point_vector, nbP, nbL, L, chosen_landmarks, delta); - landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta); - nbL = chosen_landmarks.size(); - std::cout << "Number of landmarks is " << nbL << std::endl; - //int width = (int)pow(nbL, 1.0/dim); landmark_choice_bcc(point_vector, nbP, width, L, chosen_landmarks); - for (auto i: chosen_landmarks) - { - ok = ok && (std::count(chosen_landmarks.begin(),chosen_landmarks.end(),i) == 1); - if (!ok) break; - } - - } - int bl = nbL, curr_min = bl; - write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_pointset",point_vector); //write_points("landmarks/initial_landmarks",L); - //for (int i = 0; i < 1; i++) - for (int i = 0; bl > 0; i++) + for (int i = 0; i < 11; i++) + //for (int i = 0; bl > 0; i++) { - std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; - bl=landmark_perturbation(point_vector, nbL, L, chosen_landmarks); - if (bl < curr_min) - curr_min=bl; - write_points("landmarks/landmarks0",L); + //std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; + double delta = pow(10, -(1.0*i)/2); + std::cout << "delta = " << delta << std::endl; + L = {}; chosen_landmarks = {}; + landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta); + int nbL = chosen_landmarks.size(); + std::cout << "Number of landmarks = " << nbL << std::endl; + landmark_perturbation(point_vector, nbL, L, chosen_landmarks); + //write_points("landmarks/landmarks0",L); } } diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index 06bf5a9f..69ef5fbf 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -335,7 +335,7 @@ void generate_points_grid(Point_Vector& W, int width, int D) int cell_i = i; for (int l = 0; l < D; ++l) { - point.push_back(0.01*(cell_i%width)); + point.push_back((2.0/width)*(cell_i%width)); cell_i /= width; } W.push_back(point); @@ -516,7 +516,7 @@ void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, // while (!res.second) // { do chosen_landmark = rand.get_int(0,nbP); - while (std::count(landmarks_ind.begin(),landmarks_ind.end(),chosen_landmark)!=0); + 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]; @@ -672,7 +672,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< count_badlinks++; //std::cout << u << " "; Point_d& l = landmarks[u]; - Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + Fuzzy_sphere fs(l, sqrt(lambda), 0, traits); std::vector curr_perturb; L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); for (int i: curr_perturb) diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp index b3b84b1f..f78bcdab 100644 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -267,6 +267,8 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::cout << "Enter (D+1) nearest landmarks\n"; //std::cout << "Size of the tree is " << L.size() << std::endl; int D = W[0].size(); + clock_t start, end; + start = clock(); for (int i = 0; i < nbP; i++) { //std::cout << "Entered witness number " << i << std::endl; @@ -294,18 +296,23 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< lambda = dist; } //std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; - - } + } //std::cout << "\n"; + end = clock(); + std::cout << "WL matrix construction on " << nbL << " landmarks took " << (double)(end-start)/CLOCKS_PER_SEC << "s.\n"; + std::string out_file = "wl_result"; - write_wl(out_file,WL); + //write_wl(out_file,WL); //******************** Constructng a witness complex std::cout << "Entered witness complex construction\n"; Witness_complex<> witnessComplex; witnessComplex.setNbL(nbL); + start = clock(); witnessComplex.witness_complex(WL); + end = clock(); + std::cout << "Witness complex construction on " << nbL << " landmarks took " << (double)(end-start)/CLOCKS_PER_SEC << "s.\n"; //******************** Making a set of bad link landmarks std::cout << "Entered bad links\n"; std::set< int > perturbL; @@ -355,7 +362,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< std::cout << "lambda=" << lambda << std::endl; // Write the WL matrix in a file - + /* char buffer[100]; int i = sprintf(buffer,"stree_result.txt"); @@ -366,8 +373,9 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< witnessComplex.st_to_file(ofs); ofs.close(); } + */ //witnessComplex.write_badlinks("badlinks"); - write_edges_gnuplot("landmarks/edges", witnessComplex, landmarks); + //write_edges_gnuplot("landmarks/edges", witnessComplex, landmarks); return count_badlinks; } @@ -417,7 +425,7 @@ int main (int argc, char * const argv[]) file_name.erase(0, last_slash_idx + 1); } write_points("landmarks/initial_pointset",point_vector); - write_points("landmarks/initial_landmarks",L); + //write_points("landmarks/initial_landmarks",L); //for (int i = 0; bl != 0; i++) for (int i = 0; i < 1; i++) { diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index 74aae875..550c9392 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -431,6 +431,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< else std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n"; */ + */ //******************** Making a set of bad link landmarks std::cout << "Entered bad links\n"; std::set< int > perturbL; @@ -446,7 +447,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< count_badlinks++; //std::cout << u << " "; Point_d& l = landmarks[u]; - Fuzzy_sphere fs(l, sqrt(lambda)*3, 0, traits); + Fuzzy_sphere fs(l, sqrt(lambda), 0, traits); std::vector curr_perturb; L.search(std::insert_iterator>(curr_perturb,curr_perturb.begin()),fs); for (int i: curr_perturb) @@ -577,8 +578,8 @@ int main (int argc, char * const argv[]) write_points("landmarks/initial_pointset",point_vector); write_points("landmarks/initial_landmarks",L); - //for (int i = 0; bl > 0; i++) - for (int i = 0; i < 1; i++) + for (int i = 0; bl > 0; i++) + //for (int i = 0; i < 1; i++) { std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; bl=landmark_perturbation(point_vector, L, chosen_landmarks); -- cgit v1.2.3 From f44d72644fc31bc863d0c81d0ae80630ebc506e5 Mon Sep 17 00:00:00 2001 From: skachano Date: Thu, 14 Jan 2016 16:17:51 +0000 Subject: Added example stuff git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@968 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: a769b41f7d2d784a2740c149da08dac8124e3df9 --- src/Witness_complex/example/CMakeLists.txt | 11 +- src/Witness_complex/example/generators.h | 134 +++++++++++++++++++++ .../example/witness_complex_sphere.cpp | 126 +++++++++++++++++++ .../include/gudhi/Witness_complex.h | 2 +- .../test/simple_witness_complex.cpp | 5 +- .../test/witness_complex_points.cpp | 9 +- 6 files changed, 271 insertions(+), 16 deletions(-) create mode 100644 src/Witness_complex/example/generators.h create mode 100644 src/Witness_complex/example/witness_complex_sphere.cpp (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index 4473078a..b5770ba6 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -4,9 +4,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) - - add_executable( witness_complex_bench witness_complex_bench.cpp ) - add_test( witness_complex_bench_from_bunny &{CMAKE_CURRENT_BINARY_DIR}/witness_complex_bench ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) if(CGAL_FOUND) if (NOT CGAL_VERSION VERSION_LESS 4.6.0) @@ -21,12 +18,12 @@ if(CGAL_FOUND) message(STATUS "Eigen3 use file: ${EIGEN3_USE_FILE}.") include_directories (BEFORE "../../include") - add_executable ( witness_complex_bench2 witness_complex_bench2.cpp ) - target_link_libraries(witness_complex_bench2 ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + 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 Alpha shapes feature.") + 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 Alpha shapes feature. Version 4.6.0 is required.") + 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 + +#include +#include + +typedef CGAL::Epick_d K; +typedef K::FT FT; +typedef K::Point_d Point_d; +typedef std::vector Point_Vector; +typedef CGAL::Random_points_in_cube_d Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d Random_point_iterator; + +/** + * \brief Rock age method of reading off file + * + */ +inline void +off_reader_cust ( std::string file_name , std::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; + // 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 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 rp(dim,1); + for (int i = 0; i < nbP; i++) + W.push_back(*rp++); +} + + +#endif 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 . + */ +#define BOOST_PARAMETER_MAX_ARITY 12 + + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include "generators.h" + +using namespace Gudhi; + +typedef std::vector< Vertex_handle > typeVectorVertex; +//typedef std::vector< std::vector > 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 > 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 > 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/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 791d0e45..915e445c 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -27,7 +27,7 @@ #include #include #include -#include "gudhi/distance_functions.h" +#include #include #include #include diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp index c7d85a4d..ea38b5c0 100644 --- a/src/Witness_complex/test/simple_witness_complex.cpp +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -22,9 +22,8 @@ #include #include -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Simplex_tree.h" -#include "gudhi/Witness_complex.h" +#include +#include using namespace Gudhi; diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp index a86bf493..0a50101d 100644 --- a/src/Witness_complex/test/witness_complex_points.cpp +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -22,11 +22,10 @@ #include #include -//#include "gudhi/graph_simplicial_complex.h" -#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" +#include +#include +#include +#include using namespace Gudhi; -- cgit v1.2.3 From 500d18704bc26166983ea005278bd187e2d3020c Mon Sep 17 00:00:00 2001 From: skachano Date: Fri, 15 Jan 2016 08:53:03 +0000 Subject: Fixed namespaces git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@970 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f71d4840583611d3829d8c4f1049e5d2b01ce71d --- src/Witness_complex/example/witness_complex_from_file.cpp | 1 + src/Witness_complex/example/witness_complex_sphere.cpp | 1 + .../include/gudhi/Landmark_choice_by_furthest_point.h | 8 ++++++++ .../include/gudhi/Landmark_choice_by_random_point.h | 8 ++++++++ src/Witness_complex/include/gudhi/Witness_complex.h | 9 +++++---- src/Witness_complex/test/simple_witness_complex.cpp | 5 ++--- src/Witness_complex/test/witness_complex_points.cpp | 1 + 7 files changed, 26 insertions(+), 7 deletions(-) (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp index 6add4e0a..72caf30b 100644 --- a/src/Witness_complex/example/witness_complex_from_file.cpp +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -33,6 +33,7 @@ #include "gudhi/reader_utils.h" using namespace Gudhi; +using namespace Gudhi::witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; typedef std::vector< std::vector > Point_Vector; diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index 1fa0fc42..d73445b9 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -38,6 +38,7 @@ #include "generators.h" using namespace Gudhi; +using namespace Gudhi::witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; //typedef std::vector< std::vector > Point_Vector; 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 index 6ac59ae9..050286f2 100644 --- a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -23,6 +23,10 @@ #ifndef GUDHI_LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ #define GUDHI_LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ +namespace Gudhi { + +namespace witness_complex { + /** * \class Landmark_choice_by_furthest_point * \brief The class `Landmark_choice_by_furthest_point` allows to construct the matrix @@ -96,4 +100,8 @@ public: }; +} + +} + #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 index fa822591..038deff6 100644 --- a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h @@ -23,6 +23,10 @@ #ifndef GUDHI_LANDMARK_CHOICE_BY_RANDOM_POINT_H_ #define GUDHI_LANDMARK_CHOICE_BY_RANDOM_POINT_H_ +namespace Gudhi { + +namespace witness_complex { + /** * \class Landmark_choice_by_random_point * \brief The class `Landmark_choice_by_random_point` allows to construct the matrix @@ -81,4 +85,8 @@ public: }; +} + +} + #endif diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 915e445c..8938e59d 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -44,7 +44,8 @@ namespace Gudhi { - + namespace witness_complex { + /** \class Witness_complex \brief Constructs the witness complex for the given set of witnesses and landmarks. @@ -256,13 +257,13 @@ namespace Gudhi { return false; } } - return true; // Arrive here if the not_witnessed check failed all the time + return true; } -}; //class Witness_complex - + }; //class Witness_complex + } //namespace witness_complex } // namespace Guhdi diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp index ea38b5c0..86f5bcd1 100644 --- a/src/Witness_complex/test/simple_witness_complex.cpp +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -26,11 +26,10 @@ #include using namespace Gudhi; +using namespace Gudhi::witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; typedef Witness_complex> WitnessComplex; -//typedef std::pair typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; int main (int argc, char * const argv[]) { @@ -50,5 +49,5 @@ int main (int argc, char * const argv[]) 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)); + assert(witnessComplex.is_witness_complex(knn, false)); } diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp index 0a50101d..e9fac9b8 100644 --- a/src/Witness_complex/test/witness_complex_points.cpp +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -29,6 +29,7 @@ using namespace Gudhi; +using namespace Gudhi::witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; typedef Witness_complex> WitnessComplex; -- cgit v1.2.3 From d13a8867ca368b1f56be3ba151d2042728fb4754 Mon Sep 17 00:00:00 2001 From: skachano Date: Wed, 20 Jan 2016 10:43:59 +0000 Subject: Almost all Vincent's remarks git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@982 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 32564e3521e0807e26bd2c07d682751a3d47e6da --- .../example/witness_complex_sphere.cpp | 34 ++-------------------- .../test/simple_witness_complex.cpp | 4 +-- .../test/witness_complex_points.cpp | 4 +-- 3 files changed, 7 insertions(+), 35 deletions(-) (limited to 'src/Witness_complex/example/witness_complex_sphere.cpp') diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index d73445b9..9dc458d4 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -41,36 +41,9 @@ using namespace Gudhi; using namespace Gudhi::witness_complex; typedef std::vector< Vertex_handle > typeVectorVertex; -//typedef std::vector< std::vector > 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) @@ -86,15 +59,14 @@ void write_data( Data_range & data, std::string filename ) int main (int argc, char * const argv[]) { - if (argc != 3) + if (argc != 2) { std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL \n"; + << " nbL \n"; return 0; } - std::string file_name = argv[1]; - int nbL = atoi(argv[2]); + int nbL = atoi(argv[1]); clock_t start, end; // Construct the Simplex Tree diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp index 86f5bcd1..7735ca6f 100644 --- a/src/Witness_complex/test/simple_witness_complex.cpp +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Vincent Rouvreau + * Author(s): Siargey Kachanovich * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2016 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 diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp index e9fac9b8..3850bd82 100644 --- a/src/Witness_complex/test/witness_complex_points.cpp +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Vincent Rouvreau + * Author(s): Siargey Kachanovich * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2016 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 -- cgit v1.2.3