From 68752bc0ce16dbff783b5f84a2d02a10b7d05a4e Mon Sep 17 00:00:00 2001 From: skachano Date: Thu, 9 Jun 2016 12:52:35 +0000 Subject: Added everything missing git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1265 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: c148ce9136bb0786acc1c9ae49827cb3958326c6 --- src/Witness_complex/example/bench_rwit.cpp | 616 +++++++++++++++++++++++++++++ 1 file changed, 616 insertions(+) create mode 100644 src/Witness_complex/example/bench_rwit.cpp (limited to 'src/Witness_complex/example/bench_rwit.cpp') diff --git a/src/Witness_complex/example/bench_rwit.cpp b/src/Witness_complex/example/bench_rwit.cpp new file mode 100644 index 00000000..f83f9f42 --- /dev/null +++ b/src/Witness_complex/example/bench_rwit.cpp @@ -0,0 +1,616 @@ +/* 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) 2016 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include "Landmark_choice_random_knn.h" +#include "Landmark_choice_sparsification.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "generators.h" +#include "output.h" +#include "output_tikz.h" + +using namespace Gudhi; +using namespace Gudhi::witness_complex; +using namespace Gudhi::persistent_cohomology; + +typedef std::vector Point_Vector; +typedef Simplex_tree STree; +//typedef Simplex_tree<> STree; +typedef A0_complex< STree> A0Complex; +typedef STree::Simplex_handle Simplex_handle; + +typedef A0_complex< STree > SRWit; +typedef Relaxed_witness_complex< STree > WRWit; + +/** + * \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 rips(Point_Vector & points, double alpha2, int dim_max, STree& st) +{ + Graph_t prox_graph = compute_proximity_graph(points, sqrt(alpha2) + , euclidean_distance >); + // Construct the Rips complex in a Simplex Tree + // insert the proximity graph in the simplex tree + st.insert_graph(prox_graph); + // expand the graph until dimension dim_max + st.expansion(dim_max); +} + +void output_experiment_information(char * const file_name) +{ + std::cout << "Enter a valid experiment number. Usage: " + << file_name << " exp_no options\n"; + std::cout << "Experiment description:\n" + << "0 nbP nbL dim alpha limD mu_epsilon: " + << "Build persistence diagram on relaxed witness complex " + << "built from a point cloud on (dim-1)-dimensional sphere " + << "consisting of nbP witnesses and nbL landmarks. " + << "The maximal relaxation is alpha and the limit on simplicial complex " + << "dimension is limD.\n"; + std::cout << "1 file_name nbL alpha limD: " + << "Build persistence diagram on relaxed witness complex " + << "build from a point cloud stored in a file and nbL landmarks. " + << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n"; +} + +void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT mu_epsilon = 0.1, + std::string mesh_filename = "witness") +{ + clock_t start, end; + STree simplex_tree; + + // Choose landmarks + std::vector > knn; + std::vector > distances; + start = clock(); + //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances); + + std::vector landmarks; + Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses + landmarks, // aka landmarks + alpha2, + limD, + knn, + distances); + end = clock(); + double time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Choice of " << nbL << " landmarks took " + << time << " s. \n"; + // Compute witness complex + start = clock(); + A0Complex rw(distances, + knn, + simplex_tree, + nbL, + alpha2, + limD); + end = clock(); + time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Witness complex for " << nbL << " landmarks took " + << time << " s. \n"; + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + //std::cout << simplex_tree << "\n"; + + // Compute the persistence diagram of the complex + simplex_tree.set_dimension(limD); + persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh(simplex_tree, true); + int p = 3; + pcoh.init_coefficients( p ); //initilizes the coefficient field for homology + start = clock(); + pcoh.compute_persistent_cohomology( alpha2/10 ); + end = clock(); + time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Persistence diagram took " + << time << " s. \n"; + pcoh.output_diagram(); + + int chi = 0; + for (auto sh: simplex_tree.complex_simplex_range()) + chi += 1-2*(simplex_tree.dimension(sh)%2); + std::cout << "Euler characteristic is " << chi << std::endl; + + Gudhi::witness_complex::Dim_lists simplices(simplex_tree, limD); + + // std::vector simplices; + simplices.collapse(); + simplices.output_simplices(); + + STree collapsed_tree; + for (auto sh: simplices) { + std::vector vertices; + for (int v: collapsed_tree.simplex_vertex_range(sh)) + vertices.push_back(v); + collapsed_tree.insert_simplex(vertices); + } + std::vector landmarks_ind(nbL); + for (unsigned i = 0; i != distances.size(); ++i) { + if (distances[i][0] == 0) + landmarks_ind[knn[i][0]] = i; + } + //write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplices, false, true); + write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, mesh_filename+"_before_collapse.mesh"); + + collapsed_tree.set_dimension(limD); + persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh2(collapsed_tree, true); + pcoh2.init_coefficients( p ); //initilizes the coefficient field for homology + pcoh2.compute_persistent_cohomology( alpha2/10 ); + pcoh2.output_diagram(); + + chi = 0; + for (auto sh: simplices) + chi += 1-2*(simplex_tree.dimension(sh)%2); + std::cout << "Euler characteristic is " << chi << std::endl; + write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, mesh_filename+"_after_collapse.mesh"); +} + +void rips_experiment(Point_Vector & points, double threshold, int dim_max) +{ + typedef STree ST; + clock_t start, end; + ST st; + + // Compute the proximity graph of the points + start = clock(); + rips(points, threshold, dim_max, st); + end = clock(); + + double time = static_cast(end - start) / CLOCKS_PER_SEC; + std::cout << "Rips complex took " + << time << " s. \n"; + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + //std::cout << " and has dimension " << st.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + persistent_cohomology::Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + int p = 3; + double min_persistence = -1; //threshold/5; + pcoh.init_coefficients(p); + pcoh.compute_persistent_cohomology(min_persistence); + pcoh.output_diagram(); +} + + +int experiment0 (int argc, char * const argv[]) +{ + if (argc != 8) { + std::cerr << "Usage: " << argv[0] + << " 0 nbP nbL dim alpha limD mu_epsilon\n"; + return 0; + } + /* + boost::filesystem::path p; + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + + int nbP = atoi(argv[2]); + int nbL = atoi(argv[3]); + int dim = atoi(argv[4]); + double alpha = atof(argv[5]); + int limD = atoi(argv[6]); + double mu_epsilon = atof(argv[7]); + + // Read the point file + Point_Vector point_vector; + generate_points_sphere(point_vector, nbP, dim); + std::cout << "Successfully generated " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + rw_experiment(point_vector, nbL, alpha, limD); + return 0; +} + +int experiment1 (int argc, char * const argv[]) +{ + if (argc != 8) { + std::cerr << "Usage: " << argv[0] + << " 1 file_name nbL alpha mu_epsilon limD experiment_name\n"; + return 0; + } + /* + boost::filesystem::path p; + for (; argc > 2; --argc, ++argv) + p /= argv[1]; + */ + + std::string file_name = argv[2]; + int nbL = atoi(argv[3]), limD = atoi(argv[6]); + double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]); + std::string experiment_name = argv[7]; + + // Read the point file + Point_Vector point_vector; + read_points_cust(file_name, point_vector); + std::cout << "The file contains " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + STree simplex_tree; + std::vector > knn; + std::vector > distances; + std::vector landmarks; + Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses + landmarks, // aka landmarks + alpha2, + limD, + knn, + distances); + WRWit rw(distances, + knn, + simplex_tree, + nbL, + 0, + limD); + std::vector landmarks_ind(nbL); + for (unsigned i = 0; i != distances.size(); ++i) { + if (distances[i][0] == 0) + landmarks_ind[knn[i][0]] = i; + } + + write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, experiment_name+"_witness_complex.mesh"); + + int chi = 0; + for (auto sh: simplex_tree.complex_simplex_range()) + chi += 1-2*(simplex_tree.dimension(sh)%2); + std::cout << "Euler characteristic is " << chi << std::endl; + + rw_experiment(point_vector, nbL, alpha2, limD, mu_epsilon, experiment_name); + + // ok = false; + // while (!ok) { + // std::cout << "Rips complex: parameters threshold, limD.\n"; + // std::cout << "Enter threshold: "; + // std::cin >> alpha; + // std::cout << "Enter limD: "; + // std::cin >> limD; + // std::cout << "Start Rips complex...\n"; + // rips_experiment(point_vector, alpha, limD); + // std::cout << "Is the result correct? [y/n]: "; + // char answer; + // std::cin >> answer; + // switch (answer) { + // case 'n': + // ok = false; break; + // default : + // ok = true; break; + // } + // } + return 0; +} + +/******************************************************************************************** + * Length of the good interval experiment + *******************************************************************************************/ + +struct Pers_endpoint { + double alpha; + bool start; + int dim; + Pers_endpoint(double alpha_, bool start_, int dim_) + : alpha(alpha_), start(start_), dim(dim_) + {} +}; + +/* +struct less_than_key { + inline bool operator() (const MyStruct& struct1, const MyStruct& struct2) { + return (struct1.key < struct2.key); + } +}; +*/ + +double good_interval_length(const std::vector & desired_homology, STree & simplex_tree, double alpha2) +{ + int nbL = simplex_tree.num_vertices(); + int p = 3; + persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh(simplex_tree, true); + pcoh.init_coefficients( p ); //initilizes the coefficient field for homology + pcoh.compute_persistent_cohomology( -1 ); + std::ofstream out_stream("pers_diag.tmp"); + pcoh.output_diagram(out_stream); + out_stream.close(); + std::ifstream in_stream("pers_diag.tmp", std::ios::in); + std::string line; + std::vector pers_endpoints; + while (getline(in_stream, line)) { + int p, dim; + double alpha_start, alpha_end; + std::istringstream iss(line); + iss >> p >> dim >> alpha_start >> alpha_end; + if (alpha_start != alpha_end) { + if (alpha_end < alpha_start) + alpha_end = alpha2; + pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim)); + pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim)); + } + } + std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl; + in_stream.close(); + std::sort(pers_endpoints.begin(), + pers_endpoints.end(), + [](const Pers_endpoint & p1, const Pers_endpoint & p2){ + return p1.alpha < p2.alpha;} + ); + write_barcodes("pers_diag.tmp", alpha2); + /* + for (auto p: pers_endpoints) { + std::cout << p.alpha << " " << p.dim << " " << p.start << "\n"; + } + */ + std::vector current_homology(nbL-1,0); + current_homology[0] = 1; // for the compulsary "0 0 inf" entry + double good_start = 0, good_end = 0; + double sum_intervals = 0; + int num_pieces = 0; + bool interval_in_process = (desired_homology == current_homology); + for (auto p: pers_endpoints) { + /* + std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start + << " ["; + for (int v: current_homology) + std::cout << v << " "; + std::cout << "]\n"; + */ + if (p.start) + current_homology[p.dim]++; + else + current_homology[p.dim]--; + if (interval_in_process) { + good_end = p.alpha; + sum_intervals += good_end - good_start; + std::cout << "good_start = " << good_start + << ", good_end = " << good_end << "\n"; + + Gudhi::witness_complex::Dim_lists simplices(simplex_tree, nbL-1, (good_end - good_start)/2); + //simplices.collapse(); + //simplices.output_simplices(); + interval_in_process = false; + //break; + } + else if (desired_homology == current_homology) { + interval_in_process = true; + good_start = p.alpha; + num_pieces++; + } + } + std::cout << "Number of good homology intervals: " << num_pieces << "\n"; + return sum_intervals; +} + + + +void run_comparison(std::vector > const & knn, + std::vector > const & distances, + Point_Vector & points, + unsigned nbL, + double alpha2, + std::vector& desired_homology) +{ + clock_t start, end; + STree simplex_tree; + + alpha2 = 0.55; + std::cout << "alpha2 = " << alpha2 << "\n"; + start = clock(); + SRWit srwit(distances, + knn, + simplex_tree, + nbL, + alpha2, + nbL-1); + end = clock(); + std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl; + simplex_tree.set_dimension(nbL-1); + + std::cout << "Good homology interval length for SRWit is " + << good_interval_length(desired_homology, simplex_tree, alpha2) << "\n"; + std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; + + STree simplex_tree2; + alpha2 = 0.35; + std::cout << "alpha2 = " << alpha2 << "\n"; + start = clock(); + WRWit wrwit(distances, + knn, + simplex_tree2, + nbL, + alpha2, + nbL-1); + end = clock(); + std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl; + simplex_tree.set_dimension(nbL-1); + + std::cout << "Good homology interval length for WRWit is " + << good_interval_length(desired_homology, simplex_tree2, alpha2) << "\n"; + std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; + + STree simplex_tree3; + std::cout << "Enter alpha2 for Rips: "; + std::cin >> alpha2; + std::cout << "\n"; + std::cout << "points.size() = " << points.size() << "\n"; + start = clock(); + rips(points, + 0.2, + nbL-1, + simplex_tree3); + end = clock(); + std::cout << "Rips.size = " << simplex_tree3.num_simplices() << std::endl; + simplex_tree.set_dimension(nbL-1); + + std::cout << "Good homology interval length for WRWit is " + << good_interval_length(desired_homology, simplex_tree3, alpha2) << "\n"; + std::cout << "Time: " << static_cast(end - start) / CLOCKS_PER_SEC << " s. \n"; +} + +int experiment2(int argc, char * const argv[]) +{ + for (unsigned d = 3; d < 4; d++) { + // Sphere S^d + Point_Vector point_vector; + unsigned N = 1; + double alpha2 = 2.4 - 0.4*d; + switch (d) { + case 1: alpha2 = 2.2; break; + case 2: alpha2 = 1.7; break; + case 3: alpha2 = 1.5; break; + case 4: alpha2 = 1.4; break; + default: alpha2 = 1.4; break; + } + std::cout << "alpha2 = " << alpha2 << "\n"; + unsigned nbL = 20; + std::vector desired_homology(nbL-1,0); + desired_homology[0] = 1; desired_homology[d] = 1; + + + for (unsigned i = 1; i <= N; ++i) { + unsigned nbW = 1000*i;//, nbL = 20; + double mu_epsilon = 1/sqrt(nbL); + std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl; + generate_points_sphere(point_vector, i*1000, d+1); + std::vector landmarks; + Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + + std::vector > knn; + std::vector > distances; + + std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; + Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses + landmarks, // aka landmarks + alpha2, + nbL-1, + knn, + distances); + run_comparison(knn, distances, point_vector, nbL, alpha2, desired_homology); + } + } + /* + { + // SO(3) + Point_Vector point_vector; + double alpha2 = 0.6; + std::cout << "alpha2 = " << alpha2 << "\n"; + unsigned nbL = 150; + std::vector desired_homology(nbL-1,0); + desired_homology[0] = 1; desired_homology[1] = 1; desired_homology[2] = 1; //Kl + // desired_homology[0] = 1; desired_homology[3] = 1; //SO3 + + double mu_epsilon = 1/sqrt(nbL); + if (argc < 3) std::cerr << "No file name indicated!\n"; + read_points_cust(argv[2], point_vector); + int nbW = point_vector.size(); + std::cout << "Running test SO(3), |W|=" << nbW << ", |L|=" << nbL << std::endl; + std::vector landmarks; + Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks); + + std::vector > knn; + std::vector > distances; + + std::cout << "|L| after sparsification: " << landmarks.size() << "\n"; + Gudhi::witness_complex::build_distance_matrix(point_vector, // aka witnesses + landmarks, // aka landmarks + alpha2, + nbL-1, + knn, + distances); + run_comparison(knn, distances, point_vector, nbL, alpha2, desired_homology); + } + */ + return 0; +} + +int experiment3(int argc, char * const argv[]) +{ + // COLLAPSES EXPERIMENT + + return 0; +} + +int main (int argc, char * const argv[]) +{ + if (argc == 1) { + output_experiment_information(argv[0]); + return 1; + } + switch (atoi(argv[1])) { + case 0 : + return experiment0(argc, argv); + break; + case 1 : + return experiment1(argc, argv); + break; + case 2 : + return experiment2(argc, argv); + break; + default : + output_experiment_information(argv[0]); + return 1; + } +} -- cgit v1.2.3