diff options
Diffstat (limited to 'src/Witness_complex/example/witness_complex_cube.cpp')
| -rw-r--r-- | src/Witness_complex/example/witness_complex_cube.cpp | 484 |
1 files changed, 435 insertions, 49 deletions
diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp index b6051a5f..a9a2959b 100644 --- a/src/Witness_complex/example/witness_complex_cube.cpp +++ b/src/Witness_complex/example/witness_complex_cube.cpp @@ -47,12 +47,15 @@ #include <CGAL/Kd_tree.h> #include <CGAL/Euclidean_distance.h> #include <CGAL/Kernel_d/Sphere_d.h> +#include <CGAL/Kernel_d/Hyperplane_d.h> +#include <CGAL/enum.h> #include <CGAL/Kernel_d/Vector_d.h> #include <CGAL/point_generators_d.h> #include <CGAL/constructions_d.h> #include <CGAL/Fuzzy_sphere.h> #include <CGAL/Random.h> +#include <CGAL/Timer.h> #include <CGAL/Delaunay_triangulation.h> @@ -66,7 +69,10 @@ using namespace Gudhi; typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; typedef K::Point_d Point_d; -//typedef CGAL::Cartesian_d<double> K; +typedef K::Vector_d Vector_d; +typedef K::Oriented_side_d Oriented_side_d; +typedef K::Has_on_positive_side_d Has_on_positive_side_d; + //typedef CGAL::Point_d<K> Point_d; typedef K::FT FT; typedef CGAL::Search_traits< @@ -105,7 +111,12 @@ typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; typedef Delaunay_triangulation::Facet Facet; -typedef CGAL::Sphere_d<K> Sphere_d; +typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; +typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; +//typedef CGAL::Sphere_d<K> Sphere_d; +typedef K::Sphere_d Sphere_d; +typedef K::Hyperplane_d Hyperplane_d; + bool toric=false; @@ -230,6 +241,170 @@ void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, landmark_count++; } +bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) +{ + for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) + if (*v_it == v) + return true; + return false; +} + +bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, bool is_infinite, const Point_d& p, FT delta) +{ + if (!is_infinite) + // FINITE CASE + { + Sphere_d cs(vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + { + //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); + //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) + if (dist2 >= r*r && dist2 <= r*r+delta*delta) + return true; + } + } + } + else + // INFINITE CASE + { + Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); + while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) + v++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); + Vector_d orth_v = facet_plane.orthogonal_vector(); + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (!t.is_infinite(v_it)) + if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) + { + std::vector<FT> coords; + Point_d p = v_it->point(); + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + if (!p_is_inside && p_delta_is_inside) + return true; + } + } + return false; +} + + +bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells) +{ + Euclidean_distance ed; + std::vector<Point_d> vertices; + if (!t.is_infinite(c)) + { + // if the cell is finite, we look if the protection is violated + for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) + vertices.push_back((*v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d center_cs = cs.center(); + FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); + FT dist2 = ed.transformed_distance(center_cs, p); + // if the new point is inside the protection ball of a non conflicting simplex + //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) + if (dist2 >= r*r && dist2 <= r*r+delta*delta) + return true; + c->tds_data().mark_visited(); + marked_cells.push_back(c); + // if the new point is inside the circumscribing ball : continue violation searching on neughbours + if (dist2 < r*r) + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells)) + return true; + } + // if the new point is outside the protection sphere + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is guaranteed to be finite + vertices.clear(); vertices.push_back(p); + for (int i = 0; i < D+1; ++i) + if (i != index) + vertices.push_back(parent_cell->vertex(i)->point()); + Delaunay_vertex vertex_to_check; + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, false, vertex_to_check->point(), delta)) + return true; + } + } + else + { + // Inside of the convex hull is + side. Outside is - side. + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!t.is_infinite(*vh_it)) + vertices.push_back((*vh_it)->point()); + Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); + while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) + v_it++; + Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); + //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); + Vector_d orth_v = facet_plane.orthogonal_vector(); + std::vector<FT> coords; + auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); + for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) + coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); + Point_d p_delta = Point_d(coords); + bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); + bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); + + if (!p_is_inside && p_delta_is_inside) + return true; + //if the cell is infinite we look at the neighbours regardless + c->tds_data().mark_visited(); + marked_cells.push_back(c); + if (p_is_inside) + for (int i = 0; i < D+1; ++i) + { + Full_cell_handle next_c = c->neighbor(i); + if (next_c->tds_data().is_clear() && + is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells)) + return true; + } + else + { + // facet f is on the border of the conflict zone : check protection of simplex {p,f} + // the new simplex is finite if the parent cell is finite + vertices.clear(); vertices.push_back(p); + bool new_simplex_is_finite = false; + for (int i = 0; i < D+1; ++i) + if (i != index) + { + if (t.is_infinite(parent_cell->vertex(i))) + new_simplex_is_finite = true; + else + vertices.push_back(parent_cell->vertex(i)->point()); + } + Delaunay_vertex vertex_to_check; + for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) + if (!vertex_is_in_full_cell(*vh_it, parent_cell)) + { + vertex_to_check = *vh_it; break; + } + if (new_cell_is_violated(t, vertices, new_simplex_is_finite, vertex_to_check->point(), delta)) + return true; + } + } + return false; +} + bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) { Euclidean_distance ed; @@ -238,6 +413,22 @@ bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT de Delaunay_triangulation::Facet ft; Delaunay_triangulation::Full_cell_handle c; Delaunay_triangulation::Locate_type lt; + std::vector<Full_cell_handle> marked_cells; + c = t.locate(p, lt, f, ft, v); + bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells); + for (Full_cell_handle fc : marked_cells) + fc->tds_data().clear(); + return violation_existing_cells; +} + +bool old_is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) +{ + 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)) @@ -245,7 +436,7 @@ bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT de std::vector<Point_d> vertices; for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) vertices.push_back((*v_it)->point()); - Sphere_d cs(D, vertices.begin(), vertices.end()); + Sphere_d cs( vertices.begin(), vertices.end()); Point_d center_cs = cs.center(); FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); FT dist2 = ed.transformed_distance(center_cs, p); @@ -253,39 +444,88 @@ bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT de if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) return true; } + t.insert(p, c); return false; } +void write_delaunay_mesh(Delaunay_triangulation& t, const Point_d& p) +{ + std::ofstream ofs ("delaunay.mesh", std::ofstream::out); + int nbV = t.number_of_vertices()+1; + ofs << "MeshVersionFormatted 1\nDimension 2\n"; + ofs << "Vertices\n" << nbV << "\n"; + int ind = 1; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + { + if (t.is_infinite(v_it)) + continue; + for (auto coord = v_it->point().cartesian_begin(); coord != v_it->point().cartesian_end(); ++coord) + ofs << *coord << " "; + ofs << "508\n"; + index_of_vertex[v_it] = ind++; + } + for (auto coord = p.cartesian_begin(); coord != p.cartesian_end(); ++coord) + ofs << *coord << " "; + ofs << "208\n"; + /* + int nbFacets = 0; + for (auto ft_it = t.finite_facets_begin(); ft_it != t.finite_facets_end(); ++ft_it) + nbFacets++; + ofs << "\nEdges\n" << nbFacets << "\n\n"; + for (auto ft_it = t.facets_begin(); ft_it != t.facets_end(); ++ft_it) + { + if (t.is_infinite(ft_it)) + continue; + for (auto vh_it = ft_it->vertices_begin(); vh_it != ft_it->vertices_end(); ++vh_it) + ofs << index_of_vertex[*vh_it] << " "; + } + */ + ofs << "Triangles " << t.number_of_finite_full_cells()+1 << "\n"; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) + ofs << index_of_vertex[*vh_it] << " "; + ofs << "508\n"; + } + ofs << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; + ofs << "End\n"; + ofs.close(); +} + bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) { + // Verification part Euclidean_distance ed; - 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) - { + if (!t.is_infinite(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) + if (!vertex_is_in_full_cell(v_it, fc_it)) { std::vector<Point_d> vertices; for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) vertices.push_back((*fc_v_it)->point()); - Sphere_d cs(D, vertices.begin(), vertices.end()); + Sphere_d cs( vertices.begin(), vertices.end()); Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point())); + FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(0)->point())); FT dist2 = ed.transformed_distance(center_cs, v_it->point()); //if the new point is inside the protection ball of a non conflicting simplex - if (dist2 <= (r+delta)*(r+delta)) - return false; + //std::cout << "Dist^2 = " << dist2 << " (r+delta)*(r+delta) = " << (r+delta)*(r+delta) << " r^2 = " << r*r <<"\n"; + //if (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r) + if (dist2 <= r*r+delta*delta && dist2 >= r*r) + { + write_delaunay_mesh(t, v_it->point()); + std::cout << "Problematic vertex " << *v_it << " "; + std::cout << "Problematic cell " << *fc_it << "\n"; + std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", r^2+delta^2 = " << r*r+delta*delta << "\n"; + return false; + } } - } + return true; } @@ -295,33 +535,65 @@ void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks.push_back(W[landmarks_ind[j]]); } -void landmark_choice_by_delaunay(Point_Vector& W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT delta) +void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) { - 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) + // Store vertex indices in a map + int ind = 0; //index of a vertex + std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; + for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) + if (t.is_infinite(v_it)) + continue; + else + index_of_vertex[v_it] = ind++; + // Write full cells as vectors in full_cells + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) { - 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); + if (t.is_infinite(fc_it)) + continue; + std::vector<int> cell; + for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it) + cell.push_back(index_of_vertex[*v_it]); + full_cells.push_back(cell); } - 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); +} + +FT sampling_radius(Delaunay_triangulation& t) +{ + FT epsilon2 = 4.0; + for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) + { + if (t.is_infinite(fc_it)) + continue; + Point_Vector vertices; + for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it) + vertices.push_back((*fc_v_it)->point()); + Sphere_d cs( vertices.begin(), vertices.end()); + Point_d csc = cs.center(); + bool in_cube = true; + for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi) + if (*xi > 1.0 || *xi < -1.0) + { + in_cube = false; break; + } + if (!in_cube) + continue; + FT 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<int>& landmarks_ind, FT delta) +void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT delta, std::vector<std::vector<int>>& full_cells) { - int D = W[0].size(); + unsigned D = W[0].size(); Torus_distance td; Euclidean_distance ed; Delaunay_triangulation t(D); @@ -335,32 +607,106 @@ void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& temp_vector.push_back(i); unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed)); + //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end()); for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it) index_list.push_front(*it); } - // add the first D+1 vertices to form one non-empty cell + for (unsigned pos1 = 0; pos1 < D+1; ++pos1) + { + std::vector<FT> point; + for (unsigned i = 0; i < pos1; ++i) + point.push_back(-1); + if (pos1 != D) + point.push_back(1); + for (unsigned i = pos1+1; i < D; ++i) + point.push_back(0); + assert(point.size() == D); + W[index_list.front()] = Point_d(point); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); + index_list.pop_front(); + } + // add the first D+1 vertices to form one finite cell + /* for (int i = 0; i <= D+1; ++i) { + t.insert(W[index_list.front()]); insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); index_list.pop_front(); } + */ + /* + { + std::vector<FT> coords; + for (int i = 0; i < D; ++i) + coords.push_back(-1); + W[index_list.front()] = Point_d(coords); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); + index_list.pop_front(); + for (int i = 0; i < D; ++i) + { + coords.clear(); + for (int j = 0; j < D; ++j) + if (i == j) + coords.push_back(1); + else + coords.push_back(-1); + W[index_list.front()] = Point_d(coords); + insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); + index_list.pop_front(); + } + } + */ + //std::cout << t; + //assert(t.number_of_vertices() == D+1); + //assert(landmarks_ind.size() == D+1); + //assert(W[landmarks_ind[0]][0] == 0); // add other vertices if they don't violate protection std::list<int>::iterator list_it = index_list.begin(); while (list_it != index_list.end()) - if (!is_violating_protection(W[*list_it], t, D, delta)) - { + { + if (!is_violating_protection(W[*list_it], t, D, delta)) + { // If no conflicts then insert in every copy of T^3 + is_violating_protection(W[*list_it], t, D, delta); insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); index_list.erase(list_it); list_it = index_list.begin(); + //std::cout << "index_list_size() = " << index_list.size() << "\n"; } - else - list_it++; + else + { + list_it++; + //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; + } + //write_delaunay_mesh(t, W[*list_it]); + } fill_landmarks(W, landmarks, landmarks_ind); + fill_full_cell_vector(t, full_cells); + if (triangulation_is_protected(t, delta)) + std::cout << "Triangulation is ok\n"; + else + std::cout << "Triangulation is BAD!! T_T しくしく!\n"; + write_delaunay_mesh(t, Point_d(std::vector<FT>({0,0}))); + //std::cout << t << std::endl; } +template <typename T> +void print_vector(std::vector<T> v) +{ + std::cout << "["; + if (!v.empty()) + { + std::cout << *(v.begin()); + for (auto it = v.begin()+1; it != v.end(); ++it) + { + std::cout << ","; + std::cout << *it; + } + } + std::cout << "]"; +} -int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind) +int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind, std::vector<std::vector<int>>& full_cells) { //******************** Preface: origin point int D = W[0].size(); @@ -426,6 +772,20 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std witnessComplex.setNbL(nbL); witnessComplex.witness_complex(WL); + //******************** Verifying if all full cells are in the complex + + int in=0, not_in=0; + for (auto cell : full_cells) + { + //print_vector(cell); + if (witnessComplex.find(cell) != witnessComplex.null_simplex()) + in++; + else + not_in++; + } + std::cout << "Out of all the cells in Delaunay triangulation:\n" << in << " are in the witness complex\n" << + not_in << " are not.\n"; + //******************** Making a set of bad link landmarks /* std::cout << "Entered bad links\n"; @@ -476,6 +836,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std landmarks[u] = Point_d(point); } std::cout << "lambda=" << lambda << std::endl; + */ char buffer[100]; int i = sprintf(buffer,"stree_result.txt"); @@ -486,7 +847,9 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std witnessComplex.st_to_file(ofs); ofs.close(); } + write_edges("landmarks/edges", witnessComplex, landmarks); + /* return count_badlinks; */ return 0; @@ -495,22 +858,27 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std int main (int argc, char * const argv[]) { - if (argc != 3) + if (argc != 4) { std::cerr << "Usage: " << argv[0] - << " nbP dim\n"; + << " nbP dim delta\n"; return 0; } int nbP = atoi(argv[1]); int dim = atoi(argv[2]); + double delta = atof(argv[3]); 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<int> chosen_landmarks; //write_points("landmarks/initial_pointset",point_vector); //write_points("landmarks/initial_landmarks",L); + CGAL::Timer timer; + /* for (int i = 0; i < 11; i++) //for (int i = 0; bl > 0; i++) { @@ -518,11 +886,29 @@ int main (int argc, char * const argv[]) 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); + std::vector<std::vector<int>> full_cells; + timer.start(); + landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, full_cells); + timer.stop(); + FT epsilon2 = point_sampling_radius_by_delaunay(L); + std::cout << "Final epsilon = " << epsilon2 << ". Ratio = " << epsilon/epsilon2 << std::endl; + write_points("landmarks/initial_landmarks",L); int nbL = chosen_landmarks.size(); - std::cout << "Number of landmarks = " << nbL << std::endl; - landmark_perturbation(point_vector, nbL, L, chosen_landmarks); + std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl; + landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); + timer.reset(); //write_points("landmarks/landmarks0",L); } - + */ + std::vector<std::vector<int>> full_cells; + timer.start(); + landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, full_cells); + timer.stop(); + FT epsilon2 = point_sampling_radius_by_delaunay(L); + std::cout << "Final epsilon = " << epsilon2 << ". Ratio = " << epsilon/epsilon2 << std::endl; + write_points("landmarks/initial_landmarks",L); + int nbL = chosen_landmarks.size(); + std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl; + landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); + timer.reset(); } |