/* 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): Pawel Dlotko
*
* Copyright (C) 2015 INRIA (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 .
*/
//stat part:
#include
#include
#include
#include
#include
//persistence part:
#include
#include
#include
#include
#include
using Persistence_landscape = Gudhi::Gudhi_stat::Persistence_landscape;
typedef int Vertex_handle;
//typedef double Filtration_value;
//if this variable is -1, then the infinite interals are ignored. If not, they infinite values are replaced with what_to_replace_infinite_intervals_with:
double what_to_replace_infinite_intervals_with = -1;
class compute_persistence_landscape_of_a_point_cloud_in_certain_dimension
{
public:
compute_persistence_landscape_of_a_point_cloud_in_certain_dimension( std::vector< std::vector< double > >& points_ , int dimension , double threshold_ , int coeficient_field_ = 11 , double min_persistence_ = 0 ):dim( dimension ),points(points_),threshold(threshold_),coeficient_field(coeficient_field_),min_persistence(min_persistence_){}
//This function takes a vector of indices (numbers_to_sample). It will select the points from this->points having those indices, construct Rips complex and persistence intervals based on this.
//Then it will filter the intervals to find only those in the dimension this->dim, and construct a persistence landascape based on this. Thie will be the result of the procedure.
Persistence_landscape operator()( std::vector< size_t > numbers_to_sample )
{
bool dbg = false;
//take the subsampled points:
std::vector< std::vector< double > > points_in_subsample;
points_in_subsample.reserve( numbers_to_sample.size() );
for ( size_t i = 0 ; i != numbers_to_sample.size() ; ++i )
{
points_in_subsample.push_back( this->points[ numbers_to_sample[i] ] );
}
using Stree = Gudhi::Simplex_tree;
using Filtration_value = Stree::Filtration_value;
using Rips_complex = Gudhi::rips_complex::Rips_complex;
//construct a Rips complex based on it and compute its persistence:
Rips_complex rips_complex(points_in_subsample, this->threshold, Euclidean_distance());
// Construct the Rips complex in a Simplex Tree
Stree st;
// expand the graph until dimension dim_max
rips_complex.create_complex(st, this->dim + 1);
// Compute the persistence diagram of the complex
Gudhi::persistent_cohomology::Persistent_cohomology pcoh(st);
// initializes the coefficient field for homology
pcoh.init_coefficients( this->coeficient_field );
pcoh.compute_persistent_cohomology(this->min_persistence);
auto persistence_pairs = pcoh.get_persistent_pairs();
//From the persistence take only this in the dimension this->dim:
if ( dbg )std::cerr << "Here are the persistence pairs :\n";
std::vector< std::pair< double,double > > persistence_in_fixed_dimension;
for ( size_t i = 0 ; i != persistence_pairs.size() ; ++i )
{
if ( st.dimension( std::get<0>(persistence_pairs[i]) ) == this->dim )
{
double birth = st.filtration( std::get<0>(persistence_pairs[i]) );
double death = st.filtration( std::get<1>(persistence_pairs[i]) );
if ( std::get<1>(persistence_pairs[i]) != st.null_simplex() )
{
//finite interval
persistence_in_fixed_dimension.push_back( std::pair( birth , death ) );
if (dbg){std::cout << "birth : " << birth << " , death : " << death << std::endl;}
}
else
{
//infinite interval
if ( what_to_replace_infinite_intervals_with != -1 )
{
persistence_in_fixed_dimension.push_back( std::pair( birth , what_to_replace_infinite_intervals_with ) );
if (dbg){std::cout << "birth : " << birth << " , death : " << what_to_replace_infinite_intervals_with << std::endl;}
}
}
}
}
if ( dbg )std::cerr << "Persistence pairs computed \n";
//Construct and return the persistence landscape:
return Persistence_landscape( persistence_in_fixed_dimension );
}
private:
int dim;
std::vector< std::vector< double > >& points;
double threshold;
int coeficient_field;
double min_persistence;
};
class distance_between_landscapes
{
public:
distance_between_landscapes( double exponent_ ):exponent(exponent_){}
double operator()( const Persistence_landscape& first , const Persistence_landscape& second )
{
return first.distance( second, this->exponent );
}
private:
double exponent;
};
int main( int argc , char** argv )
{
std::cout << "The parameters of this program are : " << std::endl;
std::cout << "(1) a name of a file with points," << std:: endl;
std::cout << "(2) a number of repetitions of bootstrap (integer)," << std::endl;
std::cout << "(3) a size of subsample (integer, smaller than the number of points. " << std::endl;
std::cout << "(4) An real value p such that L^p distance is going to be computed. \n";
std::cout << "(5) A dimension of persistence that is to be taken into account (positive integer) \n";
std::cout << "(6) A maximal diameter to which complex is to be grown (positive integer) \n";
std::cout << "(d) a quantile (real number between 0 and 1. If you do not know what to set, set it to 0.95." << std::endl;
if ( argc != 8 )
{
std::cerr << "Wrong number of parameters, the program will now terminate.\n";
return 1;
}
const char* filename = argv[1];
size_t number_of_repetitions_of_bootstrap = (size_t)atoi( argv[2] );
size_t size_of_subsample = (size_t)atoi( argv[3] );
double p = atoi( argv[4] );
int dimension = atoi( argv[5] );
double threshold = atof( argv[6] );
double quantile = atof( argv[7] );
std::cout << "Now we will read points from the file : " << filename << " and then perform " << number_of_repetitions_of_bootstrap << " times the bootstrap on it by choosing subsample of a size " << size_of_subsample << std::endl;
std::vector< std::vector< double > > points = Gudhi::Gudhi_stat::read_numbers_from_file_line_by_line( filename );
std::cout << "Read : " << points.size() << " points.\n";
distance_between_landscapes distance( p );//L^p distance.
compute_persistence_landscape_of_a_point_cloud_in_certain_dimension characteristic_fun( points , dimension , threshold );
//and now we can run the real bootstrap.
//template < typename PointCloudCharacteristics , typename CharacteristicFunction , typename DistanceBetweenPointsCharacteristics >
//In this case, the PointCloudCharacteristics is just a vector of numbers of points (in a order fixed on points vector).
//CharacteristicFunction is just identity, transforming std::vector< size_t > to itself.
//DistanceBetweenPointsCharacteristics is the place were all happens. This class hace the information about the coordinates of the points, and allows to compute a Hausdorff distance between
//the collection of all points, and the subsample.
double result = Gudhi::Gudhi_stat::bootstrap<
Persistence_landscape , //PointCloudCharacteristics, persistence landascapes constructed based on vector of
//pairs of birth--death values in a cartain dimension.
compute_persistence_landscape_of_a_point_cloud_in_certain_dimension , //CharacteristicFunction, in this case, we will need to compute persistence in a certain dimension.
distance_between_landscapes //DistanceBetweenPointsCharacteristics. In this case
>
( points.size() , characteristic_fun , distance , number_of_repetitions_of_bootstrap , size_of_subsample , quantile );
std::cout << "result of bootstrap : " << result << std::endl;
return 0;
}