diff options
Diffstat (limited to 'src/Witness_complex')
16 files changed, 0 insertions, 7892 deletions
diff --git a/src/Witness_complex/example/Torus_distance.h b/src/Witness_complex/example/Torus_distance.h deleted file mode 100644 index 5ae127df..00000000 --- a/src/Witness_complex/example/Torus_distance.h +++ /dev/null @@ -1,209 +0,0 @@ -#ifndef GUDHI_TORUS_DISTANCE_H_ -#define GUDHI_TORUS_DISTANCE_H_ - -#include <math.h> - -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/Epick_d.h> - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> 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<FT,D>& 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<FT,D>& r, - std::vector<FT>& 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<FT,D>& 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<FT,D>& r, - std::vector<FT>& 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/protected_sets/output_tikz.h b/src/Witness_complex/example/protected_sets/output_tikz.h deleted file mode 100644 index edfd9a5f..00000000 --- a/src/Witness_complex/example/protected_sets/output_tikz.h +++ /dev/null @@ -1,67 +0,0 @@ -#ifndef OUTPUT_TIKZ_H -#define OUTPUT_TIKZ_H - -#include <vector> -#include <string> -#include <algorithm> -#include <fstream> -#include <cmath> - -void write_tikz_plot(std::vector<FT> data, std::string filename) -{ - int n = data.size(); - FT vmax = *(std::max_element(data.begin(), data.end())); - //std::cout << std::log10(vmax) << " " << std::floor(std::log10(vmax)); - - FT order10 = pow(10,std::floor(std::log10(vmax))); - int digit = std::floor( vmax / order10) + 1; - if (digit == 4 || digit == 6) digit = 5; - if (digit > 6) digit = 10; - FT plot_max = digit*order10; - std::cout << plot_max << " " << vmax; - FT hstep = 10.0/(n-1); - FT wstep = 10.0 / plot_max; - - std::cout << "(eps_max-eps_min)/(N-48) = " << (vmax-*data.begin())/(data.size()-48) << "\n"; - std::ofstream ofs(filename, std::ofstream::out); - - ofs << - "\\documentclass{standalone}\n" << - "\\usepackage[utf8]{inputenc}\n" << - "\\usepackage{amsmath}\n" << - "\\usepackage{tikz}\n\n" << - "\\begin{document}\n" << - "\\begin{tikzpicture}\n"; - - ofs << "\\draw[->] (0,0) -- (0,11);" << std::endl << - "\\draw[->] (0,0) -- (11,0);" << std::endl << - "\\foreach \\i in {1,...,10}" << std::endl << - "\\draw (0,\\i) -- (-0.05,\\i);" << std::endl << - "\\foreach \\i in {1,...,10}" << std::endl << - "\\draw (\\i,0) -- (\\i,-0.05);" << std::endl << std::endl << - - "\\foreach \\i in {1,...,10}" << std::endl << - "\\draw[dashed] (-0.05,\\i) -- (11,\\i);" << std::endl << std::endl << - - "\\node at (-0.5,11) {$*$}; " << std::endl << - "\\node at (11,-0.5) {$*$}; " << std::endl << - "\\node at (-0.5,-0.5) {0}; " << std::endl << - "\\node at (-0.5,10) {" << plot_max << "}; " << std::endl << - "%\\node at (10,-0.5) {2}; " << std::endl; - - ofs << "\\draw[red] (0," << wstep*data[0] << ")"; - for (int i = 1; i < n; ++i) - ofs << " -- (" << hstep*i << "," << wstep*data[i] << ")"; - ofs << ";\n"; - - ofs << - "\\end{tikzpicture}\n" << - "\\end{document}"; - - ofs.close(); - - - -} - -#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets.h b/src/Witness_complex/example/protected_sets/protected_sets.h deleted file mode 100644 index ec627808..00000000 --- a/src/Witness_complex/example/protected_sets/protected_sets.h +++ /dev/null @@ -1,597 +0,0 @@ -#ifndef PROTECTED_SETS_H -#define PROTECTED_SETS_H - -#include <algorithm> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.h> -#include <CGAL/Kernel_d/Hyperplane_d.h> -#include <CGAL/Kernel_d/Vector_d.h> - -#include <CGAL/point_generators_d.h> -#include <CGAL/constructions_d.h> - - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -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 K::Sphere_d Sphere_d; -typedef K::Hyperplane_d Hyperplane_d; - -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; -typedef Delaunay_triangulation::Facet Facet; -typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; -typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; - -typedef std::vector<Point_d> Point_Vector; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - -FT _sfty = pow(10,-14); - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// AUXILLARY FUNCTIONS -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well - * - * W is the initial point vector - * chosen_landmark is the index of the chosen point in W - * landmarks_ind is the vector of indices of already chosen points in W - * delaunay is the Delaunay triangulation - * landmark_count is the current number of chosen vertices - * torus is true iff you are working on a flat torus [-1,1]^d - * OUT: Vertex handle to the newly inserted point - */ -Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) -{ - if (!torus) - { - Delaunay_vertex v =delaunay.insert(W[chosen_landmark]); - landmarks_ind.push_back(chosen_landmark); - landmark_count++; - return v; - } - else - { - int D = W[0].size(); - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> 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; - } - v = delaunay.insert(point); - } - landmarks_ind.push_back(chosen_landmark); - landmark_count++; - return v; - } -} - -/** Small check if the vertex v is in the full cell fc - */ - -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; -} - -/** Fill chosen point vector from indices with copies if you are working on a flat torus - * - * IN: W is the point vector - * OUT: landmarks is the output vector - * IN: landmarks_ind is the vector of indices - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) -{ - if (!torus) - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); - else - { - 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); - } - } -} - -/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices - * - * IN: t is the Delaunay triangulation - * OUT: full_cells is the output vector - */ - -void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) -{ - // 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) - { - 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; - 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); - } -} - -//////////////////////////////////////////////////////////////////////////////////////////////////////////// -// IS VIOLATED TEST -//////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Check if a newly created cell is protected from old vertices - * - * t is the Delaunay triangulation - * vertices is the vector containing the point to insert and a facet f in t - * v1 is the vertex of t, such that f and v1 form a simplex - * v2 is the vertex of t, such that f and v2 form another simplex - * delta is the protection constant - * power_protection is true iff the delta-power protection is used - */ - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) -{ - assert(vertices.size() == vertices[0].size() || - vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 - assert(v1 != v2); - if (vertices.size() == vertices[0].size() + 1) - // 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - } - */ - // Check if the simplex is thick enough - Hyperplane_d tau_h(vertices.begin()+1, vertices.end()); - Vector_d orth_tau = tau_h.orthogonal_vector(); - /* - p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1)); - */ - //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n"; - FT orth_length = sqrt(orth_tau.squared_length()); - K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; - // Compute the altitude - FT h = 0; - for (o_it = orth_tau.cartesian_begin(), - p_it = vertices.begin()->cartesian_begin(), - s_it = (vertices.begin()+1)->cartesian_begin(); - o_it != orth_tau.cartesian_end(); - ++o_it, ++p_it, ++s_it) - h += (*o_it)*(*p_it - *s_it)/orth_length; - h = fabs(h); - // Is the center inside the box? - bool inside_the_box = true; - for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) - if (*c_it > 1.0 || *c_it < -1.0) - { - inside_the_box = false; break; - } - if (inside_the_box && h/r < theta0) - return true; - if (!t.is_infinite(v1)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - if (!t.is_infinite(v2)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && 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; - } - */ - if (!t.is_infinite(v1)) - { - std::vector<FT> coords; - Point_d p = v1->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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - if (!t.is_infinite(v2)) - { - std::vector<FT> coords; - Point_d p = v2->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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - -/** Auxillary recursive function to check if the point p violates the protection of the cell c and - * if there is a violation of an eventual new cell - * - * p is the point to insert - * t is the current triangulation - * c is the current cell (simplex) - * parent_cell is the parent cell (simplex) - * index is the index of the facet between c and parent_cell from parent_cell's point of view - * D is the dimension of the triangulation - * delta is the protection constant - * marked_cells is the vector of all visited cells containing p in their circumscribed ball - * power_protection is true iff you are working with delta-power protection - * - * OUT: true iff inserting p hasn't produced any violation so far - */ - -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, bool power_protection, FT theta0) -{ - 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - // if the new point is inside the circumscribing ball : continue violation searching on neighbours - //if (dist2 < r*r) - //if (dist2 < (5*r+delta)*(5*r+delta)) - if (dist2 < r*r) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, power_protection, theta0)) - 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 = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, 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) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - // If we work with power protection, we just ignore any conflicts - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - if (p_is_inside) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, power_protection, theta0)) - 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); - for (int i = 0; i < D+1; ++i) - if (i != index) - if (!t.is_infinite(parent_cell->vertex(i))) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) - return true; - } - } - //c->tds_data().clear_visited(); - //marked_cells.pop_back(); - return false; -} - -/** Checks if inserting the point p in t will make conflicts - * - * p is the point to insert - * t is the current triangulation - * D is the dimension of triangulation - * delta is the protection constant - * power_protection is true iff you are working with delta-power protection - * OUT: true iff inserting p produces a violation of delta-protection. - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) -{ - 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; - 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, power_protection, theta0); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// -//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////////////////////// -// LANDMARK CHOICE PROCEDURE -/////////////////////////////////////////////////////////////////////// - -/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. - * - * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d - * IN: nbP is the size of W - * OUT: landmarks is the output vector for the points - * OUT: landmarks_ind is the output vector for the indices of the selected points in W - * IN: delta is the constant of protection - * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -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, bool torus, bool power_protection, FT theta0) -{ - bool return_ = true; - unsigned D = W[0].size(); - Torus_distance td; - Euclidean_distance ed; - Delaunay_triangulation t(D); - CGAL::Random rand; - int landmark_count = 0; - std::list<int> index_list; - // shuffle the list of indexes (via a vector) - { - std::vector<int> 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)); - //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); - } - if (!torus) - 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, torus); - index_list.pop_front(); - } - else if (D == 2) - { - for (int i = 0; i < 4; ++i) - for (int j = 0; j < 2; ++j) - { - W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0}); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); - index_list.pop_front(); - W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j}); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); - index_list.pop_front(); - } - } - else - std::cout << "No torus starter available for dim>2\n"; - 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, power_protection, theta0)) - { - // If no conflicts then insert in every copy of T^3 - - insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus); - if (return_) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - index_list.erase(list_it++); - /* - // PIECE OF CODE FOR DEBUGGING PURPOSES - - Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - if (triangulation_is_protected(t, delta)) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - { //THAT'S WHERE SOMETHING'S WRONG - t.remove(inserted_v); - landmarks_ind.pop_back(); - landmark_count--; - write_delaunay_mesh(t, W[*list_it], is2d); - is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore - } - */ - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - list_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //if (list_it != index_list.end()) - // write_delaunay_mesh(t, W[*list_it], is2d); - } - fill_landmarks(W, landmarks, landmarks_ind, torus); - 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, W[0], is2d); - //std::cout << t << std::endl; -} - -#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp b/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp deleted file mode 100644 index f3df3f1e..00000000 --- a/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp +++ /dev/null @@ -1,610 +0,0 @@ -#ifndef PROTECTED_SETS_H -#define PROTECTED_SETS_H - -#include <algorithm> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.h> -#include <CGAL/Kernel_d/Hyperplane_d.h> -#include <CGAL/Kernel_d/Vector_d.h> - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -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 K::Sphere_d Sphere_d; -typedef K::Hyperplane_d Hyperplane_d; - -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; -typedef Delaunay_triangulation::Facet Facet; -typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; -typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; - -typedef std::vector<Point_d> Point_Vector; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - -FT _sfty = pow(10,-14); - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// AUXILLARY FUNCTIONS -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well - * - * W is the initial point vector - * chosen_landmark is the index of the chosen point in W - * landmarks_ind is the vector of indices of already chosen points in W - * delaunay is the Delaunay triangulation - * landmark_count is the current number of chosen vertices - * torus is true iff you are working on a flat torus [-1,1]^d - * OUT: Vertex handle to the newly inserted point - */ -Delaunay_vertex insert_delaunay_landmark_with_copies(Point_d& p, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) -{ - if (!torus) - { - Delaunay_vertex v =delaunay.insert(p); - landmark_count++; - return v; - } - else - { - int D = W[0].size(); - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> point; - int cell_i = i; - for (int l = 0; l < D; ++l) - { - point.push_back(p[l] + 2.0*(cell_i%3-1)); - cell_i /= 3; - } - v = delaunay.insert(point); - } - landmark_count++; - return v; - } -} - -/** Small check if the vertex v is in the full cell fc - */ - -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; -} - -/** Fill chosen point vector from indices with copies if you are working on a flat torus - * - * IN: W is the point vector - * OUT: landmarks is the output vector - * IN: landmarks_ind is the vector of indices - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) -{ - if (!torus) - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); - else - { - 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); - } - } -} - -/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices - * - * IN: t is the Delaunay triangulation - * OUT: full_cells is the output vector - */ - -void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) -{ - // 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) - { - 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; - 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); - } -} - -//////////////////////////////////////////////////////////////////////////////////////////////////////////// -// IS VIOLATED TEST -//////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Check if a newly created cell is protected from old vertices - * - * t is the Delaunay triangulation - * vertices is the vector containing the point to insert and a facet f in t - * v1 is the vertex of t, such that f and v1 form a simplex - * v2 is the vertex of t, such that f and v2 form another simplex - * delta is the protection constant - * power_protection is true iff the delta-power protection is used - */ - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) -{ - assert(vertices.size() == vertices[0].size() || - vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 - assert(v1 != v2); - if (vertices.size() == vertices[0].size() + 1) - // 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - } - */ - // Check if the simplex is thick enough - Hyperplane_d tau_h(vertices.begin()+1, vertices.end()); - Vector_d orth_tau = tau_h.orthogonal_vector(); - /* - p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1)); - */ - //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n"; - FT orth_length = sqrt(orth_tau.squared_length()); - K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; - // Compute the altitude - FT h = 0; - for (o_it = orth_tau.cartesian_begin(), - p_it = vertices.begin()->cartesian_begin(), - s_it = (vertices.begin()+1)->cartesian_begin(); - o_it != orth_tau.cartesian_end(); - ++o_it, ++p_it, ++s_it) - h += (*o_it)*(*p_it - *s_it)/orth_length; - h = fabs(h); - // Is the center inside the box? - bool inside_the_box = true; - for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) - if (*c_it > 1.0 || *c_it < -1.0) - { - inside_the_box = false; break; - } - if (inside_the_box && h/r < theta0) - return true; - if (!t.is_infinite(v1)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - if (!t.is_infinite(v2)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && 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; - } - */ - if (!t.is_infinite(v1)) - { - std::vector<FT> coords; - Point_d p = v1->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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - if (!t.is_infinite(v2)) - { - std::vector<FT> coords; - Point_d p = v2->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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - -/** Auxillary recursive function to check if the point p violates the protection of the cell c and - * if there is a violation of an eventual new cell - * - * p is the point to insert - * t is the current triangulation - * c is the current cell (simplex) - * parent_cell is the parent cell (simplex) - * index is the index of the facet between c and parent_cell from parent_cell's point of view - * D is the dimension of the triangulation - * delta is the protection constant - * marked_cells is the vector of all visited cells containing p in their circumscribed ball - * power_protection is true iff you are working with delta-power protection - * - * OUT: true iff inserting p hasn't produced any violation so far - */ - -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, bool power_protection, FT theta0) -{ - 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - // if the new point is inside the circumscribing ball : continue violation searching on neighbours - //if (dist2 < r*r) - //if (dist2 < (5*r+delta)*(5*r+delta)) - if (dist2 < r*r) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, power_protection, theta0)) - 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 = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, 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) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - // If we work with power protection, we just ignore any conflicts - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - if (p_is_inside) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, power_protection, theta0)) - 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); - for (int i = 0; i < D+1; ++i) - if (i != index) - if (!t.is_infinite(parent_cell->vertex(i))) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) - return true; - } - } - //c->tds_data().clear_visited(); - //marked_cells.pop_back(); - return false; -} - -/** Checks if inserting the point p in t will make conflicts - * - * p is the point to insert - * t is the current triangulation - * D is the dimension of triangulation - * delta is the protection constant - * power_protection is true iff you are working with delta-power protection - * OUT: true iff inserting p produces a violation of delta-protection. - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) -{ - 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; - 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, power_protection, theta0); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - -////////////////////////////////////////////////////////////////////// -// INITIALIZATION -////////////////////////////////////////////////////////////////////// - -void initialize(Search_Tree& W, Delaunay& t, int D, int width, bool torus) -{ - if (!torus) - std::cout << "Non-toric case is not supported\n"; - else - { - if (D == 2) - { - FT stepx = 2.0/width; - FT stepy = sqrt(3)/width; - for (int i = 0; i < width; ++i) - for (int j = 0; j < floor(2*width/sqrt(3)); ++j) - { - insert_delaunay_landmark_with_copies(Point_d(step*i,)) - } - } - else (D == 3) - { - - } - else std::cout << "T^d with d>3 not supported"; - } -} - -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// -//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// - -/////////////////////////////////////////////////////////////////////// -// LANDMARK CHOICE PROCEDURE AS IN PAPER -/////////////////////////////////////////////////////////////////////// - -/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. - * - * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d - * IN: nbP is the size of W - * OUT: landmarks is the output vector for the points - * OUT: landmarks_ind is the output vector for the indices of the selected points in W - * IN: delta is the constant of protection - * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -template<class Search_Tree> -void protected_delaunay_refinement(Search_Tree& W, int nbP, Point_Vector& landmarks, FT delta, bool torus, bool power_protection, FT theta0) -{ - bool return_ = true; - unsigned D = W[0].size(); - Torus_distance td; - Euclidean_distance ed; - Delaunay_triangulation t(D); - CGAL::Random rand; - int landmark_count = 0; - //std::list<int> index_list; - // shuffle the list of indexes (via a vector) - // { - // std::vector<int> 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)); - // //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); - // } - if (torus) - if (D == 2) - // \T^2 - { - for (int i = 0; i < 4; ++i) - for (int j = 0; j < 2; ++j) - { - W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0}); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); - index_list.pop_front(); - W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j}); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus); - index_list.pop_front(); - } - } - else if (D == 3) - { - - } - //std::cout << "No torus starter available for dim>2\n"; - 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, power_protection, theta0)) - { - // If no conflicts then insert in every copy of T^3 - - insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus); - if (return_) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - index_list.erase(list_it++); - /* - // PIECE OF CODE FOR DEBUGGING PURPOSES - - Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - if (triangulation_is_protected(t, delta)) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - { //THAT'S WHERE SOMETHING'S WRONG - t.remove(inserted_v); - landmarks_ind.pop_back(); - landmark_count--; - write_delaunay_mesh(t, W[*list_it], is2d); - is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore - } - */ - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - list_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //if (list_it != index_list.end()) - // write_delaunay_mesh(t, W[*list_it], is2d); - } - fill_landmarks(W, landmarks, landmarks_ind, torus); - 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, W[0], is2d); - //std::cout << t << std::endl; -} - -#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.h b/src/Witness_complex/example/protected_sets/protected_sets_paper.h deleted file mode 100644 index 61fcc75b..00000000 --- a/src/Witness_complex/example/protected_sets/protected_sets_paper.h +++ /dev/null @@ -1,917 +0,0 @@ -#ifndef PROTECTED_SETS_H -#define PROTECTED_SETS_H - -#include <algorithm> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.h> -#include <CGAL/Kernel_d/Hyperplane_d.h> -#include <CGAL/Kernel_d/Vector_d.h> - -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Fuzzy_sphere.h> - -#include <boost/heap/fibonacci_heap.hpp> -#include <boost/heap/policies.hpp> - -#include "output_tikz.h" - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -typedef K::Line_d Line_d; -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 K::Sphere_d Sphere_d; -typedef K::Hyperplane_d Hyperplane_d; - -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; -typedef Delaunay_triangulation::Facet Facet; -typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; -typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; - -typedef std::vector<Point_d> Point_Vector; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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 CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - - -FT _sfty = pow(10,-14); - -bool experiment1, experiment2 = false; - -/* Experiment 1: epsilon as function on time **********************/ -std::vector<FT> eps_vector; - -/* Experiment 2: R/epsilon on delta *******************************/ -std::vector<FT> epsratio_vector; - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// AUXILLARY FUNCTIONS -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well - * - * W is the initial point vector - * chosen_landmark is the index of the chosen point in W - * landmarks_ind is the vector of indices of already chosen points in W - * delaunay is the Delaunay triangulation - * landmark_count is the current number of chosen vertices - * torus is true iff you are working on a flat torus [-1,1]^d - * OUT: Vertex handle to the newly inserted point - */ -Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) -{ - if (!torus) - { - Delaunay_vertex v =delaunay.insert(W[chosen_landmark]); - landmarks_ind.push_back(chosen_landmark); - landmark_count++; - return v; - } - else - { - int D = W[0].size(); - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> 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; - } - if (i == nb_cells/2) - v = delaunay.insert(point); //v = center point - else - delaunay.insert(point); - } - landmarks_ind.push_back(chosen_landmark); - landmark_count++; - return v; - } -} - -/** Small check if the vertex v is in the full cell fc - */ - -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; -} - -/** Fill chosen point vector from indices with copies if you are working on a flat torus - * - * IN: W is the point vector - * OUT: landmarks is the output vector - * IN: landmarks_ind is the vector of indices - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) -{ - if (!torus) - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); - else - { - 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); - } - } -} - -/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices - * - * IN: t is the Delaunay triangulation - * OUT: full_cells is the output vector - */ - -void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) -{ - // 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) - { - 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; - 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); - } -} - -bool sphere_intersects_cube(Point_d& c, FT r) -{ - bool in_cube = true; - // int i = 0, D = p.size(); - for (auto xi = c.cartesian_begin(); xi != c.cartesian_end(); ++xi) - // if ((*xi < 1.0 || *xi > -1.0) && - // (*xi-r < 1.0 || *xi-r > -1.0) && - // (*xi+r < 1.0 || *xi+r > -1.0)) - - if ((*xi-r < -1.0 && *xi+r < -1.0) || - (*xi-r > 1.0 && *xi+r > 1.0 )) - { - in_cube = false; break; - } - return in_cube; -} - -/** Recursive function for checking if the simplex is good, - * meaning it does not contain a k-face, which is not theta0^(k-1) thick - */ - -bool is_theta0_good(std::vector<Point_d>& vertices, FT theta0) -{ - if (theta0 > 1) - { - std::cout << "Warning! theta0 is set > 1\n"; - return false; - } - int D = vertices.size()-1; - if (D <= 1) - return true; // Edges are always good - //******** Circumscribed sphere - Euclidean_distance ed; - Sphere_d cs(vertices.begin(), vertices.end()); - FT r = sqrt(cs.squared_radius()); - for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it) - { - std::vector<Point_d> facet; - for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it) - if (f_it != v_it) - facet.push_back(*f_it); - // Compute the altitude - - if (vertices[0].size() == 3 && D == 2) - { - //Vector_d l = facet[0] - facet[1]; - FT orth_length2 = ed.transformed_distance(facet[0],facet[1]); - K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it; - FT h = 0; - // Scalar product = <sp,l> - FT scalar = 0; - for (p_it = v_it->cartesian_begin(), - s_it = facet[0].cartesian_begin(), - l_it = facet[1].cartesian_begin(); - p_it != v_it->cartesian_end(); - ++l_it, ++p_it, ++s_it) - scalar += (*l_it - *s_it)*(*p_it - *s_it); - // Gram-Schmidt for one vector - for (p_it = v_it->cartesian_begin(), - s_it = facet[0].cartesian_begin(), - l_it = facet[1].cartesian_begin(); - p_it != v_it->cartesian_end(); - ++l_it, ++p_it, ++s_it) - { - FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2; - h += hx*hx; - } - h = sqrt(h); - - if (h/(2*r) < pow(theta0, D-1)) - return false; - if (!is_theta0_good(facet, theta0)) - return false; - } - else - { - Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it); - Vector_d orth_tau = tau_h.orthogonal_vector(); - FT orth_length = sqrt(orth_tau.squared_length()); - K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; - FT h = 0; - for (o_it = orth_tau.cartesian_begin(), - p_it = v_it->cartesian_begin(), - s_it = (facet.begin())->cartesian_begin(); - o_it != orth_tau.cartesian_end(); - ++o_it, ++p_it, ++s_it) - h += (*o_it)*(*p_it - *s_it)/orth_length; - h = fabs(h); - if (h/(2*r) < pow(theta0, D-1)) - return false; - if (!is_theta0_good(facet, theta0)) - return false; - } - } - return true; -} - - -//////////////////////////////////////////////////////////////////////////////////////////////////////////// -// IS VIOLATED TEST -//////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Check if a newly created cell is protected from old vertices - * - * t is the Delaunay triangulation - * vertices is the vector containing the point to insert and a facet f in t - * v1 is the vertex of t, such that f and v1 form a simplex - * v2 is the vertex of t, such that f and v2 form another simplex - * delta is the protection constant - * power_protection is true iff the delta-power protection is used - */ - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0) -{ - assert(vertices.size() == vertices[0].size() || - vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 - assert(v1 != v2); - if (vertices.size() == vertices[0].size() + 1) - // 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - } - */ - // Check if the simplex is theta0-good - if (!is_theta0_good(vertices, theta0)) - return true; - // Is the center inside the box? (only Euclidean case) - // if (!torus) - // { - // bool inside_the_box = true; - // for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) - // if (*c_it > 1.0 || *c_it < -1.0) - // { - // inside_the_box = false; break; - // } - // if (inside_the_box && h/r < theta0) - // return true; - // } - // Check the two vertices (if not infinite) - if (!t.is_infinite(v1)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - if (!t.is_infinite(v2)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && 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; - } - */ - if (!t.is_infinite(v1)) - { - std::vector<FT> coords; - Point_d p = v1->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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - if (!t.is_infinite(v2)) - { - std::vector<FT> coords; - Point_d p = v2->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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - -/** Auxillary recursive function to check if the point p violates the protection of the cell c and - * if there is a violation of an eventual new cell - * - * p is the point to insert - * t is the current triangulation - * c is the current cell (simplex) - * parent_cell is the parent cell (simplex) - * index is the index of the facet between c and parent_cell from parent_cell's point of view - * D is the dimension of the triangulation - * delta is the protection constant - * marked_cells is the vector of all visited cells containing p in their circumscribed ball - * power_protection is true iff you are working with delta-power protection - * - * OUT: true iff inserting p hasn't produced any violation so far - */ - -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, bool power_protection, FT theta0) -{ - 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - // if the new point is inside the circumscribing ball : continue violation searching on neighbours - //if (dist2 < r*r) - //if (dist2 < (5*r+delta)*(5*r+delta)) - if (dist2 < r*r) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, power_protection, theta0)) - 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 = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, 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) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - // If we work with power protection, we just ignore any conflicts - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - if (p_is_inside) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, power_protection, theta0)) - 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); - for (int i = 0; i < D+1; ++i) - if (i != index) - if (!t.is_infinite(parent_cell->vertex(i))) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) - //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) - return true; - } - } - //c->tds_data().clear_visited(); - //marked_cells.pop_back(); - return false; -} - -/** Checks if inserting the point p in t will make conflicts - * - * p is the point to insert - * t is the current triangulation - * D is the dimension of triangulation - * delta is the protection constant - * power_protection is true iff you are working with delta-power protection - * OUT: true iff inserting p produces a violation of delta-protection. - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0) -{ - 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; - std::vector<Full_cell_handle> marked_cells; - //c = t.locate(p, lt, f, ft, v); - c = t.locate(p); - bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - - -//////////////////////////////////////////////////////////////////////// -// INITIALIZATION -//////////////////////////////////////////////////////////////////////// - -// Query for a sphere near a cite in all copies of a torus -// OUT points_inside -void torus_search(Tree& treeW, int D, Point_d cite, FT r, std::vector<int>& points_inside) -{ - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> cite_copy; - int cell_i = i; - for (int l = 0; l < D; ++l) - { - cite_copy.push_back(cite[l] + 2.0*(cell_i%3-1)); - cell_i /= 3; - } - Fuzzy_sphere fs(cite_copy, r, 0, treeW.traits()); - treeW.search(std::insert_iterator<std::vector<int>>(points_inside, points_inside.end()), fs); - } -} - - -void initialize_torus(Point_Vector& W, Tree& treeW, Delaunay_triangulation& t, FT epsilon, std::vector<int>& landmarks_ind, int& landmark_count) -{ - int D = W[0].size(); - if (D == 2) - { - int xw = 6, yw = 4; - // Triangular lattice close to regular triangles h=0.866a ~ 0.875a : 48p - for (int i = 0; i < xw; ++i) - for (int j = 0; j < yw; ++j) - { - Point_d cite1(std::vector<FT>{2.0/xw*i, 1.0/yw*j}); - std::vector<int> points_inside; - torus_search(treeW, D, cite1, epsilon, points_inside); - assert(points_inside.size() > 0); - insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), - landmarks_ind, t, landmark_count, true); - Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 1.0/yw*(j+0.5)}); - points_inside.clear(); - torus_search(treeW, D, cite2, epsilon, points_inside); - assert(points_inside.size() > 0); - insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), - landmarks_ind, t, landmark_count, true); - } - } - else if (D == 3) - { - int wd = 3; - // Body-centered cubic lattice : 54p - for (int i = 0; i < wd; ++i) - for (int j = 0; j < wd; ++j) - for (int k = 0; k < wd; ++k) - { - Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k}); - std::vector<int> points_inside; - torus_search(treeW, D, cite1, epsilon, points_inside); - assert(points_inside.size() > 0); - insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), - landmarks_ind, t, landmark_count, true); - Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)}); - points_inside.clear(); - torus_search(treeW, D, cite2, epsilon, points_inside); - assert(points_inside.size() > 0); - insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), - landmarks_ind, t, landmark_count, true); - } - } -} - -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// -//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// - -// Struct for R_max_heap elements - -struct R_max_handle -{ - FT value; - Point_d center; - - R_max_handle(FT value_, Point_d c): value(value_), center(c) - {} -}; - -struct R_max_compare -{ - bool operator()(const R_max_handle& rmh1, const R_max_handle& rmh2) const - { - return rmh1.value < rmh2.value; - } -}; - -// typedef boost::heap::fibonacci_heap<R_max_handle, boost::heap::compare<R_max_compare>> Heap; - -// void make_heap(Delaunay_triangulation& t, Heap& R_max_heap) -// { -// R_max_heap.clear(); -// 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(); -// FT r = sqrt(cs.squared_radius()); -// // A ball is in the heap, if it intersects the cube -// bool accepted = sphere_intersects_cube(csc, sqrt(r)); -// if (!accepted) -// continue; -// R_max_heap.push(R_max_handle(r, fc_it, csc)); -// } -// } - -////////////////////////////////////////////////////////////////////////////////////////////////////////// -// SAMPLING RADIUS -////////////////////////////////////////////////////////////////////////////////////////////////////////// - -R_max_handle sampling_radius(Delaunay_triangulation& t) -{ - FT epsilon2 = 0; - Point_d final_center; - Point_d control_point; - 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; - final_center = csc; - control_point = (*vertices.begin()); - } - } - return R_max_handle(sqrt(epsilon2), final_center); -} - -/////////////////////////////////////////////////////////////////////// -// LANDMARK CHOICE PROCEDURE -/////////////////////////////////////////////////////////////////////// - -/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. - * - * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d - * IN: nbP is the size of W - * OUT: landmarks is the output vector for the points - * OUT: landmarks_ind is the output vector for the indices of the selected points in W - * IN: delta is the constant of protection - * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void protected_delaunay(Point_Vector& W, - //Point_Vector& landmarks, - std::vector<int>& landmarks_ind, - FT delta, - FT epsilon, - FT alpha, - FT theta0, - //std::vector<std::vector<int>>& full_cells, - bool torus, - bool power_protection - ) -{ - //bool return_ = true; - unsigned D = W[0].size(); - int nbP = W.size(); - Torus_distance td; - Euclidean_distance ed; - Delaunay_triangulation t(D); - CGAL::Random rand; - int landmark_count = 0; - std::list<int> index_list; - //****************** Kd Tree W - STraits traits(&(W[0])); - Tree treeW(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(nbP), - typename Tree::Splitter(), - traits); - // shuffle the list of indexes (via a vector) - { - std::vector<int> 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)); - //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); - } - //******************** Initialize point set - if (!torus) - 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, torus); - index_list.pop_front(); - } - else - initialize_torus(W, treeW, t, epsilon, landmarks_ind, landmark_count); - //std::cout << "Size of treeW: " << treeW.size() << "\n"; - //std::cout << "Size of t: " << t.number_of_vertices() << "\n"; - //******************* Initialize heap for R_max - //Heap R_max_heap; - //make_heap(t, R_max_heap); - - - R_max_handle rh = sampling_radius(t); - FT epsilon0 = rh.value; - if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); - //******************** Iterative algorithm - std::vector<int> candidate_points; - torus_search(treeW, D, - rh.center, - alpha*rh.value, - candidate_points); - std::list<int>::iterator list_it; - std::vector<int>::iterator cp_it = candidate_points.begin(); - while (cp_it != candidate_points.end()) - { - if (!is_violating_protection(W[*cp_it], t, D, delta, power_protection, theta0)) - { - insert_delaunay_landmark_with_copies(W, *cp_it, landmarks_ind, t, landmark_count, torus); - //make_heap(t, R_max_heap); - rh = sampling_radius(t); - if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); - //std::cout << "rhvalue = " << rh.value << "\n"; - //std::cout << "D = " << - candidate_points.clear(); - torus_search(treeW, D, - rh.center, - alpha*rh.value, - candidate_points); - /* - // PIECE OF CODE FOR DEBUGGING PURPOSES - - Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - if (triangulation_is_protected(t, delta)) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - { //THAT'S WHERE SOMETHING'S WRONG - t.remove(inserted_v); - landmarks_ind.pop_back(); - landmark_count--; - write_delaunay_mesh(t, W[*list_it], is2d); - is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore - } - */ - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - cp_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //if (list_it != index_list.end()) - // write_delaunay_mesh(t, W[*list_it], is2d); - } - if (experiment2) epsratio_vector.push_back(rh.value/epsilon0); - std::cout << "The iteration ended when cp_count = " << candidate_points.size() << "\n"; - std::cout << "alphaRmax = " << alpha*rh.value << "\n"; - std::cout << "epsilon' = " << rh.value << "\n"; - std::cout << "nbL = " << landmarks_ind.size() << "\n"; - //fill_landmarks(W, landmarks, landmarks_ind, torus); - //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, W[0], is2d); - //std::cout << t << std::endl; -} - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Series of experiments -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -void start_experiments(Point_Vector& W, FT theta0, std::vector<int>& landmarks_ind, FT epsilon) -{ - // Experiment 1 - experiment1 = true; - protected_delaunay(W, landmarks_ind, 0.1*epsilon, epsilon, 0.5, 0, true, true); - write_tikz_plot(eps_vector,"epstime.tikz"); - experiment1 = false; - - // Experiment 2 - // experiment2 = true; - // for (FT delta = 0; delta < epsilon; delta += 0.1*epsilon) - // protected_delaunay(W, landmarks_ind, delta, epsilon, 0.5, 0, true, true); - // write_tikz_plot(epsratio_vector,"epsratio_delta.tikz"); - // experiment2 = false; - -} - -#endif diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper2.h b/src/Witness_complex/example/protected_sets/protected_sets_paper2.h deleted file mode 100644 index 04b5e3bc..00000000 --- a/src/Witness_complex/example/protected_sets/protected_sets_paper2.h +++ /dev/null @@ -1,1384 +0,0 @@ -#ifndef PROTECTED_SETS_H -#define PROTECTED_SETS_H - -#include <algorithm> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.h> -#include <CGAL/Kernel_d/Hyperplane_d.h> -#include <CGAL/Kernel_d/Vector_d.h> - -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Fuzzy_sphere.h> - -#include <boost/heap/fibonacci_heap.hpp> -#include <boost/heap/policies.hpp> - -#include "output_tikz.h" -#include "../output.h" -#include "../generators.h" - -#include <CGAL/point_generators_d.h> - - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -typedef K::Line_d Line_d; -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 K::Sphere_d Sphere_d; -typedef K::Hyperplane_d Hyperplane_d; - -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; -typedef Delaunay_triangulation::Facet Facet; -typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex; -typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; - -typedef std::vector<Point_d> Point_Vector; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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 CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; - - -FT _sfty = pow(10,-14); - -bool experiment1, experiment2, experiment3, experiment5 = false; - -/* Experiment 1: epsilon as function on time **********************/ -std::vector<FT> eps_vector; - -/* Experiment 2: R/epsilon on alpha *******************************/ -std::vector<FT> epsratio_vector; -std::vector<FT> epsslope_vector; - -/* Experiment 3: theta on delta ***********************************/ -std::vector<FT> thetamin_vector; FT curr_theta; -std::vector<FT> gammamin_vector; - -/* Statistical data ***********************************************/ -int refused_case1, refused_case2, refused_bad, refused_centers1, refused_centers2; - -void initialize_statistics() -{ - refused_case1 = 0; - refused_case2 = 0; - refused_bad = 0; - refused_centers1 = 0; - refused_centers2 = 0; -} - -void print_statistics() -{ - std::cout << " * Old simplex not protected: " << refused_case1 << "\n"; - std::cout << " * New simplex not protected: " << refused_case2 << "\n"; - std::cout << " * New simplex not good: " << refused_bad << "\n"; - std::cout << " * New-old centers too close: " << refused_centers1 << "\n"; - std::cout << " * New-new centers too close: " << refused_centers2 << "\n"; -} - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// AUXILLARY FUNCTIONS -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well - * - * W is the initial point vector - * chosen_landmark is the index of the chosen point in W - * landmarks_ind is the vector of indices of already chosen points in W - * delaunay is the Delaunay triangulation - * landmark_count is the current number of chosen vertices - * torus is true iff you are working on a flat torus [-1,1]^d - * OUT: Vertex handle to the newly inserted point - */ -Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus) -{ - if (!torus) - { - Delaunay_vertex v =delaunay.insert(W[chosen_landmark]); - landmarks_ind.push_back(chosen_landmark); - landmark_count++; - return v; - } - else - { - int D = W[0].size(); - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> 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; - } - if (i == nb_cells/2) - v = delaunay.insert(point); //v = center point - else - delaunay.insert(point); - } - landmarks_ind.push_back(chosen_landmark); - landmark_count++; - return v; - } -} - -/** Small check if the vertex v is in the full cell fc - */ - -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; -} - -/** Fill chosen point vector from indices with copies if you are working on a flat torus - * - * IN: W is the point vector - * OUT: landmarks is the output vector - * IN: landmarks_ind is the vector of indices - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus) -{ - if (!torus) - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); - else - { - 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); - } - } -} - -/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices - * - * IN: t is the Delaunay triangulation - * OUT: full_cells is the output vector - */ - -void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells) -{ - // 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) - { - 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; - 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); - } -} - -bool sphere_intersects_cube(Point_d& c, FT r) -{ - bool in_cube = true; - // int i = 0, D = p.size(); - for (auto xi = c.cartesian_begin(); xi != c.cartesian_end(); ++xi) - // if ((*xi < 1.0 || *xi > -1.0) && - // (*xi-r < 1.0 || *xi-r > -1.0) && - // (*xi+r < 1.0 || *xi+r > -1.0)) - - if ((*xi-r < -1.0 && *xi+r < -1.0) || - (*xi-r > 1.0 && *xi+r > 1.0 )) - { - in_cube = false; break; - } - return in_cube; -} - -/** Recursive function for checking if the simplex is good, - * meaning it does not contain a k-face, which is not theta0^(k-1) thick - */ - -bool is_theta0_good(std::vector<Point_d>& vertices, FT theta0) -{ - if (theta0 > 1) - { - std::cout << "Warning! theta0 is set > 1\n"; - return false; - } - int D = vertices.size()-1; - if (D <= 1) - return true; // Edges are always good - //******** Circumscribed sphere - Euclidean_distance ed; - Sphere_d cs(vertices.begin(), vertices.end()); - FT r = sqrt(cs.squared_radius()); - for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it) - { - std::vector<Point_d> facet; - for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it) - if (f_it != v_it) - facet.push_back(*f_it); - // Compute the altitude - - if (vertices[0].size() == 3 && D == 2) - { - //Vector_d l = facet[0] - facet[1]; - FT orth_length2 = ed.transformed_distance(facet[0],facet[1]); - K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it; - FT h = 0; - // Scalar product = <sp,l> - FT scalar = 0; - for (p_it = v_it->cartesian_begin(), - s_it = facet[0].cartesian_begin(), - l_it = facet[1].cartesian_begin(); - p_it != v_it->cartesian_end(); - ++l_it, ++p_it, ++s_it) - scalar += (*l_it - *s_it)*(*p_it - *s_it); - // Gram-Schmidt for one vector - for (p_it = v_it->cartesian_begin(), - s_it = facet[0].cartesian_begin(), - l_it = facet[1].cartesian_begin(); - p_it != v_it->cartesian_end(); - ++l_it, ++p_it, ++s_it) - { - FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2; - h += hx*hx; - } - h = sqrt(h); - - if (h/(2*r) < pow(theta0, D-1)) - return false; - if (!is_theta0_good(facet, theta0)) - return false; - } - else - { - Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it); - Vector_d orth_tau = tau_h.orthogonal_vector(); - FT orth_length = sqrt(orth_tau.squared_length()); - K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; - FT h = 0; - for (o_it = orth_tau.cartesian_begin(), - p_it = v_it->cartesian_begin(), - s_it = (facet.begin())->cartesian_begin(); - o_it != orth_tau.cartesian_end(); - ++o_it, ++p_it, ++s_it) - h += (*o_it)*(*p_it - *s_it)/orth_length; - h = fabs(h); - if (experiment3 && thetamin_vector[thetamin_vector.size()-1] > pow(h/(2*r), 1.0/(D-1))) - { - thetamin_vector[thetamin_vector.size()-1] = pow(h/(2*r), 1.0/(D-1)); - //std::cout << "theta=" << h/(2*r) << ", "; - } - if (h/(2*r) < pow(theta0, D-1)) - return false; - if (!is_theta0_good(facet, theta0)) - return false; - } - } - return true; -} - -/** Recursive function for checking the goodness of a simplex, - * meaning it does not contain a k-face, which is not theta0^(k-1) thick - */ - -FT theta(std::vector<Point_d>& vertices) -{ - FT curr_value = 1.0; - int D = vertices.size()-1; - if (D <= 1) - return 1; // Edges are always good - //******** Circumscribed sphere - Euclidean_distance ed; - Sphere_d cs(vertices.begin(), vertices.end()); - FT r = sqrt(cs.squared_radius()); - for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it) - { - std::vector<Point_d> facet; - for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it) - if (f_it != v_it) - facet.push_back(*f_it); - // Compute the altitude - curr_value = std::min(curr_value, theta(facet)); // Check the corresponding facet - if (vertices[0].size() == 3 && D == 2) - { - //Vector_d l = facet[0] - facet[1]; - FT orth_length2 = ed.transformed_distance(facet[0],facet[1]); - K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it; - FT h = 0; - // Scalar product = <sp,l> - FT scalar = 0; - for (p_it = v_it->cartesian_begin(), - s_it = facet[0].cartesian_begin(), - l_it = facet[1].cartesian_begin(); - p_it != v_it->cartesian_end(); - ++l_it, ++p_it, ++s_it) - scalar += (*l_it - *s_it)*(*p_it - *s_it); - // Gram-Schmidt for one vector - for (p_it = v_it->cartesian_begin(), - s_it = facet[0].cartesian_begin(), - l_it = facet[1].cartesian_begin(); - p_it != v_it->cartesian_end(); - ++l_it, ++p_it, ++s_it) - { - FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2; - h += hx*hx; - } - h = sqrt(h); - curr_value = std::min(curr_value, std::pow(h/(2*r), 1.0/(D-1))); - } - else - { - Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it); - Vector_d orth_tau = tau_h.orthogonal_vector(); - FT orth_length = sqrt(orth_tau.squared_length()); - K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it; - FT h = 0; - for (o_it = orth_tau.cartesian_begin(), - p_it = v_it->cartesian_begin(), - s_it = (facet.begin())->cartesian_begin(); - o_it != orth_tau.cartesian_end(); - ++o_it, ++p_it, ++s_it) - h += (*o_it)*(*p_it - *s_it)/orth_length; - h = fabs(h); - curr_value = std::min(curr_value, pow(h/(2*r), 1.0/(D-1))); - } - } - return curr_value; -} - -// Doubling in a way 1->2->5->10 -void double_round(int& i) -{ - FT order10 = pow(10,std::floor(std::log10(i))); - int digit = std::floor( i / order10); - std::cout << digit; - if (digit == 1) - i *= 2; - else if (digit == 2) - i = 5*i/2; - else if (digit == 5) - i *= 2; - else - std::cout << "digit not correct. digit = " << digit << std::endl; -} - -//////////////////////////////////////////////////////////////////////////////////////////////////////////// -// IS VIOLATED TEST -//////////////////////////////////////////////////////////////////////////////////////////////////////////// - -/** Check if a newly created cell is protected from old vertices - * - * t is the Delaunay triangulation - * vertices is the vector containing the point to insert and a facet f in t - * v1 is the vertex of t, such that f and v1 form a simplex - * v2 is the vertex of t, such that f and v2 form another simplex - * delta is the protection constant - * power_protection is true iff the delta-power protection is used - */ - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta0, bool power_protection, FT theta0, FT gamma0) -{ - assert(vertices.size() == vertices[0].size() || - vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1 - assert(v1 != v2); - if (vertices.size() == vertices[0].size() + 1) - // 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta)) - return true; - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) - return true; - } - } - */ - // Is the center inside the box? (only Euclidean case) - // if (!torus) - // { - // bool inside_the_box = true; - // for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it) - // if (*c_it > 1.0 || *c_it < -1.0) - // { - // inside_the_box = false; break; - // } - // if (inside_the_box && h/r < theta0) - // return true; - // } - // Check the two vertices (if not infinite) - if (!t.is_infinite(v1)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0)) - { refused_case2++; return true;} - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+r*r*delta0*delta0) - { refused_case2++; return true;} - // Check if the centers are not too close - std::vector<Point_d> sigma(vertices); - sigma[0] = v1->point(); - Sphere_d cs_sigma(sigma.begin(), sigma.end()); - Point_d csc_sigma = cs_sigma.center(); - FT r_sigma = sqrt(cs_sigma.squared_radius()); - FT dcc = sqrt(Euclidean_distance().transformed_distance(center_cs, csc_sigma)); - if (experiment3 && dcc/r < gammamin_vector[gammamin_vector.size()-1]) - gammamin_vector[gammamin_vector.size()-1] = dcc/r; - if (experiment3 && dcc/r_sigma < gammamin_vector[gammamin_vector.size()-1]) - gammamin_vector[gammamin_vector.size()-1] = dcc/r_sigma; - if (dcc < r*gamma0 || dcc < r_sigma*gamma0) - { refused_centers1++; return true; } - } - if (!t.is_infinite(v2)) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point()); - if (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0)) - { refused_case2++; return true;} - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+r*r*delta0*delta0) - { refused_case2++; return true;} - // Check if the centers are not too close - std::vector<Point_d> sigma(vertices); - sigma[0] = v2->point(); - Sphere_d cs_sigma(sigma.begin(), sigma.end()); - Point_d csc_sigma = cs_sigma.center(); - FT r_sigma = sqrt(cs_sigma.squared_radius()); - FT dcc = sqrt(Euclidean_distance().transformed_distance(center_cs, csc_sigma)); - if (experiment3 && dcc/r < gammamin_vector[gammamin_vector.size()-1]) - gammamin_vector[gammamin_vector.size()-1] = dcc/r; - if (experiment3 && dcc/r_sigma < gammamin_vector[gammamin_vector.size()-1]) - gammamin_vector[gammamin_vector.size()-1] = dcc/r_sigma; - if (dcc < r*gamma0 || dcc < r_sigma*gamma0) - { refused_centers1++; return true; } - } - // Check if the simplex is theta0-good - if (!is_theta0_good(vertices, theta0)) - { refused_bad++; 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; - } - */ - if (!t.is_infinite(v1)) - { - std::vector<FT> coords; - Point_d p = v1->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) * delta0 / 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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - if (!t.is_infinite(v2)) - { - std::vector<FT> coords; - Point_d p = v2->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) * delta0 / 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 (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - -/** Auxillary recursive function to check if the point p violates the protection of the cell c and - * if there is a violation of an eventual new cell - * - * p is the point to insert - * t is the current triangulation - * c is the current cell (simplex) - * parent_cell is the parent cell (simplex) - * index is the index of the facet between c and parent_cell from parent_cell's point of view - * D is the dimension of the triangulation - * delta is the protection constant - * marked_cells is the vector of all visited cells containing p in their circumscribed ball - * power_protection is true iff you are working with delta-power protection - * - * OUT: true iff inserting p hasn't produced any violation so far - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta0, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0, FT gamma0) -{ - 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 (!power_protection) - if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0)) - { refused_case1++; return true;} - if (power_protection) - if (dist2 >= r*r-_sfty && dist2 <= r*r+delta0*delta0*r*r) - { refused_case1++; return true;} - // if the new point is inside the circumscribing ball : continue violation searching on neighbours - //if (dist2 < r*r) - //if (dist2 < (5*r+delta)*(5*r+delta)) - if (dist2 < r*r) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, delta0, marked_cells, power_protection, theta0, gamma0)) - 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 = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta0, power_protection, theta0, gamma0)) - //if (new_cell_is_violated(t, vertices, 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) * delta0 / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - // If we work with power protection, we just ignore any conflicts - if (!power_protection && !p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - if (p_is_inside) - { - c->tds_data().mark_visited(); - marked_cells.push_back(c); - 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, delta0, marked_cells, power_protection, theta0, gamma0)) - 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); - for (int i = 0; i < D+1; ++i) - if (i != index) - if (!t.is_infinite(parent_cell->vertex(i))) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check = t.infinite_vertex(); - 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, vertex_to_check, parent_cell->vertex(index), delta0, power_protection, theta0, gamma0)) - //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) - return true; - } - } - //c->tds_data().clear_visited(); - //marked_cells.pop_back(); - return false; -} - -/** Checks if inserting the point p in t will make conflicts - * - * p is the point to insert - * t is the current triangulation - * D is the dimension of triangulation - * delta is the protection constant - * power_protection is true iff you are working with delta-power protection - * OUT: true iff inserting p produces a violation of delta-protection. - */ - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta0, bool power_protection, FT theta0, FT gamma0) -{ - 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; - std::vector<Full_cell_handle> marked_cells; - //c = t.locate(p, lt, f, ft, v); - c = t.locate(p); - bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta0, marked_cells, power_protection, theta0, gamma0); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - - -//////////////////////////////////////////////////////////////////////// -// INITIALIZATION -//////////////////////////////////////////////////////////////////////// - -// Query for a sphere near a cite in all copies of a torus -// OUT points_inside -void torus_search(Tree& treeW, int D, Point_d cite, FT r, std::vector<int>& points_inside) -{ - int nb_cells = pow(3, D); - Delaunay_vertex v; - for (int i = 0; i < nb_cells; ++i) - { - std::vector<FT> cite_copy; - int cell_i = i; - for (int l = 0; l < D; ++l) - { - cite_copy.push_back(cite[l] + 2.0*(cell_i%3-1)); - cell_i /= 3; - } - Fuzzy_sphere fs(cite_copy, r, 0, treeW.traits()); - treeW.search(std::insert_iterator<std::vector<int>>(points_inside, points_inside.end()), fs); - } -} - - -void initialize_torus(Point_Vector& W, Tree& treeW, Delaunay_triangulation& t, FT epsilon, std::vector<int>& landmarks_ind, int& landmark_count, std::vector<bool>& point_taken) -{ - initialize_statistics(); - int D = W[0].size(); - if (D == 2) - { - int xw = 6, yw = 4; - // Triangular lattice close to regular triangles h=0.866a ~ 0.875a : 48p - for (int i = 0; i < xw; ++i) - for (int j = 0; j < yw; ++j) - { - Point_d cite1(std::vector<FT>{2.0/xw*i, 2.0/yw*j}); - std::vector<int> points_inside; - torus_search(treeW, D, cite1, epsilon, points_inside); - //std::cout << "i=" << i << ", j=" << j << " "; print_vector(points_inside); std::cout << "\n"; - std::vector<int>::iterator p_it = points_inside.begin(); - while (p_it != points_inside.end() && point_taken[*p_it]) - ++p_it; - assert(p_it != points_inside.end()); - //W[*p_it] = cite1; // debug purpose - insert_delaunay_landmark_with_copies(W, *p_it, - landmarks_ind, t, landmark_count, true); - point_taken[*p_it] = true; - - Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 2.0/yw*(j+0.5)}); - points_inside.clear(); - torus_search(treeW, D, cite2, epsilon, points_inside); - //std::cout << "i=" << i << ", j=" << j << " "; print_vector(points_inside); std::cout << "\n"; - p_it = points_inside.begin(); - while (p_it != points_inside.end() && point_taken[*p_it]) - ++p_it; - assert(p_it != points_inside.end()); - //W[*p_it] = cite2; // debug purpose - insert_delaunay_landmark_with_copies(W, *p_it, - landmarks_ind, t, landmark_count, true); - point_taken[*p_it] = true; - } - } - else if (D == 3) - { - int wd = 3; - // Body-centered cubic lattice : 54p - for (int i = 0; i < wd; ++i) - for (int j = 0; j < wd; ++j) - for (int k = 0; k < wd; ++k) - { - Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k}); - std::vector<int> points_inside; - torus_search(treeW, D, cite1, epsilon, points_inside); - std::vector<int>::iterator p_it = points_inside.begin(); - while (p_it != points_inside.end() && point_taken[*p_it]) - ++p_it; - assert(p_it != points_inside.end()); - insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), - landmarks_ind, t, landmark_count, true); - point_taken[*p_it] = true; - - Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)}); - points_inside.clear(); - torus_search(treeW, D, cite2, epsilon, points_inside); - p_it = points_inside.begin(); - while (p_it != points_inside.end() && point_taken[*p_it]) - ++p_it; - assert(p_it != points_inside.end()); - insert_delaunay_landmark_with_copies(W, *(points_inside.begin()), - landmarks_ind, t, landmark_count, true); - point_taken[*p_it] = true; - } - } - //write_mesh -} - -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// -//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!// -/////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////// - -// Struct for R_max_heap elements - -struct R_max_handle -{ - FT value; - Point_d center; - - R_max_handle(FT value_, Point_d c): value(value_), center(c) - {} -}; - -struct R_max_compare -{ - bool operator()(const R_max_handle& rmh1, const R_max_handle& rmh2) const - { - return rmh1.value < rmh2.value; - } -}; - -// typedef boost::heap::fibonacci_heap<R_max_handle, boost::heap::compare<R_max_compare>> Heap; - -// void make_heap(Delaunay_triangulation& t, Heap& R_max_heap) -// { -// R_max_heap.clear(); -// 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(); -// FT r = sqrt(cs.squared_radius()); -// // A ball is in the heap, if it intersects the cube -// bool accepted = sphere_intersects_cube(csc, sqrt(r)); -// if (!accepted) -// continue; -// R_max_heap.push(R_max_handle(r, fc_it, csc)); -// } -// } - -////////////////////////////////////////////////////////////////////////////////////////////////////////// -// SAMPLING RADIUS -////////////////////////////////////////////////////////////////////////////////////////////////////////// - -R_max_handle sampling_radius(Delaunay_triangulation& t) -{ - FT epsilon2 = 0; - Point_d final_center; - Point_d control_point; - 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; - final_center = csc; - control_point = (*vertices.begin()); - } - } - return R_max_handle(sqrt(epsilon2), final_center); -} - -FT sampling_fatness(Delaunay_triangulation& t) -{ - FT curr_theta = 1.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 theta_f = theta(vertices); - curr_theta = std::min(curr_theta, theta_f); - //std::cout << "theta(sigma) = " << theta_f << "\n"; - } - return curr_theta; -} - -// Generate an epsilon sample for a given epsilon -void generate_epsilon_sample_torus(Point_Vector& W, FT epsilon, int dim, Delaunay_triangulation& t) -{ - W.clear(); - t.clear(); - int point_count = 0; - std::vector<int> point_ind; - // std::vector<FT> coords; - FT curr_eps = 2*dim; - // Initialize - // for (int i = 0; i < dim; ++i) - // coords.push_back(-1); - // R_max_handle rmh(2*sqrt(dim), Point_d(coords)); - // int N = dim; std::floor(std::pow(1/epsilon,dim)); - // std::cout << N << "\n"; - typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; - Random_cube_iterator rp(dim, 1.0); - W.push_back(*rp++); - insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); - curr_eps = sampling_radius(t).value; - while (curr_eps > epsilon) - { - - W.push_back(*rp++); - insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); - - Point_d c = sampling_radius(t).center; - W.push_back(c); - insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); - curr_eps = sampling_radius(t).value; - - std::cout << "curr_eps = " << curr_eps << "\n"; - } - // Iterate and insert in a torus - // while (rmh.value > epsilon) - // { - // W.push_back(rmh.center); - // insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true); - // rmh = sampling_radius(t); - // //std::cout << rmh.value; - // } -} - -/////////////////////////////////////////////////////////////////////// -// LANDMARK CHOICE PROCEDURE -/////////////////////////////////////////////////////////////////////// - -/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call. - * - * IN: W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d - * IN: nbP is the size of W - * OUT: landmarks is the output vector for the points - * OUT: landmarks_ind is the output vector for the indices of the selected points in W - * IN: delta is the constant of protection - * OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation - * IN: torus is true iff you are working on a flat torus [-1,1]^d - */ - -void protected_delaunay(Point_Vector& W, - //Point_Vector& landmarks, - std::vector<int>& landmarks_ind, - FT alpha, - FT epsilon, - FT delta0, - FT theta0, - FT gamma0, - //std::vector<std::vector<int>>& full_cells, - bool torus, - bool power_protection - ) -{ - //bool return_ = true; - unsigned D = W[0].size(); - int nbP = W.size(); - //FT beta = 1/(1-alpha); - //FT Ad = pow((4*alpha + 8*beta)/alpha, D); - //FT theta0 = 1/Ad; - //FT delta0 = pow(1/Ad,D); - Torus_distance td; - Euclidean_distance ed; - Delaunay_triangulation t(D); - std::vector<bool> point_taken(nbP,false); - CGAL::Random rand; - int landmark_count = 0; - std::list<int> index_list; - //****************** Kd Tree W - STraits traits(&(W[0])); - Tree treeW(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(nbP), - typename Tree::Splitter(), - traits); - // shuffle the list of indexes (via a vector) - { - std::vector<int> 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)); - //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); - } - //******************** Initialize point set - if (!torus) - 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, torus); - index_list.pop_front(); - } - else - initialize_torus(W, treeW, t, epsilon, landmarks_ind, landmark_count, point_taken); - //std::cout << "Size of treeW: " << treeW.size() << "\n"; - //std::cout << "Size of t: " << t.number_of_vertices() << "\n"; - //******************* Initialize heap for R_max - //Heap R_max_heap; - //make_heap(t, R_max_heap); - - - R_max_handle rh = sampling_radius(t); - FT epsilon0 = rh.value; - if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); - //******************** Iterative algorithm - std::vector<int> candidate_points; - torus_search(treeW, D, - rh.center, - alpha*rh.value, - candidate_points); - std::list<int>::iterator list_it; - std::vector<int>::iterator cp_it = candidate_points.begin(); - while (cp_it != candidate_points.end()) - { - if (!point_taken[*cp_it] && !is_violating_protection(W[*cp_it], t, D, delta0, power_protection, theta0, gamma0)) - { - Delaunay_vertex v = insert_delaunay_landmark_with_copies(W, *cp_it, landmarks_ind, t, landmark_count, torus); - { - // Simple check if the new cells don't have centers too close one to another - std::vector<Full_cell_handle> inc_cells; - std::back_insert_iterator<std::vector<Full_cell_handle>> out(inc_cells); - t.tds().incident_full_cells(v, out); - - std::vector<Sphere_d> spheres; - for (auto i_it = inc_cells.begin(); i_it != inc_cells.end(); ++i_it) - { - std::vector<Point_d> vertices; - for (auto v_it = (*i_it)->vertices_begin(); v_it != (*i_it)->vertices_end(); ++v_it) - vertices.push_back((*v_it)->point()); - spheres.push_back(Sphere_d(vertices.begin(), vertices.end())); - } - for (auto s_it = spheres.begin(); s_it != spheres.end(); ++s_it) - for (auto t_it = s_it+1; t_it != spheres.end(); ++t_it) - { - FT ddc2 = ed.transformed_distance(s_it->center(),t_it->center()); - if (ddc2 < gamma0*gamma0*s_it->squared_radius() || - ddc2 < gamma0*gamma0*t_it->squared_radius()) - { refused_centers2++; } - } - } - - //std::cout << *cp_it << ",\n"; - //make_heap(t, R_max_heap); - point_taken[*cp_it] = true; - rh = sampling_radius(t); - if (experiment1) eps_vector.push_back(pow(1/rh.value,D)); - //std::cout << "rhvalue = " << rh.value << "\n"; - //std::cout << "D = " << - candidate_points.clear(); - torus_search(treeW, D, - rh.center, - alpha*rh.value, - candidate_points); - cp_it = candidate_points.begin(); - /* - // PIECE OF CODE FOR DEBUGGING PURPOSES - - Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - if (triangulation_is_protected(t, delta)) - { - index_list.erase(list_it); - list_it = index_list.begin(); - } - else - { //THAT'S WHERE SOMETHING'S WRONG - t.remove(inserted_v); - landmarks_ind.pop_back(); - landmark_count--; - write_delaunay_mesh(t, W[*list_it], is2d); - is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore - } - */ - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - cp_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //if (list_it != index_list.end()) - // write_delaunay_mesh(t, W[*list_it], is2d); - } - - if (experiment2) epsratio_vector.push_back(rh.value/epsilon0); - if (experiment2) epsslope_vector.push_back( (pow(1/rh.value,D)-pow(1/epsilon0,D))/(landmarks_ind.size() - 48) ); - std::cout << "The iteration ended when cp_count = " << candidate_points.size() << "\n"; - std::cout << "alphaRmax = " << alpha*rh.value << "\n"; - std::cout << "epsilon' = " << rh.value << "\n"; - std::cout << "nbL = " << landmarks_ind.size() << "\n"; - print_statistics(); - //print_vector(landmarks_ind); std::cout << std::endl; - //std::sort(landmarks_ind.begin(), landmarks_ind.end()); - print_vector(landmarks_ind); std::cout << std::endl; - if (experiment3) thetamin_vector[thetamin_vector.size()-1] = sampling_fatness(t); - std::cout << "theta = " << sampling_fatness(t) << "\n"; - //fill_landmarks(W, landmarks, landmarks_ind, torus); - //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, W[0], true); - //std::cout << t << std::endl; -} - -void run_experiment5(Point_Vector& W, - int D, - FT alpha, - FT epsilon, - FT delta0, - FT theta0, - FT gamma0, - //std::vector<std::vector<int>>& full_cells, - bool torus, - bool power_protection - ) -{ - // INITIALIZATION - Delaunay_triangulation t(D); - std::vector<int> landmarks_ind; - int landmark_count = 0; - initialize_statistics(); - if (D == 2) - { - int xw = 6, yw = 4; - // Triangular lattice close to regular triangles h=0.866a ~ 0.875a : 48p - for (int i = 0; i < xw; ++i) - for (int j = 0; j < yw; ++j) - { - Point_d cite1(std::vector<FT>{2.0/xw*i, 2.0/yw*j}); - W.push_back(cite1); // debug purpose - insert_delaunay_landmark_with_copies(W, W.size()-1, - landmarks_ind, t, landmark_count, true); - - Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 2.0/yw*(j+0.5)}); - W.push_back(cite2); // debug purpose - insert_delaunay_landmark_with_copies(W, W.size()-1, - landmarks_ind, t, landmark_count, true); - } - } - else if (D == 3) - { - int wd = 3; - // Body-centered cubic lattice : 54p - for (int i = 0; i < wd; ++i) - for (int j = 0; j < wd; ++j) - for (int k = 0; k < wd; ++k) - { - Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k}); - W.push_back(cite1); // debug purpose - insert_delaunay_landmark_with_copies(W, W.size()-1, - landmarks_ind, t, landmark_count, true); - - Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)}); - W.push_back(cite2); // debug purpose - insert_delaunay_landmark_with_copies(W, W.size()-1, - landmarks_ind, t, landmark_count, true); - } - } - - // ITERATIONS - R_max_handle rh = sampling_radius(t); - Point_d rp = *(Random_point_iterator(D, alpha*rh.value)); - int death_count = 0; - std::cout << "death count " << death_count << " rp = " << rp << "\n"; - while (death_count < 100) - { - std::vector<FT> coords; - for (auto c_it = rh.center.cartesian_begin(), - r_it = rp.cartesian_begin(); - c_it != rh.center.cartesian_end(); - ++c_it, ++r_it) - coords.push_back(*c_it + *r_it); - Point_d new_p(coords); - if (!is_violating_protection(new_p, t, D, delta0, power_protection, theta0, gamma0)) - { - W.push_back(new_p); - insert_delaunay_landmark_with_copies(W, W.size()-1, landmarks_ind, t, landmark_count, torus); - rh = sampling_radius(t); - rp = *(Random_point_iterator(D, alpha*rh.value)); - death_count = 0; - std::cout << "death count " << death_count << " rp = " << rp << "\n"; - } - else - { - rp = *(Random_point_iterator(D, alpha*rh.value)); - death_count++; - std::cout << "death count " << death_count << " rp = " << rp << "\n"; - } - //Point_d new_p = (*rp++) + Vector_d; - } -} - -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Series of experiments -/////////////////////////////////////////////////////////////////////////////////////////////////////////// - -void start_experiments(Point_Vector& W, FT alpha, std::vector<int>& landmarks_ind, FT epsilon) -{ - int experiment_no = 1; - FT delta0 = 0.1; - FT theta0 = 0.1; - FT gamma0 = 0.01; - std::string suffix; - //std::cout << "ようこそジプシー我が神秘の部屋へ:\n"; - while (experiment_no != 0) - { - std::cout << "Enter experiment no (0 to exit): "; - std::cin >> experiment_no; - switch (experiment_no) - { - case 1: - // Experiment 1 - experiment1 = true; - eps_vector = {}; - std::cout << "Enter delta0: "; std::cin >> delta0; - std::cout << "Enter theta0: "; std::cin >> theta0; - std::cout << "Enter gamma0: "; std::cin >> gamma0; - protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true); - write_tikz_plot(eps_vector,"epstime.tikz"); - experiment1 = false; - break; - - case 2: - // Experiment 2 - suffix = ""; - experiment2 = true; - epsratio_vector = {0}; - epsslope_vector = {0}; - std::cout << "File name suffix: "; - std::cin >> suffix; - for (FT alpha = 0.01; alpha < 0.999; alpha += 0.01) - { - landmarks_ind.clear(); - std::cout << "Test for alpha = " << alpha << "\n"; - protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true); - } - write_tikz_plot(epsratio_vector,"epsratio_alpha." + suffix + ".tex"); - write_tikz_plot(epsslope_vector,"epsslope_alpha." + suffix + ".tex"); - experiment2 = false; - break; - - case 3: - // Experiment 3 - experiment3 = true; - thetamin_vector = {}; - gammamin_vector = {}; - theta0 = 0; - gamma0 = 0; - for (FT delta0 = 0; delta0 < 0.999; delta0 += 0.05) - { - landmarks_ind.clear(); - thetamin_vector.push_back(1.0); //0.7489 fatness of the initialization - gammamin_vector.push_back(10); - std::cout << "Test for delta0 = " << delta0 << "\n"; - protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true); - } - write_tikz_plot(thetamin_vector,"thetamin_delta.tex"); - write_tikz_plot(gammamin_vector,"gammamin_delta.tex"); - experiment3 = false; - break; - - // case 4: - // // Experiment 4 - // { - // int dim; - // std::cout << "Enter dimension: "; - // std::cin >> dim; - // Delaunay_triangulation t(dim); - // // for (FT eps = 0.7; eps < 1.1; eps += 0.1) - // // { - // // generate_epsilon_sample_torus(W, eps, dim, t); - // // for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - // // { - // // if (t.is_infinite(v_it)) - // // continue; - // // bool in_cube = true; - // // for (auto xi = v_it->cartesian_begin(); xi != v_it->cartesian_end(); ++xi) - // // if (*xi > 1.0 || *xi < -1.0) - // // { - // // in_cube = false; break; - // // } - // // if (!in_cube) - // // continue; - // // for (auto t.tds().incident_full_cells()) - // // } - // // std::cout << "eps = " << eps << ", real epsilon = " << sampling_radius(t).value << "\n"; - // // } - // // } - // break; - - - case 5: - // Experiment 5 - experiment5 = true; - // std::cout << "Enter dimension: "; - // std::cin >> dim; - - landmarks_ind.clear(); - W.clear(); - run_experiment5(W, alpha, epsilon, delta0, theta0, gamma0, true, true); - experiment5 = false; - break; - } - - } - -} - -#endif diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp deleted file mode 100644 index e448c55d..00000000 --- a/src/Witness_complex/example/witness_complex_cube.cpp +++ /dev/null @@ -1,590 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -// Avoiding the max arity issue with CGAL -#ifndef BOOST_PARAMETER_MAX_ARITY -# define BOOST_PARAMETER_MAX_ARITY 12 -#endif - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <iterator> -#include <chrono> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -#include "Torus_distance.h" -#include "generators.h" -#include "output.h" -//#include "protected_sets/protected_sets.h" -#include "protected_sets/protected_sets_paper2.h" - -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#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> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -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< - FT, Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d> Traits_base; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -//typedef K::Equal_d Equal_d; -//typedef CGAL::Random_points_in_cube_d<CGAL::Point_d<CGAL::Cartesian_d<FT> > > Random_cube_iterator; - -typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation; -typedef Delaunay_triangulation::Facet Facet; -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; - -/*////////////////////////////////////// - * GLOBAL VARIABLES ******************** - *////////////////////////////////////// - -//NA bool toric=false; -bool power_protection = true; -bool grid_points = true; -bool is2d = true; -//FT _sfty = pow(10,-14); -bool torus = false; - - -bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) -{ - std::cout << "Start protection verification\n"; - Euclidean_distance ed; - // Fill the map Vertices -> Numbers - std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex; - int ind = 0; - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - { - if (t.is_infinite(v_it)) - continue; - index_of_vertex[v_it] = ind++; - } - for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) - if (!t.is_infinite(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( vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(0)->point())); - 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 - if (!vertex_is_in_full_cell(v_it, fc_it)) - { - FT dist2 = ed.transformed_distance(center_cs, v_it->point()); - //if the new point is inside the protection ball of a non conflicting simplex - //std::cout << "Dist^2 = " << dist2 << " (r+delta)*(r+delta) = " << (r+delta)*(r+delta) << " r^2 = " << r*r <<"\n"; - if (!power_protection) - if (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r) - { - write_delaunay_mesh(t, v_it->point(), is2d); - // Output the problems - std::cout << "Problematic vertex " << index_of_vertex[v_it] << " "; - std::cout << "Problematic cell "; - for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) - if (!t.is_infinite(*vh_it)) - std::cout << index_of_vertex[*vh_it] << " "; - std::cout << "\n"; - std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", (r+delta)^2 = " << (r+delta)*(r+delta) << "\n"; - return false; - } - if (power_protection) - if (dist2 <= r*r+delta*delta && dist2 >= r*r) - { - write_delaunay_mesh(t, v_it->point(), is2d); - 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; -} - -////////////////////////////////////////////////////////////////////////////////////////////////////////// -// SAMPLING RADIUS -////////////////////////////////////////////////////////////////////////////////////////////////////////// - -FT sampling_radius(Delaunay_triangulation& t, FT epsilon0) -{ - FT epsilon2 = 0; - Point_d control_point; - 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; - control_point = (*vertices.begin()); - } - } - if (epsilon2 < epsilon0*epsilon0) - { - std::cout << "ACHTUNG! E' < E\n"; - std::cout << "eps = " << epsilon0 << " eps' = " << sqrt(epsilon2) << "\n"; - write_delaunay_mesh(t, control_point, is2d); - } - return sqrt(epsilon2); -} - -FT point_sampling_radius_by_delaunay(Point_Vector& points, FT epsilon0) -{ - Delaunay_triangulation t(points[0].size()); - t.insert(points.begin(), points.end()); - return sampling_radius(t, epsilon0); -} - -// A little script to make a tikz histogram of epsilon distribution -// Returns the average epsilon -FT epsilon_histogram(Delaunay_triangulation& t, int n) -{ - FT epsilon_max = 0; //sampling_radius(t,0); - FT sum_epsilon = 0; - int count_simplices = 0; - std::vector<int> histo(n+1, 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 r = sqrt(Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin()))); - if (r > epsilon_max) - epsilon_max = r; - sum_epsilon += r; - count_simplices++; - histo[floor(r/epsilon_max*n)]++; - } - std::ofstream ofs ("histogram.tikz", std::ofstream::out); - FT barwidth = 20.0/n; - int max_value = *(std::max_element(histo.begin(), histo.end())); - std::cout << max_value << std::endl; - FT ten_power = pow(10, ceil(log10(max_value))); - FT max_histo = ten_power; - if (max_value/ten_power < 2) - max_histo = 0.2*ten_power; - if (max_value/ten_power < 5) - max_histo = 0.5*ten_power; - std::cout << ceil(log10(max_value)) << std::endl << max_histo << std::endl; - FT unitht = max_histo/10.0; - - ofs << "\\draw[->] (0,0) -- (0,11);\n" << - "\\draw[->] (0,0) -- (21,0);\n" << - "\\foreach \\i in {1,...,10}\n" << - "\\draw (0,\\i) -- (-0.1,\\i);\n" << - "\\foreach \\i in {1,...,20}\n" << - "\\draw (\\i,0) -- (\\i,-0.1);\n" << - - "\\node at (-1,11) {$\\epsilon$};\n" << - "\\node at (22,-1) {$\\epsilon/\\epsilon_{max}$};\n" << - "\\node at (-0.5,-0.5) {0};\n" << - "\\node at (-0.5,10) {" << max_histo << "};\n" << - "\\node at (20,-0.5) {1};\n"; - - - for (int i = 0; i < n; ++i) - ofs << "\\draw (" << barwidth*i << "," << histo[i]/unitht << ") -- (" - << barwidth*(i+1) << "," << histo[i]/unitht << ") -- (" - << barwidth*(i+1) << ",0) -- (" << barwidth*i << ",0) -- cycle;\n"; - - ofs.close(); - - //return sum_epsilon/count_simplices; - return epsilon_max; -} - -FT epsilon_histogram_by_delaunay(Point_Vector& points, int n) -{ - Delaunay_triangulation t(points[0].size()); - t.insert(points.begin(), points.end()); - return epsilon_histogram(t, n); -} - - -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(); - std::vector<FT> orig_vector; - for (int i=0; i<D; i++) - orig_vector.push_back(0); - Point_d origin(orig_vector); - - //******************** Constructing a WL matrix - int nbP = W.size(); - Euclidean_distance ed; - FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); - std::vector<Point_d> 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<double> 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 <int> > WL(nbP); - - //********************** Neighbor search in a Kd tree - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(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); - - //******************** 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"; - 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<int> curr_perturb; - L.search(std::insert_iterator<std::vector<int>>(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<FT> 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; - */ - return 0; -} - -int main (int argc, char * const argv[]) -{ - power_protection = true;//false; - grid_points = false;//true; - torus = true; - - if (argc != 4) - { - std::cerr << "Usage: " << argv[0] - << " nbP dim delta\n"; - return 0; - } - int nbP = atoi(argv[1]); - int dim = atoi(argv[2]); - double theta0 = atof(argv[3]); - //double delta = atof(argv[3]); - - is2d = (dim == 2); - - std::cout << "Let the carnage begin!\n"; - Point_Vector point_vector; - if (grid_points) - { - generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim, torus); - nbP = (int)pow((int)pow(nbP, 1.0/dim), dim); - } - else - generate_points_random_box(point_vector, nbP, dim); - FT epsilon = point_sampling_radius_by_delaunay(point_vector, 0); - //FT epsilon = epsilon_histogram_by_delaunay(point_vector,50); - 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; - - int n = 1; - std::vector<FT> values(n,0); - std::vector<FT> time(n,0); - - //FT step = 0.001; - //FT delta = 0.01*epsilon; - //FT alpha = 0.5; - //FT step = atof(argv[3]); - - start_experiments(point_vector, theta0, chosen_landmarks, epsilon); - - // for (int i = 0; i < n; i++) - // //for (int i = 0; bl > 0; i++) - // { - // //std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n"; - // //double delta = pow(10, -(1.0*i)/2); - // //delta = step*i*epsilon; - // //theta0 = step*i; - // std::cout << "delta/epsilon = " << delta/epsilon << std::endl; - // std::cout << "theta0 = " << theta0 << std::endl; - // // Averaging the result - // int sum_values = 0; - // int nb_iterations = 1; - // std::vector<std::vector<int>> full_cells; - // for (int i = 0; i < nb_iterations; ++i) - // { - // //L = {}; - // chosen_landmarks = {}; - // //full_cells = {}; - // //timer.start(); - // //protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, epsilon, alpha, theta0, full_cells, torus, power_protection); - // protected_delaunay(point_vector, chosen_landmarks, delta, epsilon, alpha, theta0, torus, power_protection); - // //timer.stop(); - // sum_values += chosen_landmarks.size(); - // } - // //FT epsilon2 = point_sampling_radius_by_delaunay(L, epsilon); - // //std::cout << "Final epsilon = " << epsilon2 << ". Ratio = " << epsilon2/epsilon << std::endl; - // //write_points("landmarks/initial_landmarks",L); - // //std::cout << "delta/epsilon' = " << delta/epsilon2 << std::endl; - // FT nbL = (sum_values*1.0)/nb_iterations; - // //values[i] = pow((1.0*nbL)/nbP, -1.0/dim); - // values[i] = (1.0*nbL)/nbP; - // std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl; - // //landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); - // time[i] = timer.time(); - // timer.reset(); - // //write_points("landmarks/landmarks0",L); - // } - - // // OUTPUT A PLOT - // FT hstep = 20.0/(n-1); - // FT wstep = 10.0; - - // std::ofstream ofs("N'Nplot.tikz", std::ofstream::out); - // ofs << "\\draw[red] (0," << wstep*values[0] << ")"; - // for (int i = 1; i < n; ++i) - // ofs << " -- (" << hstep*i << "," << wstep*values[i] << ")"; - // ofs << ";\n"; - // ofs.close(); - /* - wstep = 0.1; - ofs = std::ofstream("time.tikz", std::ofstream::out); - ofs << "\\draw[red] (0," << wstep*time[0] << ")"; - for (int i = 1; i < n; ++i) - ofs << " -- (" << hstep*i << "," << wstep*time[i] << ")"; - ofs << ";\n"; - ofs.close(); - - - 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(); - */ -} diff --git a/src/Witness_complex/example/witness_complex_cubic_systems.cpp b/src/Witness_complex/example/witness_complex_cubic_systems.cpp deleted file mode 100644 index 2f4ee1cb..00000000 --- a/src/Witness_complex/example/witness_complex_cubic_systems.cpp +++ /dev/null @@ -1,547 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <iterator> - -#include <sys/types.h> -#include <sys/stat.h> -#include <unistd.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -#include "Torus_distance.h" - -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.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/Delaunay_triangulation.h> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -//typedef CGAL::Cartesian_d<double> K; -//typedef CGAL::Point_d<K> 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<Traits_base> Euclidean_distance; - - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -//typedef K::Equal_d Equal_d; -//typedef CGAL::Random_points_in_cube_d<CGAL::Point_d<CGAL::Cartesian_d<FT> > > Random_cube_iterator; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -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; - -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<double> 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 <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - - -void write_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<int>& 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<int>& landmarks_ind, std::vector<bool> & 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<double> 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<int>& landmarks_ind, std::vector<bool>& curr_pattern, std::vector<bool>::iterator curr_pattern_it, std::vector<bool>::iterator bool_it, std::vector<bool>::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<int>& landmarks_ind, std::vector<bool>& 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<bool> 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<int>& landmarks_ind) -{ - //******************** Preface: origin point - int D = W[0].size(); - std::vector<FT> orig_vector; - for (int i=0; i<D; i++) - orig_vector.push_back(0); - Point_d origin(orig_vector); - - //******************** Constructing a WL matrix - int nbP = W.size(); - int nbL = landmarks.size(); - Euclidean_distance ed; - FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); - std::vector<Point_d> 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<double> 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 <int> > WL(nbP); - - //********************** Neighbor search in a Kd tree - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(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<int> curr_perturb; - L.search(std::insert_iterator<std::vector<int>>(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<FT> 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<bool> face_centers(D, false); - int bl = 0; //Bad links - std::vector<std::vector<bool>> 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<int> 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<bool> 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<bool> 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<int> 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 deleted file mode 100644 index 7f8b985f..00000000 --- a/src/Witness_complex/example/witness_complex_epsilon.cpp +++ /dev/null @@ -1,55 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <vector> - -#include <CGAL/Epick_d.h> -#include <CGAL/enum.h> - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -typedef K::FT FT; -typedef K::Hyperplane_d Hyperplane_d; -typedef K::Has_on_positive_side_d Has_on_positive_side_d; - -int main () -{ - std::vector<Point_d> vertices; - Point_d v1(std::vector<FT>({-1,1})); - Point_d v2(std::vector<FT>({1,-1})); - vertices.push_back(v1); - vertices.push_back(v2); - Point_d p(std::vector<FT>({-1,-1})); - Hyperplane_d hp(vertices.begin(), vertices.end()); - //Hyperplane_d hp(vertices.begin(), vertices.end(), p, CGAL::ON_POSITIVE_SIDE); - if (Has_on_positive_side_d()(hp, p)) - std::cout << "OK\n"; - else - std::cout << "NOK\n"; - CGAL::Oriented_side side_p = K::Oriented_side_d()(hp, p); - if (side_p == CGAL::ZERO) - std::cout << "Point (-1,-1) is on the line passing through (-1,1) and (1,-1)"; - CGAL::Oriented_side side_v2 = K::Oriented_side_d()(hp, v2); - if (side_v2 != CGAL::ZERO) - std::cout << "Point (1,-1) is not on the line passing through (-1,1) and (1,-1)"; -} diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp deleted file mode 100644 index 49383154..00000000 --- a/src/Witness_complex/example/witness_complex_flat_torus.cpp +++ /dev/null @@ -1,851 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <iterator> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -//#include <boost/filesystem.hpp> - -//#include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> - -#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/point_generators_d.h> -#include <CGAL/constructions_d.h> -#include <CGAL/Fuzzy_sphere.h> -#include <CGAL/Random.h> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -//typedef CGAL::Cartesian_d<double> K; -//typedef CGAL::Point_d<K> 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<Traits_base> Euclidean_distance; - -/** - * \brief Class of distance in a flat torus in dimension 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<FT,D>& 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<FT,D>& r, - std::vector<FT>& 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<FT,D>& 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<FT,D>& r, - std::vector<FT>& 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); - } - -}; -*/ - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -//typedef K::Equal_d Equal_d; -//typedef CGAL::Random_points_in_cube_d<CGAL::Point_d<CGAL::Cartesian_d<FT> > > Random_cube_iterator; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; - -bool toric=false; - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , Point_Vector & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - Point_d p(point.begin(), point.end()); - if (point.size() != 1) - points.push_back(p); - } - in_file.close(); -} - -void generate_points_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<double> point; - int cell_i = i; - for (int l = 0; l < D; ++l) - { - point.push_back((2.0/width)*(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<double> 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++); - } -} - -/* NOT TORUS RELATED - */ -void generate_points_sphere(Point_Vector& W, int nbP, int dim) -{ - CGAL::Random_points_on_sphere_d<Point_d> rp(dim,1); - for (int i = 0; i < nbP; i++) - W.push_back(*rp++); -} -/* -void 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<double> 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 <int> > & 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<Point_d> convert_to_torus(std::vector< Point_d>& points) -{ - std::vector< Point_d > points_torus; - for (auto p: points) - { - FT theta = M_PI*p[0]; - FT phi = M_PI*p[1]; - std::vector<FT> p_torus; - p_torus.push_back((1+0.2*cos(theta))*cos(phi)); - p_torus.push_back((1+0.2*cos(theta))*sin(phi)); - p_torus.push_back(0.2*sin(theta)); - points_torus.push_back(Point_d(p_torus)); - } - return points_torus; -} - -void write_points_torus( std::string file_name, std::vector< Point_d > & points) -{ - std::ofstream ofs (file_name, std::ofstream::out); - std::vector<Point_d> points_torus = convert_to_torus(points); - for (auto w : points_torus) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); -} - -void write_points( std::string file_name, std::vector< Point_d > & points) -{ - if (toric) write_points_torus(file_name, points); - else - { - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : points) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); - } -} - - -void write_edges_torus(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - Point_Vector l_torus = convert_to_torus(landmarks); - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = l_torus[u].cartesian_begin(); it != l_torus[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = l_torus[v].cartesian_begin(); it != l_torus[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); -} - -void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - if (toric) write_edges_torus(file_name, witness_complex, landmarks); - else - { - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); - } -} - - -/** Function that chooses landmarks from W and place it in the kd-tree L. - * Note: nbL hould be removed if the code moves to Witness_complex - */ -void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - 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::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; - } - } - -/** \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<int>& 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<double> point; - std::vector<double> 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<int>& landmarks_ind) -{ - //********************Preface: origin point - int D = W[0].size(); - std::vector<FT> orig_vector; - for (int 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(); - //Point_Vector landmarks_ = landmarks; - Euclidean_distance ed; - //Equal_d ed; - //Point_d p1(std::vector<FT>({0.8,0.8})), p2(std::vector<FT>({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<Point_d> 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<double> 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 <int> > WL(nbP); - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(nb_cells*nbL), - typename Tree::Splitter(), - traits); - /*Tree2 L2(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])) ); - CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks_ext[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"; - 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; - } - 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); - /* - 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; - 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, count_bad, count_good)) - { - //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), 0, traits); - std::vector<int> curr_perturb; - L.search(std::insert_iterator<std::vector<int>>(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<int>(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 - - for (auto u: perturbL) - { - Random_point_iterator rp(D,sqrt(lambda)/8); - //std::cout << landmarks[u] << std::endl; - - std::vector<FT> 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(); - } - /* - 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<double>({0.1,0.1})), Point_d(std::vector<double>({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_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) - { - 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<std::vector< int > > WL(nbP); - //std::set<int> L; - Point_Vector L; - std::vector<int> 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); - - //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, 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(); - */ -} diff --git a/src/Witness_complex/example/witness_complex_from_off.cpp b/src/Witness_complex/example/witness_complex_from_off.cpp deleted file mode 100644 index 948f09a8..00000000 --- a/src/Witness_complex/example/witness_complex_from_off.cpp +++ /dev/null @@ -1,184 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <sys/types.h> -#include <sys/stat.h> - -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" - -using namespace Gudhi; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::vector< std::vector <double> > Point_Vector; -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , std::vector< std::vector< double > > & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -/** - * \brief Rock age method of reading off file - * - */ -inline void -off_reader_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; - // 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); - } - in_file.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 - Witness_complex<> witnessComplex; - - /* - std::cout << "Let the carnage begin!\n"; - start = clock(); - Point_Vector point_vector; - off_reader_cust(file_name, point_vector); - std::cout << "Successfully read the points\n"; - witnessComplex.setNbL(nbL); - witnessComplex.witness_complex_from_points(point_vector); - end = clock(); - std::cout << "Howdy world! The process took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - char buffer[100]; - int i = sprintf(buffer,"%s_%s_result.txt",argv[1],argv[2]); - if (i >= 0) - { - std::string out_file = (std::string)buffer; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - } - */ - std::cout << "Let the carnage begin!\n"; - Point_Vector point_vector; - off_reader_cust(file_name, point_vector); - //std::cout << "Successfully read the points\n"; - witnessComplex.setNbL(nbL); - // witnessComplex.witness_complex_from_points(point_vector); - std::vector<std::vector< int > > WL; - std::set<int> L; - start = clock(); - //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); - witnessComplex.landmark_choice_by_random_points(point_vector, point_vector.size(), L); - witnessComplex.nearest_landmarks(point_vector,L,WL); - end = clock(); - std::cout << "Landmark choice took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - // Write the WL matrix in a file - mkdir("output", S_IRWXU); - const size_t last_slash_idx = file_name.find_last_of("/"); - if (std::string::npos != last_slash_idx) - { - file_name.erase(0, last_slash_idx + 1); - } - std::string out_file = "output/"+file_name+"_"+argv[2]+".wl"; - //write_wl(out_file,WL); - start = clock(); - witnessComplex.witness_complex(WL); - // - end = clock(); - std::cout << "Howdy world! The process took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - - out_file = "output/"+file_name+"_"+argv[2]+".stree"; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - - out_file = "output/"+file_name+"_"+argv[2]+".badlinks"; - std::ofstream ofs2(out_file, std::ofstream::out); - witnessComplex.write_bad_links(ofs2); - ofs2.close(); -} diff --git a/src/Witness_complex/example/witness_complex_from_wl_matrix.cpp b/src/Witness_complex/example/witness_complex_from_wl_matrix.cpp deleted file mode 100644 index 614bb945..00000000 --- a/src/Witness_complex/example/witness_complex_from_wl_matrix.cpp +++ /dev/null @@ -1,148 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -//#include <boost/filesystem.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::vector< std::vector <double> > Point_Vector; -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust ( std::string file_name , std::vector< std::vector< double > > & points) -{ - std::ifstream in_file (file_name.c_str(),std::ios::in); - if(!in_file.is_open()) - { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while( getline ( in_file , line ) ) - { - std::vector< double > point; - std::istringstream iss( line ); - while(iss >> x) { point.push_back(x); } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -void write_wl( std::string file_name, std::vector< std::vector <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - -void read_wl( std::string file_name, std::vector< std::vector <int> > & WL) -{ - 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; - int x; - while( getline ( in_file , line ) ) - { - std::vector< int > witness; - std::istringstream iss( line ); - while(iss >> x) { witness.push_back(x); } - WL.push_back(witness); - } - in_file.close(); - -} - -int main (int argc, char * const argv[]) -{ - if (argc != 2) - { - std::cerr << "Usage: " << argv[0] - << " path_to_point_file \n"; - return 0; - } - /* - boost::filesystem::path p; - - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - std::string file_name = argv[1]; - //int nbL = atoi(argv[2]); - - clock_t start, end; - //Construct the Simplex Tree - Witness_complex<> witnessComplex; - - std::cout << "Let the carnage begin!\n"; - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - //std::cout << "Successfully read the points\n"; - // witnessComplex.witness_complex_from_points(point_vector); - std::vector<std::vector< int > > WL; - read_wl(file_name,WL); - witnessComplex.setNbL(WL[0].size()); - // Write the WL matrix in a file - std::string out_file; - write_wl(out_file,WL); - start = clock(); - witnessComplex.witness_complex(WL); - // - end = clock(); - std::cout << "Howdy world! The process took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - out_file = "output/"+file_name+"_"+argv[2]+".stree"; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - -} diff --git a/src/Witness_complex/example/witness_complex_knn_landmarks.cpp b/src/Witness_complex/example/witness_complex_knn_landmarks.cpp deleted file mode 100644 index c45bc0c1..00000000 --- a/src/Witness_complex/example/witness_complex_knn_landmarks.cpp +++ /dev/null @@ -1,210 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -#include "generators.h" -#include "output.h" -//#include <boost/filesystem.hpp> - -//#include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef CGAL::Search_traits< - FT, Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d> Traits_base; -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -typedef CGAL::Orthogonal_k_neighbor_search<STraits> 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<int, int> Point_etiquette_map; - -typedef std::vector<Point_d> Point_Vector; - -/** 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_to_tree(Point_Vector &W, int nbP, Point_etiquette_map &L_i, int nbL, std::vector< std::vector <int> > &WL) -{ - std::cout << "Enter landmark choice to kd tree\n"; - std::vector<Point_d> landmarks; - int chosen_landmark; - //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false); - Point_d* p; - srand(24660); - for (int i = 0; i < nbL; i++) - { - // while (!res.second) - // { - chosen_landmark = rand()%nbP; - p = &W[chosen_landmark]; - //L_i.emplace(chosen_landmark,i); - // } - landmarks.push_back(*p); - //std::cout << "Added landmark " << chosen_landmark << std::endl; - } - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(nbL), - typename Tree::Splitter(), - STraits((Point_d*)&(landmarks[0]))); - /*} - - -void d_nearest_landmarks(Point_Vector &W, Tree &L, Point_etiquette_map &L_i, std::vector< std::vector <int> > &WL) -{*/ - std::cout << "Enter (D+1) nearest landmarks\n"; - std::cout << "Size of the tree is " << L.size() << std::endl; -//int nbP = W.size(); - int D = W[0].size(); - 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 - K_neighbor_search search(L, w, D+1, FT(0), true, - CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,CGAL::Euclidean_distance<Traits_base>>((Point_d*)&(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 << i << " " << it->first << ": " << it->second << std::endl; - } - } -} - -int main (int argc, char * const argv[]) -{ - if (argc != 3) - { - std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL \n"; - return 0; - } - /* - boost::filesystem::path p; - - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - std::string file_name = argv[1]; - int nbL = atoi(argv[2]); - - clock_t start, end; - //Construct the Simplex Tree - Witness_complex<> witnessComplex; - - std::cout << "Let the carnage begin!\n"; - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - //std::cout << "Successfully read the points\n"; - witnessComplex.setNbL(nbL); - // witnessComplex.witness_complex_from_points(point_vector); - int nbP = point_vector.size(); - std::vector<std::vector< int > > WL(nbP); - //std::set<int> L; - Tree L; - Point_etiquette_map L_i; - start = clock(); - //witnessComplex.landmark_choice_by_furthest_points(point_vector, point_vector.size(), WL); - landmark_choice_to_tree(point_vector, nbP, L_i, nbL, WL); - //d_nearest_landmarks(point_vector, L, L_i, WL); - end = clock(); - std::cout << "Landmark choice took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - // Write the WL matrix in a file - mkdir("output", S_IRWXU); - const size_t last_slash_idx = file_name.find_last_of("/"); - if (std::string::npos != last_slash_idx) - { - file_name.erase(0, last_slash_idx + 1); - } - std::string out_file = "output/"+file_name+"_"+argv[2]+".wl"; - write_wl(out_file,WL); - start = clock(); - witnessComplex.witness_complex(WL); - // - end = clock(); - std::cout << "Howdy world! The process took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - /* - char buffer[100]; - int i = sprintf(buffer,"%s_%s_result.txt",argv[1],argv[2]); - if (i >= 0) - { - std::string out_file = (std::string)buffer; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - } - */ - - out_file = "output/"+file_name+"_"+argv[2]+".stree"; - std::ofstream ofs (out_file, std::ofstream::out); - witnessComplex.st_to_file(ofs); - ofs.close(); - - out_file = "output/"+file_name+"_"+argv[2]+".badlinks"; - std::ofstream ofs2(out_file, std::ofstream::out); - //witnessComplex.write_bad_links(ofs2); - ofs2.close(); -} diff --git a/src/Witness_complex/example/witness_complex_perturbations.cpp b/src/Witness_complex/example/witness_complex_perturbations.cpp deleted file mode 100644 index f78bcdab..00000000 --- a/src/Witness_complex/example/witness_complex_perturbations.cpp +++ /dev/null @@ -1,462 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <set> -#include <iterator> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -//#include <boost/filesystem.hpp> - -//#include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> - -#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/point_generators_d.h> -#include <CGAL/constructions_d.h> -#include <CGAL/Fuzzy_sphere.h> -#include <CGAL/Origin.h> - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef CGAL::Search_traits< - FT, Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d> Traits_base; -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -typedef CGAL::Orthogonal_k_neighbor_search<STraits> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; -//typedef K::Equal_d Equal_d; -typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; -/** - * \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 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<double> 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 <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - -void write_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<int>& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector<Point_d> landmarks; - int chosen_landmark; - //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false); - Point_d* p; - srand(24660); - for (int i = 0; i < nbL; i++) - { - // while (!res.second) - // { - chosen_landmark = rand()%nbP; - 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 landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector<Point_d> landmarks; - int chosen_landmark = 0; - //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false); - Point_d* p; - CGAL::Random rand; - for (int i = 0; i < nbL; i++) - { - // while (!res.second) - // { - do chosen_landmark = rand.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); - // } - 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<int>& landmarks_ind) -{ - //******************** Constructing a WL matrix - int nbP = W.size(); - int nbL = landmarks.size(); - //Point_Vector landmarks_ = landmarks; - Euclidean_distance ed; - //Equal_d ed; - FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); - //FT lambda = 0.1;//Euclidean_distance(); - std::vector< std::vector <int> > WL(nbP); - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(nbL), - typename Tree::Splitter(), - STraits(&(landmarks[0]))); - /*Tree2 L2(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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; - 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; - Point_d& w = W[i]; - //std::cout << "Safely constructed a point\n"; - ////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<std::ptrdiff_t,Point_d*,CGAL::Euclidean_distance<Traits_base>>(&(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 << "\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); - - //******************** 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; - int count_badlinks = 0; - std::vector< int > count_bad(D); - std::vector< int > count_good(D); - //std::cout << "Bad links around "; - for (auto u: witnessComplex.complex_vertex_range()) - if (!witnessComplex.has_good_link(u, count_bad, count_good)) - { - //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)*2, 0, STraits(&(landmarks[0]))); - L.search(std::insert_iterator<std::set<int>>(perturbL,perturbL.begin()),fs); - //L.search(std::inserter(perturbL,perturbL.begin()),fs); - //L.search(std::ostream_iterator<int>(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 << "Bad 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)/4); - //std::cout << landmarks[u] << std::endl; - - std::vector<FT> point; - for (int i = 0; i < D; i++) - { - point.push_back(W[landmarks_ind[u]][i] + (*rp)[i]); - } - 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; - // Write the WL matrix in a file - - /* - 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(); - } - */ - //witnessComplex.write_badlinks("badlinks"); - //write_edges_gnuplot("landmarks/edges", witnessComplex, landmarks); - return count_badlinks; -} - - -int main (int argc, char * const argv[]) -{ - if (argc != 3) - { - std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL \n"; - return 0; - } - /* - boost::filesystem::path p; - - for (; argc > 2; --argc, ++argv) - p /= argv[1]; - */ - std::string file_name = argv[1]; - int nbL = atoi(argv[2]); - - //clock_t start, end; - //Construct the Simplex Tree - //Witness_complex<> witnessComplex; - - std::cout << "Let the carnage begin!\n"; - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - //std::cout << "Successfully read the points\n"; - //witnessComplex.setNbL(nbL); - // witnessComplex.witness_complex_from_points(point_vector); - int nbP = point_vector.size(); - //std::vector<std::vector< int > > WL(nbP); - //std::set<int> L; - Point_Vector L; - std::vector<int> chosen_landmarks; - //Point_etiquette_map L_i; - //start = clock(); - //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++) - for (int i = 0; i < 1; 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(); - - /* - 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_protected_delaunay.cpp b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp deleted file mode 100644 index 77a167a5..00000000 --- a/src/Witness_complex/example/witness_complex_protected_delaunay.cpp +++ /dev/null @@ -1,604 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <iterator> -#include <chrono> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -#include "Torus_distance.h" - -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> -#include <CGAL/Kernel_d/Sphere_d.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/Delaunay_triangulation.h> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::Point_d Point_d; -//typedef CGAL::Cartesian_d<double> K; -//typedef CGAL::Point_d<K> 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<Traits_base> Euclidean_distance; - - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -//typedef K::Equal_d Equal_d; -//typedef CGAL::Random_points_in_cube_d<CGAL::Point_d<CGAL::Cartesian_d<FT> > > Random_cube_iterator; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -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; - -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<double> 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<double> 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 <int> > & WL) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : WL) - { - for (auto l: w) - ofs << l << " "; - ofs << "\n"; - } - ofs.close(); -} - - -void write_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<int>& 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<int>& 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<FT> 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++; -} - - - - -//////////////////////////////////////////////////////////////////////// -// OLD CODE VVVVVVVV -//////////////////////////////////////////////////////////////////////// - - -/* -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<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()); - 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<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()); - 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<int>& 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<int>& 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<int>& 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<int> index_list; - // shuffle the list of indexes (via a vector) - { - std::vector<int> 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<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 (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<int>::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<int>& landmarks_ind) -{ - //******************** Preface: origin point - int D = W[0].size(); - std::vector<FT> orig_vector; - for (int i=0; i<D; i++) - orig_vector.push_back(0); - Point_d origin(orig_vector); - - //******************** Constructing a WL matrix - int nbP = W.size(); - Euclidean_distance ed; - FT lambda = ed.transformed_distance(landmarks[0],landmarks[1]); - std::vector<Point_d> 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<double> 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 <int> > WL(nbP); - - //********************** Neighbor search in a Kd tree - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(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<int> curr_perturb; - L.search(std::insert_iterator<std::vector<int>>(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<FT> 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<int> 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 deleted file mode 100644 index bf3015fa..00000000 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ /dev/null @@ -1,457 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Siargey Kachanovich - * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - -#include <iostream> -#include <fstream> -#include <ctime> -#include <utility> -#include <algorithm> -#include <set> -#include <iterator> - -#include <sys/types.h> -#include <sys/stat.h> -//#include <stdlib.h> - -//#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/Witness_complex.h" -#include "gudhi/reader_utils.h" -#include "generators.h" -#include "output.h" -//#include <boost/filesystem.hpp> - -//#include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Cartesian_d.h> -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/property_map.h> -#include <CGAL/Epick_d.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Kd_tree.h> -#include <CGAL/Euclidean_distance.h> - -#include <CGAL/Kernel_d/Vector_d.h> -#include <CGAL/point_generators_d.h> -#include <CGAL/constructions_d.h> -#include <CGAL/Fuzzy_sphere.h> -#include <CGAL/Random.h> - - -#include <boost/tuple/tuple.hpp> -#include <boost/iterator/zip_iterator.hpp> -#include <boost/iterator/counting_iterator.hpp> -#include <boost/range/iterator_range.hpp> - -using namespace Gudhi; -//using namespace boost::filesystem; - -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef CGAL::Search_traits< - FT, Point_d, - typename K::Cartesian_const_iterator_d, - typename K::Construct_cartesian_const_iterator_d> Traits_base; -typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance; - -typedef std::vector< Vertex_handle > typeVectorVertex; - -//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; - -typedef CGAL::Search_traits_adapter< - std::ptrdiff_t, Point_d*, Traits_base> STraits; -//typedef K TreeTraits; -//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter; -//typedef CGAL::Kd_tree<STraits> Kd_tree; -typedef CGAL::Orthogonal_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> 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<int, int> Point_etiquette_map; -typedef CGAL::Kd_tree<STraits> Tree2; - -typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere; - -typedef std::vector<Point_d> Point_Vector; - -//typedef K::Equal_d Equal_d; - -bool toric=false; - -std::vector<Point_d> convert_to_torus(std::vector< Point_d>& points) -{ - std::vector< Point_d > points_torus; - for (auto p: points) - { - FT theta = M_PI*p[0]; - FT phi = M_PI*p[1]; - std::vector<FT> p_torus; - p_torus.push_back((1+0.2*cos(theta))*cos(phi)); - p_torus.push_back((1+0.2*cos(theta))*sin(phi)); - p_torus.push_back(0.2*sin(theta)); - points_torus.push_back(Point_d(p_torus)); - } - return points_torus; -} - -/** Function that chooses landmarks from W and place it in the kd-tree L. - * Note: nbL hould be removed if the code moves to Witness_complex - */ -void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind) -{ - std::cout << "Enter landmark choice to kd tree\n"; - //std::vector<Point_d> landmarks; - int chosen_landmark; - //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false); - Point_d* p; - CGAL::Random rand; - for (int i = 0; i < nbL; i++) - { - // while (!res.second) - // { - do chosen_landmark = rand.get_int(0,nbP); - while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end()); - //rand++; - //std::cout << "Chose " << chosen_landmark << std::endl; - p = &W[chosen_landmark]; - //L_i.emplace(chosen_landmark,i); - // } - landmarks.push_back(*p); - landmarks_ind.push_back(chosen_landmark); - //std::cout << "Added landmark " << chosen_landmark << std::endl; - } - } - -/** \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<int>& landmarks_ind) -{ - assert(W[0].size() == 4); //4-dimensionality required - FT phi = (1+sqrt(5))/2; - FT phi_1 = FT(1)/phi; - std::vector<FT> 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<FT>({a,b,c,d}))); - // 8 vertices - for (FT a = -0.5; a < 1; a += 1) - { - landmarks.push_back(Point_d(std::vector<FT>({a,0,0,0}))); - landmarks.push_back(Point_d(std::vector<FT>({0,a,0,0}))); - landmarks.push_back(Point_d(std::vector<FT>({0,0,a,0}))); - landmarks.push_back(Point_d(std::vector<FT>({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<FT>({a,b,c,0}))); - landmarks.push_back(Point_d(std::vector<FT>({b,a,0,c}))); - landmarks.push_back(Point_d(std::vector<FT>({c,0,a,b}))); - landmarks.push_back(Point_d(std::vector<FT>({0,c,b,a}))); - landmarks.push_back(Point_d(std::vector<FT>({a,c,0,b}))); - landmarks.push_back(Point_d(std::vector<FT>({a,0,b,c}))); - landmarks.push_back(Point_d(std::vector<FT>({c,b,0,a}))); - landmarks.push_back(Point_d(std::vector<FT>({0,b,a,c}))); - landmarks.push_back(Point_d(std::vector<FT>({b,0,c,a}))); - landmarks.push_back(Point_d(std::vector<FT>({0,a,c,b}))); - landmarks.push_back(Point_d(std::vector<FT>({b,c,a,0}))); - landmarks.push_back(Point_d(std::vector<FT>({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<int>& landmarks_ind) -{ - //********************Preface: origin point - clock_t start, end; - int D = W[0].size(); - std::vector<FT> orig_vector; - for (int 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(); - Euclidean_distance ed; - 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 <int> > WL(nbP); - Tree L(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(nbL), - typename Tree::Splitter(), - traits); - /*Tree2 L2(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(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; - start = clock(); - 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, FT(0), true, - //CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])) ); - CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])) ); - //std::cout << "Safely found nearest landmarks\n"; - for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it) - { - //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; - //Point_etiquette_map::iterator itm = L_i.find(it->first); - //assert(itm != L_i.end()); - //std::cout << "Entered KNN_it with point at distance " << it->second << "\n"; - WL[i].push_back(it->first); - //std::cout << "ITFIRST " << it->first << std::endl; - //std::cout << i << " " << it->first << ": " << it->second << std::endl; - } - 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"; - end = clock(); - std::cout << "Landmark choice for " << nbL << " landmarks took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \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); - start = clock(); - witnessComplex.witness_complex(WL); - // - end = clock(); - std::cout << "Howdy world! The process took " - << (double)(end-start)/CLOCKS_PER_SEC << " s. \n"; - //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; - 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)) - { - //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), 0, traits); - std::vector<int> curr_perturb; - L.search(std::insert_iterator<std::vector<int>>(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<int>(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 - - for (auto u: perturbL) - { - Random_point_iterator rp(D,sqrt(lambda)/8*nbL/count_badlinks); - //std::cout << landmarks[u] << std::endl; - - std::vector<FT> 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<double>({0.1,0.1})), Point_d(std::vector<double>({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<std::vector< int > > WL(nbP); - //std::set<int> L; - Point_Vector L; - std::vector<int> 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); - //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); - //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(); - - /* - 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(); - */ -} |