From 380c105e640335deb148feb80b81bbc6fdd39ff3 Mon Sep 17 00:00:00 2001 From: skachano Date: Fri, 29 May 2015 09:18:41 +0000 Subject: witness: Tweaking for the tests git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@597 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f357efa7f0967c337be027f37c236f3ad2864204 --- src/Witness_complex/example/CMakeLists.txt | 3 + .../example/witness_complex_flat_torus.cpp | 223 ++++++++++++++++++--- .../example/witness_complex_perturbations.cpp | 83 +++++++- .../include/gudhi/Witness_complex.h | 36 +++- 4 files changed, 294 insertions(+), 51 deletions(-) (limited to 'src/Witness_complex') diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt index e3768138..09167cbe 100644 --- a/src/Witness_complex/example/CMakeLists.txt +++ b/src/Witness_complex/example/CMakeLists.txt @@ -83,6 +83,9 @@ if(CGAL_FOUND) add_executable ( witness_complex_flat_torus witness_complex_flat_torus.cpp ) target_link_libraries(witness_complex_flat_torus ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) add_test(witness_complex_flat_torus ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_flat_torus ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + 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) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") endif() diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp index 1ce17fde..ff995a4f 100644 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp @@ -75,28 +75,33 @@ 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 Euclidean_distance { +/* + class Torus_distance { public: - typedef FT FT; - typedef Point_d Point_d; - typedef Point_d Query_item; - typedef typename CGAL::Dynamic_dimension_tag D; + 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 <= (1-coord)*(1-coord)) + if (coord*coord <= (box_length-coord)*(box_length-coord)) distance += coord*coord; else - distance += (1-coord)*(1-coord); + distance += (box_length-coord)*(box_length-coord); } return distance; } @@ -113,7 +118,7 @@ public: if((*qit) < r.min_coord(i)) { dist1 = (r.min_coord(i)-(*qit)); - dist2 = (1 - r.max_coord(i)+(*qit)); + dist2 = (box_length - r.max_coord(i)+(*qit)); if (dist1 < dist2) distance += dist1*dist1; else @@ -121,7 +126,7 @@ public: } else if ((*qit) > r.max_coord(i)) { - dist1 = (1 - (*qit)+r.min_coord(i)); + dist1 = (box_length - (*qit)+r.min_coord(i)); dist2 = ((*qit) - r.max_coord(i)); if (dist1 < dist2) distance += dist1*dist1; @@ -140,12 +145,13 @@ public: 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 = (1 - r.max_coord(i)+(*qit)); + dist2 = (box_length - r.max_coord(i)+(*qit)); if (dist1 < dist2) { dists[i] = dist1; @@ -155,16 +161,18 @@ public: { dists[i] = dist2; distance += dist2*dist2; + //std::cout << "Good stuff1\n"; } } else if ((*qit) > r.max_coord(i)) { - dist1 = (1 - (*qit)+r.min_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 { @@ -175,7 +183,7 @@ public: }; return distance; } - + FT max_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle& r) const { FT distance=FT(0); @@ -184,14 +192,14 @@ public: qe = construct_it(q,1); for(unsigned int i = 0;qit != qe; i++, qit++) { - if (FT(1) <= (r.min_coord(i)+r.max_coord(i))) - if ((r.max_coord(i)+r.min_coord(i)-FT(1))/FT(2.0) <= (*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 ((FT(1)-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + 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 @@ -200,6 +208,7 @@ public: return distance; } + FT max_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle& r, std::vector& dists) const { @@ -209,8 +218,8 @@ public: qe = construct_it(q,1); for(unsigned int i = 0;qit != qe; i++, qit++) { - if (FT(1) <= (r.min_coord(i)+r.max_coord(i))) - if ((r.max_coord(i)+r.min_coord(i)-FT(1))/FT(2.0) <= (*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); @@ -222,7 +231,7 @@ public: distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i)); } else - if ((FT(1)-r.max_coord(i)-r.min_coord(i))/FT(2.0) <= (*qit) || + 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))); @@ -237,7 +246,14 @@ public: } 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; } @@ -247,7 +263,7 @@ public: } }; - +*/ typedef std::vector< Vertex_handle > typeVectorVertex; @@ -259,7 +275,7 @@ typedef CGAL::Search_traits_adapter< //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 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; @@ -275,7 +291,7 @@ 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; - +bool toric=true; /** * \brief Customized version of read_points @@ -319,6 +335,14 @@ void generate_points_random_box(Point_Vector& W, int nbP, int dim) } } +/* 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) { @@ -359,6 +383,97 @@ void write_wl( std::string file_name, std::vector< std::vector > & WL) } +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 @@ -397,21 +512,41 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::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(); + std::vector 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)); + 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); Tree L(boost::counting_iterator(0), - boost::counting_iterator(nbL), + boost::counting_iterator(nb_cells*nbL), typename Tree::Splitter(), - STraits(&(landmarks[0]))); + traits); /*Tree2 L2(boost::counting_iterator(0), boost::counting_iterator(nbL), typename Tree::Splitter(), @@ -424,10 +559,14 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //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 + ////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])) ); + CGAL::Distance_adapter(&(landmarks_ext[0])) ); //std::cout << "Safely found nearest landmarks\n"; for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) { @@ -435,7 +574,7 @@ 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); + WL[i].push_back((it->first)%nbL); //std::cout << "ITFIRST " << it->first << std::endl; //std::cout << i << " " << it->first << ": " << it->second << std::endl; } @@ -463,6 +602,8 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< int count_badlinks = 0; //std::cout << "Bad links around "; for (auto u: witnessComplex.complex_vertex_range()) + { + std::cout << "Vertex " << u << " "; if (!witnessComplex.has_good_link(u)) { //std::cout << "Landmark " << u << " start!" << std::endl; @@ -470,12 +611,16 @@ 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, STraits(&(landmarks[0]))); - L.search(std::insert_iterator>(perturbL,perturbL.begin()),fs); + 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 @@ -523,7 +668,18 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< witnessComplex.st_to_file(ofs); ofs.close(); } - + /* + i = sprintf(buffer,"badlinks.txt"); + if (i >= 0) + { + std::string out_file = (std::string)buffer; + std::ofstream ofs (out_file, std::ofstream::out); + witnessComplex.write_bad_links(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; } @@ -585,13 +741,16 @@ int main (int argc, char * const argv[]) } } int bl = nbL, curr_min = bl; + write_points("landmarks/initial_pointset",point_vector); + //write_points("landmarks/initial_landmarks",L); - 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); if (bl < curr_min) curr_min=bl; + write_points("landmarks/landmarks0",L); } //end = clock(); diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp index 514b115b..07b1eb8d 100644 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ b/src/Witness_complex/example/witness_complex_perturbations.cpp @@ -157,11 +157,44 @@ void write_wl( std::string file_name, std::vector< std::vector > & WL) ofs.close(); } +void write_points( std::string file_name, std::vector< Point_d > & WL) +{ + std::ofstream ofs (file_name, std::ofstream::out); + for (auto w : WL) + { + for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) + ofs << *it << " "; + ofs << "\n"; + } + ofs.close(); +} + +void write_edges_gnuplot(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"; @@ -183,7 +216,32 @@ 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; + //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) { @@ -285,13 +343,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector< //std::cout << "landmark[0][0] after" << landmarks[0][0] << std::endl; std::cout << "lambda=" << lambda << std::endl; // Write the WL matrix in a file - /* - mkdir("output", S_IRWXU); - const size_t last_slash_idx = file_name.find_last_of("/"); - if (std::string::npos != last_slash_idx) - { - file_name.erase(0, last_slash_idx + 1); - } + char buffer[100]; int i = sprintf(buffer,"stree_result.txt"); @@ -303,7 +355,8 @@ 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); return count_badlinks; } @@ -345,10 +398,22 @@ int main (int argc, char * const argv[]) //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks); int bl = 1; + + mkdir("landmarks", S_IRWXU); + const size_t last_slash_idx = file_name.find_last_of("/"); + if (std::string::npos != last_slash_idx) + { + file_name.erase(0, last_slash_idx + 1); + } + 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 << " ========\n"; bl = landmark_perturbation(point_vector, L, chosen_landmarks); + std::ostringstream os(std::ostringstream::ate);; + os << "landmarks/landmarks0"; + write_points(os.str(),L); } //end = clock(); diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index bb553347..c8eaa3ef 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -642,14 +642,16 @@ private: 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); //std::cout << " dim " << d << "\n"; //Siblings* curr_sibl = root(); + //std::cout << "Currently at vertex " return (link_is_pseudomanifold(link_vertices,d)); } - + /** \brief Search and output links around vertices that are not pseudomanifolds * */ @@ -661,7 +663,7 @@ private: //int count = 0; for (auto v: complex_vertex_range()) { - //std::cout << "Vertex " << v << ":\n"; + std::cout << "Vertex " << v << ": "; std::vector< Vertex_handle > link_vertices; // Fill link_vertices for (auto u: complex_vertex_range()) @@ -670,10 +672,10 @@ private: if (u != v && find(edge) != null_simplex()) link_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); @@ -718,6 +720,7 @@ private: 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) { curr_simplex.push_back(*it); @@ -731,8 +734,13 @@ private: { //std::cout << " -> " << *it << "\n"; int d = link_dim(link_vertices, it+1, curr_d+1, curr_simplex); - if (d > final_d) - final_d = d; + if (d >= final_d) + { + final_d = d; + std::cout << d << " "; + print_vector(curr_simplex); + std::cout << std::endl; + } } /* else @@ -781,14 +789,22 @@ 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; } } // 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); + //std::vector components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); } void add_vertices(typeVectorVertex& link_vertices, -- cgit v1.2.3