/* 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;
}
}