diff options
Diffstat (limited to 'src/Witness_complex/example/witness_complex_cube.cpp')
| -rw-r--r-- | src/Witness_complex/example/witness_complex_cube.cpp | 808 |
1 files changed, 242 insertions, 566 deletions
diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp index a9a2959b..e448c55d 100644 --- a/src/Witness_complex/example/witness_complex_cube.cpp +++ b/src/Witness_complex/example/witness_complex_cube.cpp @@ -20,6 +20,11 @@ * 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> @@ -37,6 +42,10 @@ #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> @@ -106,8 +115,6 @@ 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; @@ -117,449 +124,84 @@ typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle; typedef K::Sphere_d Sphere_d; typedef K::Hyperplane_d Hyperplane_d; +/*////////////////////////////////////// + * GLOBAL VARIABLES ******************** + *////////////////////////////////////// -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(); -} +//NA bool toric=false; +bool power_protection = true; +bool grid_points = true; +bool is2d = true; +//FT _sfty = pow(10,-14); +bool torus = false; -void write_points( std::string file_name, std::vector< Point_d > & points) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto w : points) - { - for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - } - ofs.close(); -} - -void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks) -{ - std::ofstream ofs (file_name, std::ofstream::out); - for (auto u: witness_complex.complex_vertex_range()) - for (auto v: witness_complex.complex_vertex_range()) - { - typeVectorVertex edge = {u,v}; - if (u < v && witness_complex.find(edge) != witness_complex.null_simplex()) - { - for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n"; - for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it) - ofs << *it << " "; - ofs << "\n\n\n"; - } - } - ofs.close(); -} - - -void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count) -{ - delaunay.insert(W[chosen_landmark]); - landmarks_ind.push_back(chosen_landmark); - landmark_count++; -} - -bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc) -{ - for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it) - if (*v_it == v) - return true; - return false; -} - -bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, bool is_infinite, const Point_d& p, FT delta) -{ - if (!is_infinite) - // FINITE CASE - { - Sphere_d cs(vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0])); - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) - { - //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); - if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) - { - FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point()); - //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) - if (dist2 >= r*r && dist2 <= r*r+delta*delta) - return true; - } - } - } - else - // INFINITE CASE - { - Delaunay_triangulation::Vertex_iterator v = t.vertices_begin(); - while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end()) - v++; - Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE); - Vector_d orth_v = facet_plane.orthogonal_vector(); - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) - if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end()) - { - std::vector<FT> coords; - Point_d p = v_it->point(); - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - if (!p_is_inside && p_delta_is_inside) - return true; - } - } - return false; -} - - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells) -{ - Euclidean_distance ed; - std::vector<Point_d> vertices; - if (!t.is_infinite(c)) - { - // if the cell is finite, we look if the protection is violated - for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it) - vertices.push_back((*v_it)->point()); - Sphere_d cs( vertices.begin(), vertices.end()); - Point_d center_cs = cs.center(); - FT r = sqrt(ed.transformed_distance(center_cs, vertices[0])); - FT dist2 = ed.transformed_distance(center_cs, p); - // if the new point is inside the protection ball of a non conflicting simplex - //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) - if (dist2 >= r*r && dist2 <= r*r+delta*delta) - return true; - c->tds_data().mark_visited(); - marked_cells.push_back(c); - // if the new point is inside the circumscribing ball : continue violation searching on neughbours - if (dist2 < r*r) - for (int i = 0; i < D+1; ++i) - { - Full_cell_handle next_c = c->neighbor(i); - if (next_c->tds_data().is_clear() && - is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells)) - return true; - } - // if the new point is outside the protection sphere - else - { - // facet f is on the border of the conflict zone : check protection of simplex {p,f} - // the new simplex is guaranteed to be finite - vertices.clear(); vertices.push_back(p); - for (int i = 0; i < D+1; ++i) - if (i != index) - vertices.push_back(parent_cell->vertex(i)->point()); - Delaunay_vertex vertex_to_check; - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!vertex_is_in_full_cell(*vh_it, parent_cell)) - { - vertex_to_check = *vh_it; break; - } - if (new_cell_is_violated(t, vertices, false, vertex_to_check->point(), delta)) - return true; - } - } - else - { - // Inside of the convex hull is + side. Outside is - side. - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!t.is_infinite(*vh_it)) - vertices.push_back((*vh_it)->point()); - Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin(); - while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c)) - v_it++; - Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE); - //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point()); - Vector_d orth_v = facet_plane.orthogonal_vector(); - std::vector<FT> coords; - auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin(); - for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i) - coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length())); - Point_d p_delta = Point_d(coords); - bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p); - bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta); - - if (!p_is_inside && p_delta_is_inside) - return true; - //if the cell is infinite we look at the neighbours regardless - c->tds_data().mark_visited(); - marked_cells.push_back(c); - if (p_is_inside) - for (int i = 0; i < D+1; ++i) - { - Full_cell_handle next_c = c->neighbor(i); - if (next_c->tds_data().is_clear() && - is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells)) - return true; - } - else - { - // facet f is on the border of the conflict zone : check protection of simplex {p,f} - // the new simplex is finite if the parent cell is finite - vertices.clear(); vertices.push_back(p); - bool new_simplex_is_finite = false; - for (int i = 0; i < D+1; ++i) - if (i != index) - { - if (t.is_infinite(parent_cell->vertex(i))) - new_simplex_is_finite = true; - else - vertices.push_back(parent_cell->vertex(i)->point()); - } - Delaunay_vertex vertex_to_check; - for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it) - if (!vertex_is_in_full_cell(*vh_it, parent_cell)) - { - vertex_to_check = *vh_it; break; - } - if (new_cell_is_violated(t, vertices, new_simplex_is_finite, vertex_to_check->point(), delta)) - return true; - } - } - return false; -} - -bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) -{ - Euclidean_distance ed; - 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); - for (Full_cell_handle fc : marked_cells) - fc->tds_data().clear(); - return violation_existing_cells; -} - -bool old_is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta) +bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) { + std::cout << "Start protection verification\n"; 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( 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; - } - t.insert(p, c); - return false; -} - -void write_delaunay_mesh(Delaunay_triangulation& t, const Point_d& p) -{ - std::ofstream ofs ("delaunay.mesh", std::ofstream::out); - int nbV = t.number_of_vertices()+1; - ofs << "MeshVersionFormatted 1\nDimension 2\n"; - ofs << "Vertices\n" << nbV << "\n"; - int ind = 1; //index of a vertex + // 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; - for (auto coord = v_it->point().cartesian_begin(); coord != v_it->point().cartesian_end(); ++coord) - ofs << *coord << " "; - ofs << "508\n"; index_of_vertex[v_it] = ind++; } - for (auto coord = p.cartesian_begin(); coord != p.cartesian_end(); ++coord) - ofs << *coord << " "; - ofs << "208\n"; - /* - int nbFacets = 0; - for (auto ft_it = t.finite_facets_begin(); ft_it != t.finite_facets_end(); ++ft_it) - nbFacets++; - ofs << "\nEdges\n" << nbFacets << "\n\n"; - for (auto ft_it = t.facets_begin(); ft_it != t.facets_end(); ++ft_it) - { - if (t.is_infinite(ft_it)) - continue; - for (auto vh_it = ft_it->vertices_begin(); vh_it != ft_it->vertices_end(); ++vh_it) - ofs << index_of_vertex[*vh_it] << " "; - } - */ - ofs << "Triangles " << t.number_of_finite_full_cells()+1 << "\n"; - for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) - { - if (t.is_infinite(fc_it)) - continue; - for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it) - ofs << index_of_vertex[*vh_it] << " "; - ofs << "508\n"; - } - ofs << nbV << " " << nbV << " " << nbV << " " << 208 << "\n"; - ofs << "End\n"; - ofs.close(); -} - -bool triangulation_is_protected(Delaunay_triangulation& t, FT delta) -{ - // Verification part - Euclidean_distance ed; for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) if (!t.is_infinite(fc_it)) - for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it) - if (!t.is_infinite(v_it)) + { + 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)) - { - 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())); - 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 (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r) - if (dist2 <= r*r+delta*delta && dist2 >= r*r) - { - write_delaunay_mesh(t, v_it->point()); - std::cout << "Problematic vertex " << *v_it << " "; - std::cout << "Problematic cell " << *fc_it << "\n"; - std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", r^2+delta^2 = " << r*r+delta*delta << "\n"; - return false; - } - } - + 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; } -void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind) -{ - for (unsigned j = 0; j < landmarks_ind.size(); ++j) - landmarks.push_back(W[landmarks_ind[j]]); -} - -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; - 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); - } -} +////////////////////////////////////////////////////////////////////////////////////////////////////////// +// SAMPLING RADIUS +////////////////////////////////////////////////////////////////////////////////////////////////////////// -FT sampling_radius(Delaunay_triangulation& t) +FT sampling_radius(Delaunay_triangulation& t, FT epsilon0) { - FT epsilon2 = 4.0; + 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)) @@ -578,134 +220,106 @@ FT sampling_radius(Delaunay_triangulation& t) if (!in_cube) continue; FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin())); - if (epsilon2 > r2) - epsilon2 = r2; + 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 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); + return sampling_radius(t, epsilon0); } -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) +// A little script to make a tikz histogram of epsilon distribution +// Returns the average epsilon +FT epsilon_histogram(Delaunay_triangulation& t, int n) { - 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); - } - for (unsigned pos1 = 0; pos1 < D+1; ++pos1) - { - std::vector<FT> point; - for (unsigned i = 0; i < pos1; ++i) - point.push_back(-1); - if (pos1 != D) - point.push_back(1); - for (unsigned i = pos1+1; i < D; ++i) - point.push_back(0); - assert(point.size() == D); - W[index_list.front()] = Point_d(point); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - } - // add the first D+1 vertices to form one finite cell - /* - for (int i = 0; i <= D+1; ++i) - { - t.insert(W[index_list.front()]); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - } - */ - /* - { - std::vector<FT> coords; - for (int i = 0; i < D; ++i) - coords.push_back(-1); - W[index_list.front()] = Point_d(coords); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - for (int i = 0; i < D; ++i) - { - coords.clear(); - for (int j = 0; j < D; ++j) - if (i == j) - coords.push_back(1); - else - coords.push_back(-1); - W[index_list.front()] = Point_d(coords); - insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count); - index_list.pop_front(); - } - } - */ - //std::cout << t; - //assert(t.number_of_vertices() == D+1); - //assert(landmarks_ind.size() == D+1); - //assert(W[landmarks_ind[0]][0] == 0); - // add other vertices if they don't violate protection - std::list<int>::iterator list_it = index_list.begin(); - while (list_it != index_list.end()) + 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 (!is_violating_protection(W[*list_it], t, D, delta)) - { - // If no conflicts then insert in every copy of T^3 - is_violating_protection(W[*list_it], t, D, delta); - insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count); - index_list.erase(list_it); - list_it = index_list.begin(); - //std::cout << "index_list_size() = " << index_list.size() << "\n"; - } - else - { - list_it++; - //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n"; - } - //write_delaunay_mesh(t, W[*list_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)]++; } - fill_landmarks(W, landmarks, landmarks_ind); - fill_full_cell_vector(t, full_cells); - if (triangulation_is_protected(t, delta)) - std::cout << "Triangulation is ok\n"; - else - std::cout << "Triangulation is BAD!! T_T しくしく!\n"; - write_delaunay_mesh(t, Point_d(std::vector<FT>({0,0}))); - //std::cout << t << std::endl; + 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; } -template <typename T> -void print_vector(std::vector<T> v) +FT epsilon_histogram_by_delaunay(Point_Vector& points, int n) { - std::cout << "["; - if (!v.empty()) - { - std::cout << *(v.begin()); - for (auto it = v.begin()+1; it != v.end(); ++it) - { - std::cout << ","; - std::cout << *it; - } - } - std::cout << "]"; + 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 @@ -764,7 +378,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std } } std::string out_file = "wl_result"; - write_wl(out_file,WL); + //write_wl(out_file,WL); //******************** Constructng a witness complex std::cout << "Entered witness complex construction\n"; @@ -787,7 +401,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std 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; @@ -814,7 +428,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std if (count_bad[i] != 0) std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl; std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl; - */ + //*********************** Perturb bad link landmarks /* for (auto u: perturbL) @@ -848,16 +462,19 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std ofs.close(); } - write_edges("landmarks/edges", witnessComplex, landmarks); + //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] @@ -866,40 +483,98 @@ int main (int argc, char * const argv[]) } int nbP = atoi(argv[1]); int dim = atoi(argv[2]); - double delta = atof(argv[3]); + double theta0 = atof(argv[3]); + //double delta = atof(argv[3]); + + is2d = (dim == 2); std::cout << "Let the carnage begin!\n"; Point_Vector point_vector; - generate_points_random_box(point_vector, nbP, dim); - FT epsilon = point_sampling_radius_by_delaunay(point_vector); + 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(); /* - for (int i = 0; i < 11; 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); - std::cout << "delta = " << delta << std::endl; - L = {}; chosen_landmarks = {}; - 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(); - //write_points("landmarks/landmarks0",L); - } - */ + 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); @@ -909,6 +584,7 @@ int main (int argc, char * const argv[]) 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); + //landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells); timer.reset(); + */ } |