diff options
23 files changed, 1174 insertions, 454 deletions
diff --git a/src/Gudhi_stat/example/vector_representation.cpp b/src/Gudhi_stat/example/vector_representation.cpp index 835b52f9..bb4b5d01 100644 --- a/src/Gudhi_stat/example/vector_representation.cpp +++ b/src/Gudhi_stat/example/vector_representation.cpp @@ -57,16 +57,16 @@ int main( int argc , char** argv ) persistence2.push_back( std::make_pair(6,10) ); //create two persistence vectors based on persistence1 and persistence2: - Vector_distances_in_diagram<euclidean_distance<double> > v1 = Vector_distances_in_diagram<euclidean_distance<double> >( persistence1 , std::numeric_limits<size_t>::max() ); - Vector_distances_in_diagram<euclidean_distance<double> > v2 = Vector_distances_in_diagram<euclidean_distance<double> >( persistence2 , std::numeric_limits<size_t>::max() ); + Vector_distances_in_diagram<Euclidean_distance<double> > v1 = Vector_distances_in_diagram<Euclidean_distance<double> >( persistence1 , std::numeric_limits<size_t>::max() ); + Vector_distances_in_diagram<Euclidean_distance<double> > v2 = Vector_distances_in_diagram<Euclidean_distance<double> >( persistence2 , std::numeric_limits<size_t>::max() ); //writing to a stream: std::cout << "v1 : " << v1 << std::endl; std::cout << "v2 : " << v2 << std::endl; //averages: - Vector_distances_in_diagram<euclidean_distance<double> > average; - std::vector< Vector_distances_in_diagram<euclidean_distance<double> >* > to_average; + Vector_distances_in_diagram<Euclidean_distance<double> > average; + std::vector< Vector_distances_in_diagram<Euclidean_distance<double> >* > to_average; to_average.push_back( &v1 ); to_average.push_back( &v2 ); average.compute_average( to_average ); @@ -94,7 +94,7 @@ int main( int argc , char** argv ) return 1; } - Vector_distances_in_diagram< euclidean_distance<double> > p( argv[1] , 100 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p( argv[1] , 100 ); cout << "This is a vector corresponding to the input persistence diagram : \n"; cout << p << endl; @@ -102,7 +102,7 @@ int main( int argc , char** argv ) if ( argc == 3 ) { - Vector_distances_in_diagram< euclidean_distance<double> > p_prime( argv[2] , 100); + Vector_distances_in_diagram< Euclidean_distance<double> > p_prime( argv[2] , 100); cout << "p_prime : " <<p_prime << endl; @@ -111,7 +111,7 @@ int main( int argc , char** argv ) to_average.push_back( (Abs_Topological_data_with_averages*)(&p) ); to_average.push_back( (Abs_Topological_data_with_averages*)(&p_prime) ); - Vector_distances_in_diagram< euclidean_distance<double> > average; + Vector_distances_in_diagram< Euclidean_distance<double> > average; average.compute_average( to_average ); cout << "Here is an average : " << average << endl; diff --git a/src/Gudhi_stat/include/gudhi/Hausdorff_distances.h b/src/Gudhi_stat/include/gudhi/Hausdorff_distances.h new file mode 100644 index 00000000..de39ad7d --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Hausdorff_distances.h @@ -0,0 +1,142 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef HAUSDORFF_DISTANCES_H +#define HAUSDORFF_DISTANCES_H + +#include <cmath> +#include <limits> +#include <vector> +#include <cstdlib> +#include <iostream> + + +/** + * This file contains various implementations of Hausrodff distances between collections of points. It contains various implementations that can be used for specific case. +**/ + +/** + * The implementation below works for a case of a metric space given by a distance matrix, and a subspace of a metric space. Our task is to find a Hausdorff distance between the subspace and the whole space. + * The input is a distance matrix (lower triangular part) and a vector of bools indicating a subspace (elementss set to true belongs to the space, the elements set to false do not). +**/ +class Hausdorff_distance_between_subspace_and_the_whole_metric_space +{ +public: + Hausdorff_distance_between_subspace_and_the_whole_metric_space( const std::vector< std::vector<double> >& distance_matrix ):distance_matrix(distance_matrix){} + double operator()( const std::vector< bool >& is_subspace ) + { + double maximal_distance = -std::numeric_limits<double>::max(); + for ( size_t j = 0 ; j != this->distance_matrix.size() ; ++j ) + { + double minimal_distance = std::numeric_limits<double>::max(); + for ( size_t i = 0 ; i != this->distance_matrix.size() ; ++i ) + { + if ( !is_subspace[i] )continue; + double distance = 0; + if ( i < j ) + { + distance = this->distance_matrix[i][j]; + } + else + { + if ( i > j )distance = this->distance_matrix[j][i]; + } + if ( distance > minimal_distance )minimal_distance = distance; + } + if ( maximal_distance < minimal_distance )maximal_distance = minimal_distance; + } + return maximal_distance; + }//Hausdorff_distance_between_subspace_and_the_whole_metric_space + + double operator()( const std::vector< size_t >& subspace , const std::vector< size_t >& space = std::vector< size_t >() ) + { + bool dbg = false; + if ( dbg ) + { + std::cerr << "Calling double operator()( const std::vector< size_t >& subspace , const std::vector< size_t >& space = std::vector< size_t >() ) method \n"; + std::cerr << "subspace.size() : " << subspace.size() << std::endl; + } + double maximal_distance = -std::numeric_limits<double>::max(); + for ( size_t j = 0 ; j != this->distance_matrix.size() ; ++j ) + { + double minimal_distance = std::numeric_limits<double>::max(); + for ( size_t i = 0 ; i != subspace.size() ; ++i ) + { + double distance = 0; + if ( subspace[i] < j ) + { + distance = this->distance_matrix[ j ][ subspace[i] ]; + } + else + { + if ( subspace[i] > j )distance = this->distance_matrix[ subspace[i] ][ j ]; + } + if ( distance < minimal_distance )minimal_distance = distance; + } + if ( maximal_distance < minimal_distance )maximal_distance = minimal_distance; + } + if ( dbg ) + { + std::cerr << "maximal_distance : " << maximal_distance << std::endl; + getchar(); + } + return maximal_distance; + }//Hausdorff_distance_between_subspace_and_the_whole_metric_space +private: + const std::vector< std::vector<double> >& distance_matrix; +}; + +template <typename PointType , typename distanceFunction > +std::vector< std::vector<double> > compute_all_to_all_distance_matrix_between_points( const std::vector< PointType >& points ) +{ + bool dbg = false; + std::vector< std::vector<double> > result; + result.reserve( points.size() ); + distanceFunction f; + for ( size_t i = 0 ; i != points.size() ; ++i ) + { + std::vector<double> this_row; + this_row.reserve( i ); + for ( size_t j = 0 ; j != i ; ++j ) + { + double distance = f( points[i] , points[j] ); + if ( dbg ){std::cerr << "The distance between point : " << i << " and the point : " << j << " is :" << distance << std::endl;} + this_row.push_back( distance ); + } + result.push_back( this_row ); + } + return result; +}//compute_all_to_all_distance_matrix_between_points + +template <typename T> +class identity +{ +public: + T& operator()( T& org ) + { + return org; + } +}; + + + +#endif diff --git a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Real_valued_topological_data.h b/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Real_valued_topological_data.h deleted file mode 100644 index 09c2cf85..00000000 --- a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Real_valued_topological_data.h +++ /dev/null @@ -1,54 +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): 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 <http://www.gnu.org/licenses/>. - */ - -#ifndef Abs_Real_valued_topological_data_H_ -#define Abs_Real_valued_topological_data_H_ - -#include <gudhi/abstract_classes/Abs_Topological_data.h> - -namespace Gudhi -{ -namespace Gudhi_stat -{ - - - -/** -* This is specialization of a topological data class allows computing various real-valued characterizations of topological data. -**/ - -class Abs_Real_valued_topological_data : public virtual Abs_Topological_data -{ -public: - Abs_Real_valued_topological_data():number_of_functions_for_projections_to_reals(0){} - Abs_Real_valued_topological_data( size_t number_of_functions_ ):number_of_functions_for_projections_to_reals(number_of_functions_){} - size_t number_of_projections_to_R(){return this->number_of_functions_for_projections_to_reals;}; - virtual double project_to_R( int number_of_function ) = 0; - virtual ~Abs_Real_valued_topological_data(){} -protected: - size_t number_of_functions_for_projections_to_reals; -}; - -}//namespace Gudhi_stat -}//namespace Gudhi - -#endif diff --git a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data.h b/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data.h deleted file mode 100644 index e38ae2df..00000000 --- a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data.h +++ /dev/null @@ -1,47 +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): 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 <http://www.gnu.org/licenses/>. - */ - -#ifndef ABS_Topological_data_H_ -#define ABS_Topological_data_H_ - - - -namespace Gudhi -{ -namespace Gudhi_stat -{ - -/** -* This is an abstract container to store topological information. Most typically, this information will be some representation of persistent homology. -**/ - -class Abs_Topological_data -{ -//public: -// Abs_Topological_data(){}; -//protected: -}; - -}//namespace Gudhi_stat -}//namespace Gudhi - -#endif diff --git a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_averages.h b/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_averages.h deleted file mode 100644 index cc2cd00f..00000000 --- a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_averages.h +++ /dev/null @@ -1,62 +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): 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 <http://www.gnu.org/licenses/>. - */ - -#ifndef Abs_Topological_data_with_averages_H_ -#define Abs_Topological_data_with_averages_H_ -#include <gudhi/abstract_classes/Abs_Topological_data.h> -#include <vector> - - -namespace Gudhi -{ -namespace Gudhi_stat -{ - -/** -* This is an abstract container to store topological information. Most typically, this information will be some representation of persistent homology. -**/ - -class Abs_Topological_data_with_averages : public virtual Abs_Topological_data -{ -public: - Abs_Topological_data_with_averages(){}; - - //options: - //1) Function returns void, take std::vector< Abs_Topological_data_with_averages* > and change the object underneeth to the average of this vector. - //2) Function returns Abs_Topological_data*, and compute the average of all the objects in the vector, plus (*this). - //virtual Abs_Topological_data* compute_average( std::vector< Abs_Topological_data_with_averages* > to_average ) = 0; - //At the moment I will try to implement option (1). - virtual void compute_average( const std::vector< Abs_Topological_data_with_averages* >& to_average ) = 0; - virtual ~Abs_Topological_data_with_averages(){} -protected: -}; - - -}//namespace Gudhi_stat -}//namespace Gudhi - - - - - -#endif diff --git a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_distances.h b/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_distances.h deleted file mode 100644 index 5b683268..00000000 --- a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_distances.h +++ /dev/null @@ -1,52 +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): 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 <http://www.gnu.org/licenses/>. - */ - -#ifndef Abs_Topological_data_with_distances_H_ -#define Abs_Topological_data_with_distances_H_ - -//#include "Abs_Topological_data.h" -#include <gudhi/abstract_classes/Abs_Topological_data.h> - - -namespace Gudhi -{ -namespace Gudhi_stat -{ - - -/** -* This is an abstract container to store topological information. Most typically, this information will be some representation of persistent homology. -**/ - -class Abs_Topological_data_with_distances : public virtual Abs_Topological_data -{ -public: - Abs_Topological_data_with_distances(){}; - virtual double distance( const Abs_Topological_data_with_distances* second , double power = 1) = 0; - virtual ~Abs_Topological_data_with_distances(){} -protected: -}; - -}//namespace Gudhi_stat -}//namespace Gudhi - -#endif diff --git a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_scalar_product.h b/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_scalar_product.h deleted file mode 100644 index ccda5e63..00000000 --- a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Topological_data_with_scalar_product.h +++ /dev/null @@ -1,50 +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): 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 <http://www.gnu.org/licenses/>. - */ - -#ifndef Abs_Topological_data_with_scalar_product_H_ -#define Abs_Topological_data_with_scalar_product_H_ - -#include <gudhi/abstract_classes/Abs_Topological_data.h> -#include <vector> - -namespace Gudhi -{ -namespace Gudhi_stat -{ - -/** -* This is an abstract container to store topological information. Most typically, this information will be some representation of persistent homology. -**/ - -class Abs_Topological_data_with_scalar_product : public virtual Abs_Topological_data -{ -public: - Abs_Topological_data_with_scalar_product(){}; - virtual double compute_scalar_product( const Abs_Topological_data_with_scalar_product* second ) = 0; - virtual ~Abs_Topological_data_with_scalar_product(){} -protected: -}; - -}//namespace Gudhi_stat -}//namespace Gudhi - -#endif diff --git a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Vectorized_topological_data.h b/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Vectorized_topological_data.h deleted file mode 100644 index 2f352428..00000000 --- a/src/Gudhi_stat/include/gudhi/abstract_classes/Abs_Vectorized_topological_data.h +++ /dev/null @@ -1,56 +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): 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 <http://www.gnu.org/licenses/>. - */ - -#ifndef Abs_Vectorized_topological_data_H_ -#define Abs_Vectorized_topological_data_H_ - -#include <gudhi/abstract_classes/Abs_Topological_data.h> -#include "Abs_Topological_data.h" -#include <vector> - -namespace Gudhi -{ -namespace Gudhi_stat -{ -using namespace std; - -/** -* This is specialization of a topological data class that allows to create vectors based on topological data. -**/ - -class Abs_Vectorized_topological_data : public virtual Abs_Topological_data -{ -public: - Abs_Vectorized_topological_data():where_to_cut(10){} - Abs_Vectorized_topological_data( size_t where_to_cut ):where_to_cut(where_to_cut){} - size_t number_of_vectorize_functions(){return number_of_functions_for_vectorization;} - virtual std::vector<double> vectorize( int number_of_function ) = 0; - virtual ~Abs_Vectorized_topological_data(){} -protected: - size_t number_of_functions_for_vectorization; - size_t where_to_cut; -}; - -}//namespace Gudhi_stat -}//namespace Gudhi - -#endif diff --git a/src/Gudhi_stat/include/gudhi/bootstrap.h b/src/Gudhi_stat/include/gudhi/bootstrap.h new file mode 100644 index 00000000..ca40f9c7 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/bootstrap.h @@ -0,0 +1,143 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef BOOTSTRAP_H +#define BOOTSTRAP_H + +//concretizations +#include <gudhi/concretizations/Vector_distances_in_diagram.h> +#include <gudhi/concretizations/Persistence_landscape.h> +#include <gudhi/concretizations/Persistence_landscape_on_grid.h> +#include <gudhi/concretizations/Persistence_heat_maps.h> + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_sort.h> +#endif + +#include <vector> +#include <algorithm> +#include <omp.h> +#include <random> +#include <ctime> + +/** +* This is a generic function to perform bootstrap. +* In this function we assume that there is a class to compute characteristic of collection of points (PointCloudCharacteristics) and that it stores coordinates of all points. It allows to compute the characteristic +* of the whole point cloud (by using CharacteristicFunction) or of it proper subset of the whole point cloud (given the list of numers of points in the subset). +* Both functionalities will be used in this implementation. +* The characteristic of point cloud, may be the poit cloud itself, its persistence diagram in a fixed dimension, or anything else. We only assume that space of points characteristics is a metric space +* and that we can compute a distance between two characteristics of collections of points by using DistanceBetweenPointsCharacteristics function. +**/ + + +template < typename PointCloudCharacteristics , typename CharacteristicFunction , typename DistanceBetweenPointsCharacteristics > +double bootstrap( size_t number_of_points , CharacteristicFunction f , DistanceBetweenPointsCharacteristics distance , size_t number_of_repetitions , size_t size_of_subsample , double quantile = 0.95 ) +{ + bool dbg = false; + + if ( size_of_subsample >= number_of_points ) + { + std::cerr << "Size of subsample is greater or equal to the number of points. The bootstrap procedure do not make sense in this case. \n"; + return 0; + } + + //initialization of a random number generator: + std::srand ( unsigned ( std::time(0) ) ); + + //we will shuffle the vector of numbers 0,1,2,...,points.size()-1 in order to pick a subset of a size size_of_subsample + std::vector< size_t > numbers_to_sample( number_of_points ); + for ( size_t i = 0 ; i != number_of_points ; ++i ) + { + numbers_to_sample[i] = i; + } + + //now we compute the characteristic od all the points: + PointCloudCharacteristics characteristic_of_all_points = f( numbers_to_sample ); + + //vector to keep the distances between characteristic_of_points and characteristic_of_subsample: + std::vector< double > vector_of_distances( number_of_repetitions , 0 ); + + + + #ifdef GUDHI_USE_TBB + tbb::parallel_for ( tbb::blocked_range<size_t>(0, number_of_repetitions), [&](const tbb::blocked_range<size_t>& range) + { + for ( size_t it_no = range.begin() ; it_no != range.end() ; ++it_no ) + #else + for ( size_t it_no = 0 ; it_no < number_of_repetitions ; ++it_no ) + #endif + { + if ( dbg ) + { + std::cout << "Still : " << number_of_repetitions-it_no << " tests to go. \n The subsampled vector consist of points number : "; + } + //do a random shuffle of vector_of_characteristics_of_poits + std::random_shuffle( numbers_to_sample.begin() , numbers_to_sample.end() ); + + //construct a vector< PointType > of a size size_of_subsample: + std::vector< size_t > subsampled_points; + subsampled_points.reserve( size_of_subsample ); + for ( size_t i = 0 ; i != size_of_subsample ; ++i ) + { + subsampled_points.push_back( numbers_to_sample[i] ); + if ( dbg )std::cout << numbers_to_sample[i] << " , "; + } + + + //now we can compute characteristic of subsampled_points: + PointCloudCharacteristics characteristic_of_subsampled_points = f( subsampled_points ); + if ( dbg )std::cout << std::endl << "Characteristic of subsampled points computed.\n"; + + //and now we compute distance between characteristic_of_points and characteristic_of_subsample. Note that subsampled points go first, and this is neded, since sometimes all points are not needed. + double dist = distance( characteristic_of_subsampled_points , characteristic_of_all_points ); + + if ( dbg ) std::cout << "The distance between characteristic of all points and the characteristic of subsample is : " << dist << std::endl; + + vector_of_distances[it_no] = dist; + } + #ifdef GUDHI_USE_TBB + } + ); + #endif + + size_t position_of_quantile = floor(quantile*vector_of_distances.size()); + if ( position_of_quantile ) --position_of_quantile; + if ( dbg ) + { + std::cout << "position_of_quantile : " << position_of_quantile << ", and here is the array : " << std::endl; + for ( size_t i = 0 ; i != vector_of_distances.size() ; ++i ) + { + std::cout << vector_of_distances[i] << std::endl; + } + std::cout << std::endl; + } + + //now we need to sort the vector_of_distances and find the quantile: + std::nth_element (vector_of_distances.begin(), vector_of_distances.begin()+position_of_quantile, vector_of_distances.end()); + + if ( dbg )std::cout << "Result : " << vector_of_distances[ position_of_quantile ] << std::endl; + + return vector_of_distances[ position_of_quantile ]; + +}//bootstrap + +#endif diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h index f846ac08..df7c97af 100644 --- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h +++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h @@ -240,7 +240,7 @@ public: * In the first line, the values min and max of the image are stored * In the next lines, we have the persistence images in a form of a bitmap image. **/ - void print_to_file( const char* filename ); + void print_to_file( const char* filename )const; /** * A function that load a heat map from file to the current object (and arase qhatever was stored in the current object before). @@ -251,7 +251,7 @@ public: /** * The procedure checks if min_, max_ and this->heat_maps sizes are the same. **/ - inline bool check_if_the_same( const Persistence_heat_maps& second ) + inline bool check_if_the_same( const Persistence_heat_maps& second )const { bool dbg = false; if ( this->heat_map.size() != second.heat_map.size() ) @@ -277,18 +277,18 @@ public: /** * Return minimal range value of persistent image. **/ - inline double gimme_min(){return this->min_;} + inline double gimme_min()const{return this->min_;} /** * Return maximal range value of persistent image. **/ - inline double gimme_max(){return this->max_;} + inline double gimme_max()const{return this->max_;} /** * Operator == to check if to persistence heat maps are the same. **/ - bool operator == ( const Persistence_heat_maps& rhs ) + bool operator == ( const Persistence_heat_maps& rhs )const { bool dbg = false; if ( !this->check_if_the_same(rhs) ) @@ -317,7 +317,7 @@ public: /** * Operator != to check if to persistence heat maps are different. **/ - bool operator != ( const Persistence_heat_maps& rhs ) + bool operator != ( const Persistence_heat_maps& rhs )const { return !( (*this) == rhs ); } @@ -326,7 +326,7 @@ public: /** * A function to generate a gnuplot script to vizualize the persistent image. **/ - void plot( const char* filename ); + void plot( const char* filename )const; //Implementations of functions for various concepts. @@ -334,11 +334,11 @@ public: /** * This function produce a vector of doubles based on a persisence heat map. It is required in a concept Vectorized_topological_data */ - std::vector<double> vectorize( int number_of_function ); + std::vector<double> vectorize( int number_of_function )const; /** * This function return the number of functions that allows vectorization of persistence heat map. It is required in a concept Vectorized_topological_data. **/ - size_t number_of_vectorize_functions() + size_t number_of_vectorize_functions()const { return this->number_of_functions_for_vectorization; } @@ -346,7 +346,7 @@ public: /** * This function is required by the Real_valued_topological_data concept. It returns various projections od the persistence heat map to a real line. **/ - double project_to_R( int number_of_function ); + double project_to_R( int number_of_function )const; /** * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept. **/ @@ -360,7 +360,7 @@ public: * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. * This function is required in Topological_data_with_distances concept. **/ - double distance( const Persistence_heat_maps& second_ , double power = 1); + double distance( const Persistence_heat_maps& second_ , double power = 1)const; /** * A function to compute averaged persistence heat map, based on vector of persistence heat maps. @@ -373,10 +373,30 @@ public: * The parameter of this functionis a const reference to an object of a class Persistence_heat_maps. * This function is required in Topological_data_with_scalar_product concept. **/ - double compute_scalar_product( const Persistence_heat_maps& second_ ); + double compute_scalar_product( const Persistence_heat_maps& second_ )const; //end of implementation of functions needed for concepts. + + /** + * The x-range of the persistence heat map. + **/ + std::pair< double , double > gimme_x_range()const + { + return std::make_pair( this->min_ , this->max_ ); + } + + /** + * The y-range of the persistence heat map. + **/ + std::pair< double , double > gimme_y_range()const + { + return this->gimme_x_range(); + } + + + + protected: //private methods std::vector< std::vector<double> > check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps ); @@ -644,11 +664,11 @@ void Persistence_heat_maps<Scalling_of_kernels>::compute_percentage_of_active( c template <typename Scalling_of_kernels> -void Persistence_heat_maps<Scalling_of_kernels>::plot( const char* filename ) +void Persistence_heat_maps<Scalling_of_kernels>::plot( const char* filename )const { ofstream out; std::stringstream ss; - ss << filename << "_Gnuplot_script"; + ss << filename << "_GnuplotScript"; out.open( ss.str().c_str() ); out << "plot '-' matrix with image" << std::endl; @@ -666,7 +686,7 @@ void Persistence_heat_maps<Scalling_of_kernels>::plot( const char* filename ) template <typename Scalling_of_kernels> -void Persistence_heat_maps<Scalling_of_kernels>::print_to_file( const char* filename ) +void Persistence_heat_maps<Scalling_of_kernels>::print_to_file( const char* filename )const { ofstream out; @@ -746,7 +766,7 @@ void Persistence_heat_maps<Scalling_of_kernels>::load_from_file( const char* fil //Concretizations of virtual methods: template <typename Scalling_of_kernels> -std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize( int number_of_function ) +std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize( int number_of_function )const { //convert this->heat_map into one large vector: size_t size_of_result = 0; @@ -770,7 +790,7 @@ std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize( int n } template <typename Scalling_of_kernels> -double Persistence_heat_maps<Scalling_of_kernels>::distance( const Persistence_heat_maps& second , double power ) +double Persistence_heat_maps<Scalling_of_kernels>::distance( const Persistence_heat_maps& second , double power )const { //first we need to check if (*this) and second are defined on the same domain and have the same dimensions: if ( !this->check_if_the_same(second) ) @@ -793,7 +813,7 @@ double Persistence_heat_maps<Scalling_of_kernels>::distance( const Persistence_h } template <typename Scalling_of_kernels> -double Persistence_heat_maps<Scalling_of_kernels>::project_to_R( int number_of_function ) +double Persistence_heat_maps<Scalling_of_kernels>::project_to_R( int number_of_function )const { double result = 0; for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) @@ -813,7 +833,7 @@ void Persistence_heat_maps<Scalling_of_kernels>::compute_average( const std::vec } template <typename Scalling_of_kernels> -double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product( const Persistence_heat_maps& second ) +double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product( const Persistence_heat_maps& second )const { //first we need to check if (*this) and second are defined on the same domain and have the same dimensions: if ( !this->check_if_the_same(second) ) diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h index 5d21ba97..b917de69 100644 --- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h +++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h @@ -61,55 +61,85 @@ public: /**
* The procedure returns a pair the first element of which is the leftmost end of the interval, and the second element of which is the rightmost end of the interval.
**/
- std::pair<double,double> min_max();
+ std::pair<double,double> min_max()const;
+
+ /**
+ * This procedure returns x-range of a given persistence diagram.
+ **/
+ std::pair< double , double > gimme_x_range()const
+ {
+ double min_ = std::numeric_limits<int>::max();
+ double max_ = -std::numeric_limits<int>::max();
+ for ( size_t i = 0 ; i != this->intervals.size() ; ++i )
+ {
+ if ( this->intervals[i].first < min_ )min_ = this->intervals[i].first;
+ if ( this->intervals[i].first > max_ )max_ = this->intervals[i].first;
+ }
+ return std::make_pair( min_ , max_ );
+ }
+
+ /**
+ * This procedure returns y-range of a given persistence diagram.
+ **/
+ std::pair< double , double > gimme_y_range()const
+ {
+ double min_ = std::numeric_limits<int>::max();
+ double max_ = -std::numeric_limits<int>::max();
+ for ( size_t i = 0 ; i != this->intervals.size() ; ++i )
+ {
+ if ( this->intervals[i].second < min_ )min_ = this->intervals[i].second;
+ if ( this->intervals[i].second > max_ )max_ = this->intervals[i].second;
+ }
+ return std::make_pair( min_ , max_ );
+ }
/**
* Procedure that compute the vector of lengths of the dominant (i.e. the longest) persistence intervals. The list is truncated at the parameter of the call where_to_cut (set by default to 100).
**/
- std::vector<double> length_of_dominant_intervals( size_t where_to_cut = 100 );
+ std::vector<double> length_of_dominant_intervals( size_t where_to_cut = 100 )const;
/**
* Procedure that compute the vector of the dominant (i.e. the longest) persistence intervals. The parameter of the procedure (set by default to 100) is the number of dominant intervals returned by the procedure.
**/
- std::vector< std::pair<double,double> > dominant_intervals( size_t where_to_cut = 100 );
+ std::vector< std::pair<double,double> > dominant_intervals( size_t where_to_cut = 100 )const;
/**
* Procedure to compute a histogram of interva's length. A histogram is a block plot. The number of blocks is determined by the first parameter of the function (set by default to 10).
* For the sake of argument let us assume that the length of the longest interval is 1 and the number of bins is 10. In this case the i-th block correspond to a range between i-1/10 and i10.
* The vale of a block supported at the interval is the number of persistence intervals of a length between x_0 and x_1.
**/
- std::vector< size_t > histograms_of_lengths( size_t number_of_bins = 10 );
+ std::vector< size_t > histograms_of_lengths( size_t number_of_bins = 10 )const;
/**
* Based on a histogram of intervals lengts computed by the function histograms_of_lengths H the procedure below computes the cumulative histogram. The i-th position of the resulting histogram
* is the sume of values of H for the positions from 0 to i.
**/
- std::vector< size_t > cumulative_histograms_of_lengths( size_t number_of_bins = 10 );
+ std::vector< size_t > cumulative_histograms_of_lengths( size_t number_of_bins = 10 )const;
/**
* In this procedure we assume that each barcode is a characteristic function of a hight equal to its length. The persistence diagram is a sum of such a functions. The procedure below construct a function being a
* sum of the characteristic functions of persistence intervals. The first two parameters are the range in which the function is to be computed and the last parameter is the number of bins in
* the discretization of the interval [_min,_max].
**/
- std::vector< double > characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins = 10 );
+ std::vector< double > characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins = 10 )const;
/**
* Cumulative version of the function characteristic_function_of_diagram
**/
- std::vector< double > cumulative_characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins = 10 );
+ std::vector< double > cumulative_characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins = 10 )const;
/**
* Compute the funtion of persistence Betti numbers. The returned value is a vector of pair. First element of each pair is a place where persistence Betti numbers change.
* Second element of each pair is the value of Persistence Betti numbers at that point.
**/
- std::vector< std::pair< double , size_t > > compute_persistent_betti_numbers();
+ std::vector< std::pair< double , size_t > > compute_persistent_betti_numbers()const;
/**
*This is a non optimal procedure that compute vector of distances from each point of diagram to its k-th nearest neighbor (k is a parameted of the program). The resulting vector is by default truncated to 10
*elements (this value can be changed by using the second parameter of the program). The points are returned in order from the ones which are farthest away from their k-th nearest neighbors.
**/
- std::vector< double > k_n_n( size_t k , size_t where_to_cut = 10 );
+ std::vector< double > k_n_n( size_t k , size_t where_to_cut = 10 )const;
/**
* Operator that send the diagram to a stream.
@@ -126,7 +156,7 @@ public: /**
* Generating gnuplot script to plot the interval.
**/
- void plot( const char* filename )
+ void plot( const char* filename , double min_x = -1 , double max_x = -1 , double min_y = -1 , double max_y = -1 ) const
{
//this program create a gnuplot script file that allows to plot persistence diagram.
std::ofstream out;
@@ -135,9 +165,18 @@ public: nameSS << filename << "_GnuplotScript";
std::string nameStr = nameSS.str();
out.open( (char*)nameStr.c_str() );
+
std::pair<double,double> min_max_values = this->min_max();
- out << "set xrange [" << min_max_values.first - 0.1*(min_max_values.second-min_max_values.first) << " : " << min_max_values.second + 0.1*(min_max_values.second-min_max_values.first) << " ]" << std::endl;
- out << "set yrange [" << min_max_values.first - 0.1*(min_max_values.second-min_max_values.first) << " : " << min_max_values.second + 0.1*(min_max_values.second-min_max_values.first) << " ]" << std::endl;
+ if ( min_x == max_x )
+ {
+ out << "set xrange [" << min_max_values.first - 0.1*(min_max_values.second-min_max_values.first) << " : " << min_max_values.second + 0.1*(min_max_values.second-min_max_values.first) << " ]" << std::endl;
+ out << "set yrange [" << min_max_values.first - 0.1*(min_max_values.second-min_max_values.first) << " : " << min_max_values.second + 0.1*(min_max_values.second-min_max_values.first) << " ]" << std::endl;
+ }
+ else
+ {
+ out << "set xrange [" << min_x << " : " << max_x << " ]" << std::endl;
+ out << "set yrange [" << min_y << " : " << max_y << " ]" << std::endl;
+ }
out << "plot '-' using 1:2 notitle \"" << filename << "\", \\" << std::endl;
out << " '-' using 1:2 notitle with lp" << std::endl;
for ( size_t i = 0 ; i != this->intervals.size() ; ++i )
@@ -152,6 +191,8 @@ public: std::cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << std::endl;
}
+
+
/**
@@ -162,7 +203,7 @@ public: /**
* Return the persistence interval at the given position. Note that intervals are not sorted with respect to their lengths.
**/
- inline std::pair< double,double > operator [](size_t i)
+ inline std::pair< double,double > operator [](size_t i)const
{
if ( i >= this->intervals.size() )throw("Index out of range! Operator [], one_d_gaussians class\n");
return this->intervals[i];
@@ -187,11 +228,11 @@ public: * This is a simple function projectig the persistence intervals to a real number. The function we use here is a sum of squared lendgths of intervals. It can be naturally interpreted as
* sum of step function, where the step hight it equal to the length of the interval.
**/
- double project_to_R( int number_of_function );
+ double project_to_R( int number_of_function )const;
/**
* The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept.
**/
- size_t number_of_projections_to_R()
+ size_t number_of_projections_to_R()const
{
return this->number_of_functions_for_projections_to_reals;
}
@@ -199,14 +240,14 @@ public: /**
* Return a familly of vectors obtained from the persistence diagram. The i-th vector consist of the lenfth of i dominant persistence intervals.
**/
- std::vector<double> vectorize( int number_of_function )
+ std::vector<double> vectorize( int number_of_function )const
{
return this->length_of_dominant_intervals( number_of_function );
}
/**
* This function return the number of functions that allows vectorization of a persisence diagram. It is required in a concept Vectorized_topological_data.
**/
- size_t number_of_vectorize_functions()
+ size_t number_of_vectorize_functions()const
{
return this->number_of_functions_for_vectorization;
}
@@ -215,7 +256,7 @@ public: *Computations of distance from the current persistnce diagram to the persistence diagram given as a parameter of this function.
*The last parameter, power, is here in case we would like to compute p=th Wasserstein distance. At the moment, for the bottleneck distances, it will be ignored.
**/
- double distance( const Persistence_intervals& second , double power = 1)
+ double distance( const Persistence_intervals& second , double power = 1) const
{
return 1;
//waiting for Francois Godi for the code. We will compute here the Bottleneck distnace.
@@ -306,9 +347,8 @@ Persistence_intervals::Persistence_intervals( const std::vector< std::pair< doub }
-std::vector< double > Persistence_intervals::length_of_dominant_intervals( size_t where_to_cut )
-{
- this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+std::vector< double > Persistence_intervals::length_of_dominant_intervals( size_t where_to_cut )const
+{
std::vector< double > result( this->intervals.size() );
for ( size_t i = 0 ; i != this->intervals.size() ; ++i )
{
@@ -329,7 +369,7 @@ bool compare( const std::pair< size_t , double >& first , const std::pair< size_ }
-std::vector< std::pair<double,double> > Persistence_intervals::dominant_intervals( size_t where_to_cut )
+std::vector< std::pair<double,double> > Persistence_intervals::dominant_intervals( size_t where_to_cut )const
{
bool dbg = false;
std::vector< std::pair< size_t , double > > position_length_vector( this->intervals.size() );
@@ -353,7 +393,7 @@ std::vector< std::pair<double,double> > Persistence_intervals::dominant_interval }//dominant_intervals
-std::vector< size_t > Persistence_intervals::histograms_of_lengths( size_t number_of_bins )
+std::vector< size_t > Persistence_intervals::histograms_of_lengths( size_t number_of_bins )const
{
bool dbg = false;
@@ -401,7 +441,7 @@ std::vector< size_t > Persistence_intervals::histograms_of_lengths( size_t numbe }
-std::vector< size_t > Persistence_intervals::cumulative_histograms_of_lengths( size_t number_of_bins )
+std::vector< size_t > Persistence_intervals::cumulative_histograms_of_lengths( size_t number_of_bins )const
{
std::vector< size_t > histogram = this->histograms_of_lengths( number_of_bins );
std::vector< size_t > result( histogram.size() );
@@ -416,7 +456,7 @@ std::vector< size_t > Persistence_intervals::cumulative_histograms_of_lengths( s }
-std::vector< double > Persistence_intervals::characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins )
+std::vector< double > Persistence_intervals::characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins )const
{
bool dbg = false;
@@ -476,7 +516,7 @@ std::vector< double > Persistence_intervals::characteristic_function_of_diagram( -std::vector< double > Persistence_intervals::cumulative_characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins )
+std::vector< double > Persistence_intervals::cumulative_characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins )const
{
std::vector< double > intsOfBars = this->characteristic_function_of_diagram( x_min , x_max , number_of_bins );
std::vector< double > result( intsOfBars.size() );
@@ -497,7 +537,7 @@ bool compare_first_element_of_pair( const std::pair< T , bool >& f, const std::p }
-std::vector< std::pair< double , size_t > > Persistence_intervals::compute_persistent_betti_numbers()
+std::vector< std::pair< double , size_t > > Persistence_intervals::compute_persistent_betti_numbers()const
{
std::vector< std::pair< double , bool > > places_where_pbs_change( 2*this->intervals.size() );
@@ -536,7 +576,7 @@ inline double compute_euclidean_distance( const std::pair< double,double > & f , }
-std::vector< double > Persistence_intervals::k_n_n( size_t k , size_t where_to_cut )
+std::vector< double > Persistence_intervals::k_n_n( size_t k , size_t where_to_cut )const
{
bool dbg = false;
if ( dbg )
@@ -650,7 +690,7 @@ std::vector< double > Persistence_intervals::k_n_n( size_t k , size_t where_to_c }
-std::pair<double,double> Persistence_intervals::min_max()
+std::pair<double,double> Persistence_intervals::min_max()const
{
double min_ = std::numeric_limits<int>::max();
double max_ = -std::numeric_limits<int>::max();
@@ -662,7 +702,7 @@ std::pair<double,double> Persistence_intervals::min_max() return std::make_pair( min_ , max_ );
}
-double Persistence_intervals::project_to_R( int number_of_function )
+double Persistence_intervals::project_to_R( int number_of_function )const
{
double result = 0;
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h index b3caaedc..fd6837c0 100644 --- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h +++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h @@ -257,6 +257,24 @@ public: } return maxValue; } + + + /** + * Computations of minimum (y) value of landscape. + **/ + double compute_minimum()const + { + double minValue = 0; + if ( this->land.size() ) + { + minValue = std::numeric_limits<int>::max(); + for ( size_t i = 0 ; i != this->land[0].size() ; ++i ) + { + if ( this->land[0][i].second < minValue )minValue = this->land[0][i].second; + } + } + return minValue; + } /** * Computations of a L^i norm of landscape, where i is the input parameter. @@ -334,7 +352,7 @@ public: * The number of projections to R is defined to the number of nonzero landscape functions. I-th projection is an integral of i-th landscape function over whole R. * This function is required by the Real_valued_topological_data concept. **/ - double project_to_R( int number_of_function ) + double project_to_R( int number_of_function )const { return this->compute_integral_of_landscape( (size_t)number_of_function ); } @@ -342,7 +360,7 @@ public: /** * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept. **/ - size_t number_of_projections_to_R() + size_t number_of_projections_to_R()const { return this->number_of_functions_for_projections_to_reals; } @@ -350,7 +368,7 @@ public: /** * This function produce a vector of doubles based on a landscape. It is required in a concept Vectorized_topological_data */ - std::vector<double> vectorize( int number_of_function ) + std::vector<double> vectorize( int number_of_function )const { //TODO, think of something smarter over here std::vector<double> v; @@ -368,7 +386,7 @@ public: /** * This function return the number of functions that allows vectorization of persistence laandscape. It is required in a concept Vectorized_topological_data. **/ - size_t number_of_vectorize_functions() + size_t number_of_vectorize_functions()const { return this->number_of_functions_for_vectorization; } @@ -433,7 +451,7 @@ public: * The parameter of this functionis a Persistence_landscape. * This function is required in Topological_data_with_distances concept. **/ - double distance( const Persistence_landscape& second , double power = 1 ) + double distance( const Persistence_landscape& second , double power = 1 )const { if ( power != -1 ) { @@ -451,14 +469,50 @@ public: * The parameter of this functionis a Persistence_landscape. * This function is required in Topological_data_with_scalar_product concept. **/ - double compute_scalar_product( const Persistence_landscape& second ) + double compute_scalar_product( const Persistence_landscape& second )const { return compute_inner_product( (*this) , second ); } //end of implementation of functions needed for concepts. - + /** + * This procedure returns x-range of a given level persistence landscape. If a default value is used, the x-range + * of 0th level landscape is given (and this range contains the ranges of all other landscapes). + **/ + std::pair< double , double > gimme_x_range( size_t level = 0 )const + { + std::pair< double , double > result; + if ( level < this->land.size() ) + { + result = std::make_pair( this->land[level][1].first , this->land[level][ this->land[level].size() - 2 ].first ); + } + else + { + result = std::make_pair( 0,0 ); + } + return result; + } + + /** + * This procedure returns y-range of a given level persistence landscape. If a default value is used, the y-range + * of 0th level landscape is given (and this range contains the ranges of all other landscapes). + **/ + std::pair< double , double > gimme_y_range( size_t level = 0 )const + { + std::pair< double , double > result; + if ( level < this->land.size() ) + { + double maxx = this->compute_maximum(); + double minn = this->compute_minimum(); + result = std::make_pair( minn , maxx ); + } + else + { + result = std::make_pair( 0,0 ); + } + return result; + } @@ -484,7 +538,7 @@ public: //a function used to create a gnuplot script for visualization of landscapes - void plot( const char* filename ,int from = -1, int to = -1 , double xRangeBegin = -1 , double xRangeEnd = -1 , double yRangeBegin = -1 , double yRangeEnd = -1 ); + void plot( const char* filename, double xRangeBegin = -1 , double xRangeEnd = -1 , double yRangeBegin = -1 , double yRangeEnd = -1,int from = -1, int to = -1 ); protected: @@ -1436,7 +1490,7 @@ double compute_inner_product( const Persistence_landscape& l1 , const Persistenc } -void Persistence_landscape::plot( const char* filename , int from, int to , double xRangeBegin , double xRangeEnd , double yRangeBegin , double yRangeEnd ) +void Persistence_landscape::plot( const char* filename, double xRangeBegin , double xRangeEnd , double yRangeBegin , double yRangeEnd , int from , int to ) { //this program create a gnuplot script file that allows to plot persistence diagram. ofstream out; diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h index 80c6628d..6ac4deca 100644 --- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h +++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h @@ -512,6 +512,56 @@ public: } /** + * This procedure returns x-range of a given level persistence landscape. If a default value is used, the x-range + * of 0th level landscape is given (and this range contains the ranges of all other landscapes). + **/ + std::pair< double , double > gimme_x_range( size_t level = 0 )const + { + return std::make_pair( this->grid_min , this->grid_max ); + //std::pair< double , double > result; + //if ( level < this->land.size() ) + //{ + // double dx = (this->grid_max - this->grid_min)/(double)this->values_of_landscapes.size(); + // size_t first_nonzero = 0; + // while ( (first_nonzero != this->values_of_landscapes.size()) && (this->values_of_landscapes[level][first_nonzero] == 0) )++first_nonzero; + // + // if ( first_nonzero == 0 ) + // { + // return std::make_pair( 0,0 );//this landscape is empty. + // } + // + // size_t last_nonzero = 0; + // while ( (last_nonzero != 0) && (this->values_of_landscapes[level][last_nonzero] == 0) )--last_nonzero; + // + // result = std::make_pair( this->grid_min +first_nonzero*dx , this->grid_max - last_nonzero*dx ); + //} + //else + //{ + // result = std::make_pair( 0,0 ); + //} + //return result; + } + + /** + * This procedure returns y-range of a persistence landscape. If a default value is used, the y-range + * of 0th level landscape is given (and this range contains the ranges of all other landscapes). + **/ + std::pair< double , double > gimme_y_range( size_t level = 0 )const + { + return this->compute_minimum_maximum(); + //std::pair< double , double > result; + //if ( level < this->land.size() ) + //{ + // result = this->compute_minimum_maximum() + //} + //else + //{ + // result = std::make_pair( 0,0 ); + //} + //return result; + } + + /** * This function computes maximal lambda for which lambda-level landscape is nonzero. **/ size_t number_of_nonzero_levels()const @@ -527,7 +577,7 @@ public: /** * Computations of a L^i norm of landscape, where i is the input parameter. **/ - double compute_norm_of_landscape( double i ) + double compute_norm_of_landscape( double i )const { std::vector< std::pair< double , double > > p; Persistence_landscape_on_grid l(p,this->grid_min,this->grid_max,this->values_of_landscapes.size()-1); @@ -788,7 +838,7 @@ public: * The number of projections to R is defined to the number of nonzero landscape functions. I-th projection is an integral of i-th landscape function over whole R. * This function is required by the Real_valued_topological_data concept. **/ - double project_to_R( int number_of_function ) + double project_to_R( int number_of_function )const { return this->compute_integral_of_landscape( (size_t)number_of_function ); } @@ -796,7 +846,7 @@ public: /** * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept. **/ - size_t number_of_projections_to_R() + size_t number_of_projections_to_R()const { return number_of_functions_for_projections_to_reals; } @@ -807,7 +857,7 @@ public: /** * This function produce a vector of doubles based on a landscape. It is required in a concept Vectorized_topological_data */ - std::vector<double> vectorize( int number_of_function ) + std::vector<double> vectorize( int number_of_function )const { //TODO, think of something smarter over here if ( ( number_of_function < 0 ) || ( (size_t)number_of_function >= this->values_of_landscapes.size() ) ) @@ -821,7 +871,7 @@ public: /** * This function return the number of functions that allows vectorization of persistence laandscape. It is required in a concept Vectorized_topological_data. **/ - size_t number_of_vectorize_functions() + size_t number_of_vectorize_functions()const { return number_of_functions_for_vectorization; } @@ -902,7 +952,7 @@ public: * The parameter of this functionis a Persistence_landscape_on_grid. * This function is required in Topological_data_with_distances concept. **/ - double distance( const Persistence_landscape_on_grid& second , double power = 1 ) + double distance( const Persistence_landscape_on_grid& second , double power = 1 )const { if ( power != -1 ) { @@ -948,15 +998,23 @@ public: /** * A function that returns values of landsapes. It can be used for vizualization **/ - std::vector< std::vector< double > > output_for_visualization() + std::vector< std::vector< double > > output_for_visualization()const { return this->values_of_landscapes; } /** - * function used to create a gnuplot script for visualization of landscapes + * function used to create a gnuplot script for visualization of landscapes. Over here we need to specify which landscapes do we want to plot. **/ - void plot( const char* filename , size_t from_ = std::numeric_limits<size_t>::max(), size_t to_ = std::numeric_limits<size_t>::max() ); + void plot( const char* filename , size_t from_ , size_t to_ )const + { + this->plot( filename , from_ , to_ ); + } + + /** + * function used to create a gnuplot script for visualization of landscapes. Over here we can restrict also x and y range of the landscape. + **/ + void plot( const char* filename, double min_x = -1 , double max_x = -1 , double min_y = -1 , double max_y = -1 , size_t from_ = std::numeric_limits<size_t>::max(), size_t to_= std::numeric_limits<size_t>::max() )const; protected: @@ -1158,7 +1216,7 @@ void Persistence_landscape_on_grid::print_to_file( const char* filename )const out.close(); } -void Persistence_landscape_on_grid::plot( const char* filename , size_t from_ , size_t to_ ) +void Persistence_landscape_on_grid::plot( const char* filename, double min_x , double max_x , double min_y , double max_y, size_t from_ , size_t to_ )const { //this program create a gnuplot script file that allows to plot persistence diagram. ofstream out; @@ -1168,9 +1226,17 @@ void Persistence_landscape_on_grid::plot( const char* filename , size_t from_ , std::string nameStr = nameSS.str(); out.open( (char*)nameStr.c_str() ); - std::pair<double,double> min_max = compute_minimum_maximum(); - out << "set xrange [" << this->grid_min << " : " << this->grid_max << "]" << endl; - out << "set yrange [" << min_max.first << " : " << min_max.second << "]" << endl; + if ( min_x == max_x ) + { + std::pair<double,double> min_max = compute_minimum_maximum(); + out << "set xrange [" << this->grid_min << " : " << this->grid_max << "]" << endl; + out << "set yrange [" << min_max.first << " : " << min_max.second << "]" << endl; + } + else + { + out << "set xrange [" << min_x << " : " << max_x << "]" << endl; + out << "set yrange [" << min_y << " : " << max_y << "]" << endl; + } size_t number_of_nonzero_levels = this->number_of_nonzero_levels(); double dx = ( this->grid_max - this->grid_min )/((double)this->values_of_landscapes.size()-1); diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h b/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h index 7211f6eb..9d02e917 100644 --- a/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h +++ b/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h @@ -33,6 +33,7 @@ //gudhi include #include <gudhi/concretizations/read_persitence_from_file.h> #include <gudhi/common_gudhi_stat.h> + using namespace std; namespace Gudhi @@ -42,16 +43,31 @@ namespace Gudhi_stat template <typename T> -struct euclidean_distance +struct Euclidean_distance { double operator() ( const std::pair< T,T >& f , const std::pair<T,T>& s ) { return sqrt( (f.first-s.first)*(f.first-s.first) + (f.second-s.second)*(f.second-s.second) ); } + double operator() ( const std::vector< T >& f , const std::vector < T >& s ) + { + if ( f.size() != s.size() ) + { + std::cerr << "Not compatible points dimensions in the procedure to compute Euclidean distance. The program will now terminate. \n"; + std::cout << f.size() << " , " << s.size() << std::endl; + throw "Not compatible points dimensions in the procedure to compute Euclidean distance. The program will now terminate. \n"; + } + double result = 0; + for ( size_t i = 0 ; i != f.size() ; ++i ) + { + result += ( f[i]-s[i] )*( f[i]-s[i] ); + } + return sqrt( result ); + } }; template <typename T> -struct maximum_distance +struct Maximum_distance { double operator() ( const std::pair< T,T >& f , const std::pair<T,T>& s ) { @@ -113,7 +129,7 @@ public: /** * This procedure gives the value of a vector on a given position. **/ - inline double vector_in_position( size_t position ) + inline double vector_in_position( size_t position )const { if ( position >= this->sorted_vector_of_distances.size() )throw("Wrong position in accessing Vector_distances_in_diagram::sorted_vector_of_distances\n"); return this->sorted_vector_of_distances[position]; @@ -127,12 +143,12 @@ public: /** * Write a vector to a file. **/ - void write_to_file( const char* filename ); + void write_to_file( const char* filename )const; /** * Write a vector to a file. **/ - void print_to_file( const char* filename ) + void print_to_file( const char* filename )const { this->write_to_file( filename ); } @@ -145,7 +161,7 @@ public: /** * Comparision operators: **/ - bool operator == ( const Vector_distances_in_diagram& second ) + bool operator == ( const Vector_distances_in_diagram& second )const { if ( this->sorted_vector_of_distances.size() != second.sorted_vector_of_distances.size() )return false; for ( size_t i = 0 ; i != this->sorted_vector_of_distances.size() ; ++i ) @@ -155,7 +171,7 @@ public: return true; } - bool operator != ( const Vector_distances_in_diagram& second ) + bool operator != ( const Vector_distances_in_diagram& second )const { return !( *this == second ); } @@ -167,11 +183,11 @@ public: /** * Compute projection to real numbers of persistence vector. This function is required by the Real_valued_topological_data concept **/ - double project_to_R( int number_of_function ); + double project_to_R( int number_of_function )const; /** * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept. **/ - size_t number_of_projections_to_R() + size_t number_of_projections_to_R()const { return this->number_of_functions_for_projections_to_reals; } @@ -179,11 +195,11 @@ public: /** * Compute a vectorization of a persistent vectors. It is required in a concept Vectorized_topological_data. **/ - std::vector<double> vectorize( int number_of_function ); + std::vector<double> vectorize( int number_of_function )const; /** * This function return the number of functions that allows vectorization of a persisence vector. It is required in a concept Vectorized_topological_data. **/ - size_t number_of_vectorize_functions() + size_t number_of_vectorize_functions()const { return this->number_of_functions_for_vectorization; } @@ -196,12 +212,12 @@ public: /** * Compute a distance of two persistent vectors. This function is required in Topological_data_with_distances concept. **/ - double distance( const Vector_distances_in_diagram& second , double power = 1); + double distance( const Vector_distances_in_diagram& second , double power = 1)const; /** * Compute a scalar product of two persistent vectors. This function is required in Topological_data_with_scalar_product concept. **/ - double compute_scalar_product( const Vector_distances_in_diagram& second ); + double compute_scalar_product( const Vector_distances_in_diagram& second )const; //end of implementation of functions needed for concepts. @@ -224,7 +240,7 @@ public: /** * For visualization use output from vectorize and build histograms. **/ - std::vector< double > output_for_visualization() + std::vector< double > output_for_visualization()const { return this->sorted_vector_of_distances; } @@ -233,10 +249,10 @@ public: /** * Create a gnuplot script to vizualize the data structure. **/ - void plot( const char* filename ) + void plot( const char* filename )const { std::stringstream gnuplot_script; - gnuplot_script << filename << "_Gnuplot_script"; + gnuplot_script << filename << "_GnuplotScript"; ofstream out; out.open( gnuplot_script.str().c_str() ); out << "set style data histogram" << std::endl; @@ -251,6 +267,23 @@ public: out.close(); std::cout << "To vizualize, open gnuplot and type: load \'" << gnuplot_script.str().c_str() << "\'" << std::endl; } + + /** + * The x-range of the persistence vector. + **/ + std::pair< double , double > gimme_x_range()const + { + return std::make_pair( 0 , this->sorted_vector_of_distances.size() ); + } + + /** + * The y-range of the persistence vector. + **/ + std::pair< double , double > gimme_y_range()const + { + if ( this->sorted_vector_of_distances.size() == 0 )return std::make_pair(0,0); + return std::make_pair( this->sorted_vector_of_distances[0] , 0); + } private: @@ -458,7 +491,7 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distances_via_vect //Implementations of functions for various concepts. template <typename F> -double Vector_distances_in_diagram<F>::project_to_R( int number_of_function ) +double Vector_distances_in_diagram<F>::project_to_R( int number_of_function )const { if ( number_of_function > this->number_of_functions_for_projections_to_reals )throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::project_to_R"; if ( number_of_function < 0 )throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::project_to_R"; @@ -508,7 +541,7 @@ void Vector_distances_in_diagram<F>::compute_average( const std::vector< Vector_ } template <typename F> -double Vector_distances_in_diagram<F>::distance( const Vector_distances_in_diagram& second_ , double power ) +double Vector_distances_in_diagram<F>::distance( const Vector_distances_in_diagram& second_ , double power )const { bool dbg = false; @@ -569,7 +602,7 @@ double Vector_distances_in_diagram<F>::distance( const Vector_distances_in_diagr } template < typename F> -std::vector<double> Vector_distances_in_diagram<F>::vectorize( int number_of_function ) +std::vector<double> Vector_distances_in_diagram<F>::vectorize( int number_of_function )const { if ( number_of_function > this->number_of_functions_for_vectorization )throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::vectorize"; if ( number_of_function < 0 )throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::vectorize"; @@ -584,7 +617,7 @@ std::vector<double> Vector_distances_in_diagram<F>::vectorize( int number_of_fun template < typename F> -void Vector_distances_in_diagram<F>::write_to_file( const char* filename ) +void Vector_distances_in_diagram<F>::write_to_file( const char* filename )const { std::ofstream out; out.open( filename ); @@ -620,7 +653,7 @@ void Vector_distances_in_diagram<F>::load_from_file( const char* filename ) } template < typename F> -double Vector_distances_in_diagram<F>::compute_scalar_product( const Vector_distances_in_diagram& second_vector ) +double Vector_distances_in_diagram<F>::compute_scalar_product( const Vector_distances_in_diagram& second_vector )const { double result = 0; for ( size_t i = 0 ; i != std::min(this->sorted_vector_of_distances.size(),second_vector.sorted_vector_of_distances.size()) ; ++i ) diff --git a/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h b/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h index 953547c7..e16c1f20 100644 --- a/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h +++ b/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h @@ -27,6 +27,8 @@ #include <iostream> #include <fstream> #include <sstream> +#include <vector> +#include <algorithm> #include <unistd.h> @@ -252,6 +254,43 @@ std::vector< std::pair< double , double > > read_gudhi_file( const char* filenam return barcode; }//read_gudhi_file + +std::vector< std::vector< double > > read_numbers_from_file_line_by_line( const char* filename ) +{ + bool dbg = false; + if ( !( access( filename, F_OK ) != -1 ) ) + { + std::cerr << "The file : " << filename << " do not exist. The program will now terminate \n"; + throw "The file from which you are trying to read the persistence landscape do not exist. The program will now terminate \n"; + } + + std::vector< std::vector< double > > result; + double number; + + std::ifstream in(filename); + std::string line; + while ( in.good() ) + { + std::getline(in,line); + std::stringstream ss; + ss << line; + + if ( dbg )std::cerr << "\n Reading line : " << line << std::endl; + + std::vector< double > this_line; + while ( ss.good() ) + { + ss >> number; + this_line.push_back( number ); + if ( dbg )std::cerr << number << " "; + } + if ( this_line.size() && in.good() ) result.push_back( this_line ); + } + in.close(); + + return result; +}//read_numbers_from_file_line_by_line + }//namespace Gudhi_stat }//namespace Gudhi diff --git a/src/Gudhi_stat/include/gudhi/multiplicative_bootstrap.h b/src/Gudhi_stat/include/gudhi/multiplicative_bootstrap.h new file mode 100644 index 00000000..fc7a8091 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/multiplicative_bootstrap.h @@ -0,0 +1,140 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef BOOTSTRAP_H +#define BOOTSTRAP_H + +//concretizations +#include <gudhi/concretizations/Vector_distances_in_diagram.h> +#include <gudhi/concretizations/Persistence_landscape.h> +#include <gudhi/concretizations/Persistence_landscape_on_grid.h> +#include <gudhi/concretizations/Persistence_heat_maps.h> + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_sort.h> +#endif + +#include <vector> +#include <algorithm> +#include <omp.h> +#include <random> +#include <ctime> + +template < typename TopologicalObject > +class difference_of_objects +{ +public: + TopologicalObject operator()( const TopologicalObject& first, const TopologicalObject& second )const + { + return first-second; + } +}; + +template < typename TopologicalObject > +class norm_of_objects +{ +public: + double operator()( const TopologicalObject& obj )const + { + TopologicalObject* empty = new TopologicalObject; + double dist = empty->distance( obj ); + delete empty; + return dist; + } +}; + + + +/** +* This is a generic function to perform multiplicative bootstrap. +**/ + + +template < typename TopologicalObject , typename OperationOnObjects , typename NormOnObjects > +double multiplicative_bootstrap( const std::vector< TopologicalObject* >& topological_objects , size_t number_of_bootstrap_operations , const OperationOnObjects& oper , const NormOnObjects& norm , double quantile = 0.95 ) +{ + bool dbg = false; + + //initialization of a random number generator: + std::random_device rd; + std::mt19937 generator( time(NULL) ); + std::normal_distribution<> norm_distribution(0.,1.); + + + //first compute an average of topological_objects + TopologicalObject average; + average.compute_average( topological_objects ); + + std::vector< double > vector_of_intermediate_characteristics( number_of_bootstrap_operations , 0 ); + + + #ifdef GUDHI_USE_TBB + tbb::parallel_for ( tbb::blocked_range<size_t>(0, number_of_bootstrap_operations), [&](const tbb::blocked_range<size_t>& range) + { + for ( size_t it_no = range.begin() ; it_no != range.end() ; ++it_no ) + #else + for ( size_t it_no = 0 ; it_no < number_of_bootstrap_operations ; ++it_no ) + #endif + { + if ( dbg ) + { + std::cout << "Still : " << number_of_bootstrap_operations-it_no << " tests to go. \n The subsampled vector consist of points number : "; + } + //and compute the intermediate characteristic: + TopologicalObject result; + for ( size_t i = 0 ; i != topological_objects.size() ; ++i ) + { + double rand_variable = norm_distribution( generator ); + result += rand_variable*oper(*(topological_objects[i]) , average); + } + result = result*(1.0/sqrt( topological_objects.size() )); + vector_of_intermediate_characteristics[it_no] = norm(result); + } + #ifdef GUDHI_USE_TBB + } + ); + #endif + + + + size_t position_of_quantile = floor(quantile*vector_of_intermediate_characteristics.size()); + if ( position_of_quantile ) --position_of_quantile; + if ( dbg ) + { + std::cout << "position_of_quantile : " << position_of_quantile << ", and here is the array : " << std::endl; + for ( size_t i = 0 ; i != vector_of_intermediate_characteristics.size() ; ++i ) + { + std::cout << vector_of_intermediate_characteristics[i] << std::endl; + } + std::cout << std::endl; + } + + //now we need to sort the vector_of_distances and find the quantile: + std::nth_element (vector_of_intermediate_characteristics.begin(), vector_of_intermediate_characteristics.begin()+position_of_quantile, vector_of_intermediate_characteristics.end()); + double result = vector_of_intermediate_characteristics[ position_of_quantile ]/(sqrt( topological_objects.size() )); + if ( dbg )std::cout << "Result : " << result << std::endl; + + return result; + +}//bootstrap + +#endif diff --git a/src/Gudhi_stat/include/gudhi/permutation_test.h b/src/Gudhi_stat/include/gudhi/permutation_test.h index 18b94eab..ac3a7f53 100644 --- a/src/Gudhi_stat/include/gudhi/permutation_test.h +++ b/src/Gudhi_stat/include/gudhi/permutation_test.h @@ -41,13 +41,16 @@ namespace Gudhi_stat template <typename Representation_of_persistence> double permutation_test( const std::vector<Representation_of_persistence*>& data_1 , const std::vector<Representation_of_persistence*>& data_2 , size_t number_of_permutations , double exponent = 1 ) { + bool dbg = true; try { - Representation_of_persistence* av_1 = new Representation_of_persistence; - av_1->compute_average( data_1 ); - Representation_of_persistence* av_2 = new Representation_of_persistence; - av_2->compute_average( data_2 ); - double initial_distance = av_1->distance( *av_2 , exponent ); + Representation_of_persistence av_1;// = new Representation_of_persistence; + av_1.compute_average( data_1 ); + Representation_of_persistence av_2;// = new Representation_of_persistence; + av_2.compute_average( data_2 ); + double initial_distance = av_1.distance( av_2 , exponent ); + + double counter = 0; for ( size_t i = 0 ; i != number_of_permutations ; ++i ) @@ -60,10 +63,22 @@ double permutation_test( const std::vector<Representation_of_persistence*>& data std::vector<Representation_of_persistence*> first_part(&all_data[0],&all_data[data_1.size()]); std::vector<Representation_of_persistence*> second_part(&all_data[data_1.size()],&all_data[all_data.size()]); - av_1->compute_average( first_part ); - av_2->compute_average( second_part ); - double distance_after_permutations = av_1->distance( *av_2 , exponent ); + + Representation_of_persistence av_1;// = new Representation_of_persistence; + Representation_of_persistence av_2;// = new Representation_of_persistence; + av_1.compute_average( first_part ); + av_2.compute_average( second_part ); + double distance_after_permutations = av_1.distance( av_2 , exponent ); + //delete av_1; + //delete av_2; if ( distance_after_permutations > initial_distance )++counter; + if ( dbg ) + { + std::cerr << "Permutation number : " << i << std::endl; + std::cerr << "distance_after_permutations : " << distance_after_permutations << std::endl; + std::cerr << "initial_distance : " << initial_distance << std::endl; + std::cerr << "counter : " << counter << std::endl; + } } return counter / (double)number_of_permutations; } diff --git a/src/Gudhi_stat/include/gudhi/topological_process.h b/src/Gudhi_stat/include/gudhi/topological_process.h index 729ef2f8..e7bd46de 100644 --- a/src/Gudhi_stat/include/gudhi/topological_process.h +++ b/src/Gudhi_stat/include/gudhi/topological_process.h @@ -28,6 +28,7 @@ #include <gudhi/concretizations/Vector_distances_in_diagram.h> #include <gudhi/concretizations/Persistence_landscape.h> #include <gudhi/concretizations/Persistence_landscape_on_grid.h> +#include <gudhi/concretizations/Persistence_heat_maps.h> #include <vector> //extras @@ -38,11 +39,18 @@ namespace Gudhi namespace Gudhi_stat { +//over here we will need a few version of construct_representation_from_file procedure, since different representations may require different parameters. This is a procedure that in my +//oppinion cannot be standarize, since construction of representation cannot. But, the remaining part of the code in my opinion is free from any details of representation. -template <typename Representation> -std::vector< Representation* > construct_representation_from_file( const char* filename ) + +//the reason I am separating the process of getting the intervals from the process of consstructing the represnetation is the following: for soem representations, we may need to ahve the same +//scale. To determine the scale, we need to know all the intervals before. So, we read the intervals first, then if needed, process them, to get the parametersof the representation, +//and only then construct the representations. +std::vector< std::vector< std::pair< double , double > > > read_persistence_pairs_from_file_with_names_of_files( const char* filename ) { bool dbg = false; + std::vector< std::vector< std::pair< double , double > > > result; + std::vector< std::string > files = readFileNames( filename ); std::cout << "Here are the filenames in the file : " << filename << std::endl; @@ -51,11 +59,10 @@ std::vector< Representation* > construct_representation_from_file( const char* f std::cout << files[i] << std::endl; } - std::vector< Representation* > result( files.size() ); for ( size_t i = 0 ; i != files.size() ; ++i ) { std::vector< std::pair< double , double > > diag = read_standard_file( files[i].c_str() ); - + result.push_back( diag ); if ( dbg ) { std::cerr << "Here is a diagram from a file : " << files[i].c_str() << std::endl; @@ -65,8 +72,134 @@ std::vector< Representation* > construct_representation_from_file( const char* f } getchar(); } - - Representation* l = new Representation( diag ); + } + return result; +} + +//When workign with time varying data, we tpically have a collection of files with intervals to produce one process. The intervals from all the files that constitute a single process can be read with +//the read_persistence_pairs_from_file_with_names_of_files procedure. +//But then, we also may need to read the persistence intervals of collection of proesses we may want to read the intervals first. This is what the procedre +std::vector< std::vector< std::vector< std::pair< double , double > > > > read_persistence_pairs_from_files_with_names_of_files( const std::vector<const char*>& filenames ) +{ + std::vector< std::vector< std::vector< std::pair< double , double > > > > result; + result.reserve( filenames.size() ); + for ( size_t file_no = 0 ; file_no != filenames.size() ; ++file_no ) + { + result.push_back( read_persistence_pairs_from_file_with_names_of_files( filenames[file_no] ) ); + } + return result; +}//read_persistence_pairs_from_files_with_names_of_files + + +std::pair< std::pair< double,double > , std::pair< double,double > > find_x_and_y_ranges_of_intervals( const std::vector< std::vector< std::pair< double , double > > >& intervals_from_file ) +{ + double min_x = std::numeric_limits< double >::max(); + double max_x = -std::numeric_limits< double >::max(); + double min_y = std::numeric_limits< double >::max(); + double max_y = -std::numeric_limits< double >::max(); + for ( size_t i = 0 ; i != intervals_from_file.size() ; ++i ) + { + for ( size_t j = 0 ; j != intervals_from_file[i].size() ; ++j ) + { + if ( min_x > intervals_from_file[i][j].first )min_x = intervals_from_file[i][j].first; + if ( max_x < intervals_from_file[i][j].first )max_x = intervals_from_file[i][j].first; + + if ( min_y > intervals_from_file[i][j].second )min_y = intervals_from_file[i][j].second; + if ( max_y < intervals_from_file[i][j].second )max_y = intervals_from_file[i][j].second; + } + } + return std::make_pair( std::make_pair( min_x,max_x ) , std::make_pair( min_y,max_y ) ); +}//find_x_and_y_ranges_of_intervals + + +std::pair< std::pair< double,double > , std::pair< double,double > > find_x_and_y_ranges_of_intervals( const std::vector< std::vector< std::vector< std::pair< double , double > > > >& intervals_from_files ) +{ + double min_x = std::numeric_limits< double >::max(); + double max_x = -std::numeric_limits< double >::max(); + double min_y = std::numeric_limits< double >::max(); + double max_y = -std::numeric_limits< double >::max(); + for ( size_t i = 0 ; i != intervals_from_files.size() ; ++i ) + { + std::pair< std::pair< double,double > , std::pair< double,double > > ranges = find_x_and_y_ranges_of_intervals( intervals_from_files[i] ); + if ( min_x > ranges.first.first ) min_x = ranges.first.first; + if ( max_x < ranges.first.second ) max_x = ranges.first.second; + if ( min_y > ranges.second.first ) min_y = ranges.second.first; + if ( max_y < ranges.second.second ) max_y = ranges.second.second; + } + return std::make_pair( std::make_pair( min_x,max_x ) , std::make_pair( min_y,max_y ) ); +} + + +template <typename Representation> +std::vector< Representation* > construct_representation_from_file( const char* filename ) +{ + std::vector< std::vector< std::pair< double , double > > > intervals_from_file = read_persistence_pairs_from_file_with_names_of_files( filename ); + std::vector< Representation* > result( intervals_from_file.size() ); + for ( size_t i = 0 ; i != intervals_from_file.size() ; ++i ) + { + Representation* l = new Representation( intervals_from_file[i] ); + result[i] = l; + } + return result; +} + +//this one can be use for Persistence_intervals.h and Persistence_landscape.h +template <typename Representation> +std::vector< Representation* > construct_representation_from_file( const std::vector< std::vector< std::pair< double , double > > >& intervals_from_file) +{ + + std::vector< Representation* > result( intervals_from_file.size() ); + for ( size_t i = 0 ; i != intervals_from_file.size() ; ++i ) + { + Representation* l = new Representation( intervals_from_file[i] ); + result[i] = l; + } + return result; +} + +//this one can be use for Persistence_heat_maps.h +template <typename Representation> +std::vector< Representation* > construct_representation_from_file( const std::vector< std::vector< std::pair< double , double > > >& intervals_from_file , + std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1), + bool erase_below_diagonal = false , + size_t number_of_pixels = 1000 , + double min_ = -1 , double max_ = -1 ) +{ + std::vector< Representation* > result( intervals_from_file.size() ); + for ( size_t i = 0 ; i != intervals_from_file.size() ; ++i ) + { + Representation* l = new Representation( intervals_from_file[i] , filter , erase_below_diagonal , number_of_pixels , min_ , max_ ); + result[i] = l; + } + return result; +} + +//this one can be use for Persistence_landscape_on_grid.h +template <typename Representation> +std::vector< Representation* > construct_representation_from_file( const std::vector< std::vector< std::pair< double , double > > >& intervals_from_file, + double grid_min_ , double grid_max_ , size_t number_of_points_ +) +{ + std::vector< Representation* > result( intervals_from_file.size() ); + for ( size_t i = 0 ; i != intervals_from_file.size() ; ++i ) + { + Representation* l = new Representation( intervals_from_file[i] , grid_min_ , grid_max_ , number_of_points_ ); + result[i] = l; + } + return result; +} + + +//this one can be use for Vector_distances_in_diagram.h +template <typename Representation> +std::vector< Representation* > construct_representation_from_file( const std::vector< std::vector< std::pair< double , double > > >& intervals_from_file, + size_t where_to_cut +) +{ + std::vector< Representation* > result( intervals_from_file.size() ); + for ( size_t i = 0 ; i != intervals_from_file.size() ; ++i ) + { + Representation* l = new Representation( intervals_from_file[i] , where_to_cut ); result[i] = l; } return result; @@ -77,9 +210,33 @@ template <typename Representation> class Topological_process { public: - Topological_process(); + Topological_process(){}; + ~Topological_process() + { + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + delete this->data[i]; + } + } + Topological_process( const Topological_process& org ) + { + this->data = std::vector< Representation* >( org.data.size() ); + for ( size_t i = 0 ; i != org.data.size() ; ++i ) + { + this->data[i] = new Representation( *org.data[i] ); + } + } + Topological_process& operator = ( const Topological_process& rhs ) + { + for ( size_t i = 0 ; i != rhs.data.size() ; ++i ) + { + this->data[i] = new Representation( *rhs.data[i] ); + } + return *this; + } + Topological_process( const std::vector< Representation* >& data_ ):data(data_){} - double distance( const Topological_process& second ) + double distance( const Topological_process& second , double exponent = 1 ) { if ( this->data.size() != second.data.size() ) { @@ -88,12 +245,12 @@ public: double result = 0; for ( size_t i = 0 ; i != this->data.size() ; ++i ) { - result += this->data[i]->distance( *second.data[i] ); + result += this->data[i]->distance( *second.data[i] , exponent ); } return result; } - void compute_average( const std::vector< Representation* >& to_average ) + void compute_average( const std::vector< Topological_process* >& to_average ) { //since we will substitute whatever data we have in this object with an average, we clear the data in this object first: this->data.clear(); @@ -101,14 +258,127 @@ public: if ( to_average.size() == 0 )return; for ( size_t i = 1 ; i != to_average.size() ; ++i ) { - if ( to_average[0].data.size() != to_average[i].data.size ) + if ( to_average[0]->data.size() != to_average[i]->data.size() ) { throw "Incompatible lengths of topological processes, the averages cannot be computed \n"; } } + + this->data.reserve( to_average[0]->data.size() ); + + for ( size_t level = 0 ; level != to_average[0]->data.size() ; ++level ) + { + //now we will be averaging the level level: + std::vector< Representation* > to_average_at_this_level; + to_average_at_this_level.reserve( to_average.size() ); + for ( size_t i = 0 ; i != to_average.size() ; ++i ) + { + to_average_at_this_level.push_back( to_average[i]->data[level] ); + } + Representation* average_representation_on_this_level = new Representation; + average_representation_on_this_level->compute_average( to_average_at_this_level ); + this->data.push_back( average_representation_on_this_level ); + } + } + + std::pair< double , double > gimme_x_range()const + { + double min_x = std::numeric_limits< double >::max(); + double max_x = -std::numeric_limits< double >::max(); + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + std::pair< double , double > xrange = this->data[i]->gimme_x_range(); + if ( min_x > xrange.first )min_x = xrange.first; + if ( max_x < xrange.second )max_x = xrange.second; + } + return std::make_pair( min_x , max_x ); + } + + std::pair< double , double > gimme_y_range()const + { + double min_y = std::numeric_limits< double >::max(); + double max_y = -std::numeric_limits< double >::max(); + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + std::pair< double , double > yrange = this->data[i]->gimme_y_range(); + if ( min_y > yrange.first )min_y = yrange.first; + if ( max_y < yrange.second )max_y = yrange.second; + } + return std::make_pair( min_y , max_y ); + } + + /** + * The procedure checks if the ranges of data are the same for all of them. + **/ + bool are_the_data_aligned()const + { + if ( this->data.size() == 0 )return true;//empty collection is aligned + std::pair< double , double > x_range = this->data[0]->gimme_x_range(); + std::pair< double , double > y_range = this->data[0]->gimme_y_range(); + for ( size_t i = 1 ; i != this->data.size() ; ++i ) + { + if ( (x_range != this->data[i]->gimme_x_range()) || (y_range != this->data[i]->gimme_y_range()) ) + { + return false; + } + } + return true; } + + //scalar products? //confidence bounds? + + void plot( const char* filename , size_t delay = 30 , double min_x = -1 , double max_x = -1 , double min_y = -1 , double max_y = -1 ) + { + std::vector< std::string > filenames; + //over here we need to + for ( size_t i = 0 ; i != this->data.size() ; ++i ) + { + std::stringstream ss; + ss << filename << "_" << i; + if ( ( min_x != max_x ) && ( min_y != max_y ) ) + { + //in this case, we set up the uniform min and max values for pciture + //this->data[i]->plot( ss.str().c_str() , min_x , max_x , min_y , max_y ); + } + else + { + //in this case, we DO NOT set up the uniform min and max values for pciture + this->data[i]->plot( ss.str().c_str() ); + } + ss << "_GnuplotScript"; + filenames.push_back( ss.str() ); + } + //and now we have to call it somehow for all the files. We will create a script that will call all the other scripts and ceate a sequence of jpg/png files. + + + std::stringstream gif_file_name; + gif_file_name << filename << ".gif"; + std::stringstream gif_gnuplot_script_file_name; + gif_gnuplot_script_file_name << filename << "_gif_gnuplot_script"; + + ofstream out; + out.open( gif_gnuplot_script_file_name.str().c_str() ); + out << "set terminal gif animate delay " << delay << std::endl; + out << "set output '" << gif_file_name.str() << "'" << std::endl; + if ( min_x != max_x ) + { + out << "set xrange [" << min_x << ":" << max_x << "]" << std::endl; + out << "set yrange [" << min_y << ":" << max_y << "]" << std::endl; + } + + for ( size_t i = 0 ; i != filenames.size() ; ++i ) + { + out << " load '" << filenames[i] << "'" << std::endl; + } + out.close(); + + std::cout << std::endl << std::endl << std::endl << "Open gnuplot terminal and type load '" << gif_gnuplot_script_file_name.str() << "' to create an animated gif \n"; + }//plot + + + std::vector< Representation* > gimme_data(){return this->data;} private: std::vector< Representation* > data; }; diff --git a/src/Gudhi_stat/test/vector_representation_test.cpp b/src/Gudhi_stat/test/vector_representation_test.cpp index 3e14bd59..2131ed29 100644 --- a/src/Gudhi_stat/test/vector_representation_test.cpp +++ b/src/Gudhi_stat/test/vector_representation_test.cpp @@ -51,10 +51,10 @@ BOOST_AUTO_TEST_CASE(check_read_write_to_files) intervals.push_back( std::make_pair(4,7) ); intervals.push_back( std::make_pair(9,10) ); intervals.push_back( std::make_pair(3,11) ); - Vector_distances_in_diagram< euclidean_distance<double> > p( intervals , -1 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p( intervals , -1 ); p.write_to_file( "test_vector_representation_write_read" ); - Vector_distances_in_diagram< euclidean_distance<double> > q; + Vector_distances_in_diagram< Euclidean_distance<double> > q; q.load_from_file( "test_vector_representation_write_read" ); BOOST_CHECK( p == q ); @@ -64,7 +64,7 @@ BOOST_AUTO_TEST_CASE(check_read_write_to_files) BOOST_AUTO_TEST_CASE(check_sortev_vector_distances_template) { - Vector_distances_in_diagram< euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); std::vector< double > sortev_vector_distances_template; sortev_vector_distances_template.push_back( 0.609968 ); @@ -184,7 +184,7 @@ BOOST_AUTO_TEST_CASE(check_sortev_vector_distances_template) BOOST_AUTO_TEST_CASE(check_projections_to_R) { - Vector_distances_in_diagram< euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); std::vector< double > proj; proj.push_back( 0 ); proj.push_back( 0.6099679993 ); @@ -300,8 +300,8 @@ BOOST_AUTO_TEST_CASE(check_projections_to_R) BOOST_AUTO_TEST_CASE(check_distance_computations) { - Vector_distances_in_diagram< euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); - Vector_distances_in_diagram< euclidean_distance<double> > p_prime( "data/file_with_diagram" , 10 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p_prime( "data/file_with_diagram" , 10 ); std::vector< std::pair<double,double> > intervals(10); intervals[0] = std::pair<double,double>( 1,2 ); intervals[1] = std::pair<double,double>( 3,4 ); @@ -313,7 +313,7 @@ BOOST_AUTO_TEST_CASE(check_distance_computations) intervals[7] = std::pair<double,double>( 15,16 ); intervals[8] = std::pair<double,double>( 17,18 ); intervals[9] = std::pair<double,double>( 19,20 ); - Vector_distances_in_diagram< euclidean_distance<double> > p_bis( intervals , 10 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p_bis( intervals , 10 ); //cerr << "p_prime.distance( (Abs_Topological_data_with_distances*)(&p_bis) , 1 ) : " << p_prime.distance( (Abs_Topological_data_with_distances*)(&p_bis) , 1 ) << endl; BOOST_CHECK( almost_equal ( p_prime.distance( p_bis , 1 ) , 1.86428 ) ); } @@ -322,7 +322,7 @@ BOOST_AUTO_TEST_CASE(check_distance_computations) BOOST_AUTO_TEST_CASE(check_compute_average) { - Vector_distances_in_diagram< euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); + Vector_distances_in_diagram< Euclidean_distance<double> > p( "data/file_with_diagram" , 100 ); //compute average std::vector< std::pair<double,double> > i1(3); i1[0] = std::pair<double,double>( 1,2 ); @@ -334,17 +334,17 @@ BOOST_AUTO_TEST_CASE(check_compute_average) i2[1] = std::pair<double,double>( 2,15); i2[2] = std::pair<double,double>( 6,17 ); - Vector_distances_in_diagram< euclidean_distance<double> > A( i1 , -1 ); - Vector_distances_in_diagram< euclidean_distance<double> > B( i1 , -1 ); + Vector_distances_in_diagram< Euclidean_distance<double> > A( i1 , -1 ); + Vector_distances_in_diagram< Euclidean_distance<double> > B( i1 , -1 ); - std::vector< Vector_distances_in_diagram< euclidean_distance<double> >* > to_average; + std::vector< Vector_distances_in_diagram< Euclidean_distance<double> >* > to_average; to_average.push_back( &A ); to_average.push_back( &B ); - Vector_distances_in_diagram< euclidean_distance<double> > average; + Vector_distances_in_diagram< Euclidean_distance<double> > average; average.compute_average( to_average ); - Vector_distances_in_diagram< euclidean_distance<double> > template_average; + Vector_distances_in_diagram< Euclidean_distance<double> > template_average; template_average.load_from_file( "data/average_of_persistence_vectors" ); BOOST_CHECK( template_average == average ); diff --git a/src/Gudhi_stat/utilities/CMakeLists.txt b/src/Gudhi_stat/utilities/CMakeLists.txt index f98bcecb..a26aa5c5 100644 --- a/src/Gudhi_stat/utilities/CMakeLists.txt +++ b/src/Gudhi_stat/utilities/CMakeLists.txt @@ -111,6 +111,9 @@ add_executable ( topological_process_2 topological_process_2.cpp ) target_link_libraries(topological_process_2 ${Boost_SYSTEM_LIBRARY}) add_executable ( Hausdorff_bootstrap Hausdorff_bootstrap.cpp ) +if (TBB_FOUND) +target_link_libraries(Hausdorff_bootstrap ${TBB_LIBRARIES}) +endif(TBB_FOUND) target_link_libraries(Hausdorff_bootstrap ${Boost_SYSTEM_LIBRARY}) @@ -120,3 +123,10 @@ target_link_libraries(Landscape_bootstrap ${TBB_LIBRARIES}) endif(TBB_FOUND) target_link_libraries(Landscape_bootstrap ${Boost_SYSTEM_LIBRARY}) + +add_executable ( Multiplicative_bootstrap Multiplicative_bootstrap.cpp ) +if (TBB_FOUND) +target_link_libraries(Multiplicative_bootstrap ${TBB_LIBRARIES}) +endif(TBB_FOUND) +target_link_libraries(Multiplicative_bootstrap ${Boost_SYSTEM_LIBRARY}) + diff --git a/src/Gudhi_stat/utilities/Hausdorff_bootstrap.cpp b/src/Gudhi_stat/utilities/Hausdorff_bootstrap.cpp index 95e1aa99..16c9b085 100644 --- a/src/Gudhi_stat/utilities/Hausdorff_bootstrap.cpp +++ b/src/Gudhi_stat/utilities/Hausdorff_bootstrap.cpp @@ -38,7 +38,7 @@ int main( int argc , char** argv ) std::cout << "(a) a name of a file with points," << std:: endl; std::cout << "(b) a number of repetitions of bootstrap (integer)," << std::endl; std::cout << "(c) a size of subsample (integer, smaller than the number of points," << std::endl; - std::cout << "(d) a squantile (real number between 0 and 1. If you do not know what to set, set it to 0.95." << std::endl; + 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 != 5 ) { std::cerr << "Wrong number of parameters, the program will now terminate.\n"; diff --git a/src/Gudhi_stat/utilities/Landscape_bootstrap.cpp b/src/Gudhi_stat/utilities/Landscape_bootstrap.cpp index e9e39e50..5b3d0609 100644 --- a/src/Gudhi_stat/utilities/Landscape_bootstrap.cpp +++ b/src/Gudhi_stat/utilities/Landscape_bootstrap.cpp @@ -143,7 +143,7 @@ int main( int argc , char** argv ) 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 squantile (real number between 0 and 1. If you do not know what to set, set it to 0.95." << std::endl; + 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"; diff --git a/src/Gudhi_stat/utilities/Multiplicative_bootstrap.cpp b/src/Gudhi_stat/utilities/Multiplicative_bootstrap.cpp new file mode 100644 index 00000000..489ad909 --- /dev/null +++ b/src/Gudhi_stat/utilities/Multiplicative_bootstrap.cpp @@ -0,0 +1,69 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + + +#include <gudhi/Hausdorff_distances.h> +#include <gudhi/multiplicative_bootstrap.h> +#include <gudhi/concretizations/read_persitence_from_file.h> +#include <gudhi/concretizations/Vector_distances_in_diagram.h> + +using namespace std; +using namespace Gudhi; +using namespace Gudhi::Gudhi_stat; + + + +int main( int argc , char** argv ) +{ + std::cout << "The parameters of this program are : " << std::endl; + std::cout << "(a) a name of a file with names of files with points," << std:: endl; + std::cout << "(b) a number of repetitions of bootstrap (integer)," << std::endl; + std::cout << "(c) 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 != 4 ) + { + std::cerr << "Wrong number of parameters, the program will now terminate.\n"; + return 1; + } + + const char* file_with_filenames = argv[1]; + size_t number_of_repetitions_of_bootstrap = (size_t)atoi( argv[2] ); + double quantile = atof( argv[3] ); + + std::vector< std::string > filenames = readFileNames( file_with_filenames ); + std::vector< Persistence_landscape* > collection_of_landscapes( filenames.size() ); + for ( size_t i = 0 ; i != filenames.size() ; ++i ) + { + std::vector< std::pair< double , double > > diag = read_standard_file( filenames[i].c_str() ); + collection_of_landscapes[i] = new Persistence_landscape( diag ); + } + + //now we can run the bootstrap: + difference_of_objects<Persistence_landscape> diff; + norm_of_objects<Persistence_landscape> norm; + + double result = + multiplicative_bootstrap< Persistence_landscape , difference_of_objects<Persistence_landscape> , norm_of_objects<Persistence_landscape> >( collection_of_landscapes , number_of_repetitions_of_bootstrap , diff , norm , quantile ); + + std::cout << "result of bootstrap : " << result << std::endl; + return 0; +} + |