diff options
Diffstat (limited to 'src/Gudhi_stat/include')
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/PSSK.h | 181 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/Persistence_heat_maps.h | 1042 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/Persistence_intervals.h | 700 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/Persistence_intervals_with_distances.h | 65 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/Persistence_landscape.h | 1498 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/Persistence_landscape_on_grid.h | 1563 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/Persistence_vectors.h | 778 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/common_persistence_representations.h | 154 | ||||
| -rw-r--r-- | src/Gudhi_stat/include/gudhi/read_persistence_from_file.h | 300 |
9 files changed, 6281 insertions, 0 deletions
diff --git a/src/Gudhi_stat/include/gudhi/PSSK.h b/src/Gudhi_stat/include/gudhi/PSSK.h new file mode 100644 index 00000000..bd142c02 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/PSSK.h @@ -0,0 +1,181 @@ +/* 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/>. + */ + +#pragma once +#ifndef PSSK_H +#define PSSK_H + +//gudhi include +#include <gudhi/Persistence_heat_maps.h> + + +namespace Gudhi +{ +namespace Persistence_representations +{ + +/** +* This is a version of a representation presented in https://arxiv.org/abs/1412.6821 +* In that paper the authors are using the representation just to compute kernel. Over here, we extend the usability by far. +* Note that the version presented here is not exact, since we are discretizing the kernel. +* The only difference with respect to the original class is the method of creation. We have full (square) image, and for every point (p,q), we add a kernel at (p,q) and the negative kernel +* at (q,p) +**/ + +class PSSK : public Persistence_heat_maps<constant_scaling_function> +{ +public: + PSSK():Persistence_heat_maps(){} + + PSSK(const std::vector< std::pair< double,double > > & interval , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1 ) + :Persistence_heat_maps() + { + this->construct( interval , filter , number_of_pixels , min_ , max_ ); + } + + + PSSK( const char* filename , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1 , unsigned dimension = std::numeric_limits<unsigned>::max() ): + Persistence_heat_maps() + { + std::vector< std::pair< double , double > > intervals_; + if ( dimension == std::numeric_limits<unsigned>::max() ) + { + intervals_ = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + intervals_ = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + this->construct( intervals_ , filter , number_of_pixels , min_ , max_ ); + } + +protected: + void construct( const std::vector< std::pair<double,double> >& intervals_ , + std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1), + size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1 ); +}; + +//if min_ == max_, then the program is requested to set up the values itself based on persistence intervals +void PSSK::construct( const std::vector< std::pair<double,double> >& intervals_ , + std::vector< std::vector<double> > filter, + size_t number_of_pixels , double min_ , double max_ ) +{ + bool dbg = false; + if ( dbg ){std::cerr << "Entering construct procedure \n";getchar();} + + if ( min_ == max_ ) + { + //in this case, we want the program to set up the min_ and max_ values by itself. + min_ = std::numeric_limits<int>::max(); + max_ = -std::numeric_limits<int>::max(); + + + for ( size_t i = 0 ; i != intervals_.size() ; ++i ) + { + if ( intervals_[i].first < min_ )min_ = intervals_[i].first; + if ( intervals_[i].second > max_ )max_ = intervals_[i].second; + } + //now we have the structure filled in, and moreover we know min_ and max_ values of the interval, so we know the range. + + //add some more space: + min_ -= fabs(max_ - min_)/100; + max_ += fabs(max_ - min_)/100; + } + + if ( dbg ) + { + std::cerr << "min_ : " << min_ << std::endl; + std::cerr << "max_ : " << max_ << std::endl; + std::cerr << "number_of_pixels : " << number_of_pixels << std::endl; + getchar(); + } + + this->min_ = min_; + this->max_ = max_; + + + + //initialization of the structure heat_map + std::vector< std::vector<double> > heat_map_; + for ( size_t i = 0 ; i != number_of_pixels ; ++i ) + { + std::vector<double> v( number_of_pixels , 0 ); + heat_map_.push_back( v ); + } + this->heat_map = heat_map_; + + if (dbg)std::cerr << "Done creating of the heat map, now we will fill in the structure \n"; + + for ( size_t pt_nr = 0 ; pt_nr != intervals_.size() ; ++pt_nr ) + { + //compute the value of intervals_[pt_nr] in the grid: + int x_grid = (int)((intervals_[pt_nr].first - this->min_)/( this->max_-this->min_ )*number_of_pixels); + int y_grid = (int)((intervals_[pt_nr].second - this->min_)/( this->max_-this->min_ )*number_of_pixels); + + if ( dbg ) + { + std::cerr << "point : " << intervals_[pt_nr].first << " , " << intervals_[pt_nr].second << std::endl; + std::cerr << "x_grid : " << x_grid << std::endl; + std::cerr << "y_grid : " << y_grid << std::endl; + } + + //x_grid and y_grid gives a center of the kernel. We want to have its lower left cordner. To get this, we need to shift x_grid and y_grid by a grid diameter. + x_grid -= filter.size()/2; + y_grid -= filter.size()/2; + //note that the numbers x_grid and y_grid may be negative. + + if ( dbg ) + { + std::cerr << "After shift : \n";; + std::cerr << "x_grid : " << x_grid << std::endl; + std::cerr << "y_grid : " << y_grid << std::endl; + std::cerr << "filter.size() : " << filter.size() << std::endl; + getchar(); + } + + + for ( size_t i = 0 ; i != filter.size() ; ++i ) + { + for ( size_t j = 0 ; j != filter.size() ; ++j ) + { + //if the point (x_grid+i,y_grid+j) is the correct point in the grid. + if ( + ((x_grid+i)>=0) && (x_grid+i<this->heat_map.size()) + && + ((y_grid+j)>=0) && (y_grid+j<this->heat_map.size()) + ) + { + if ( dbg ){std::cerr << y_grid+j << " " << x_grid+i << std::endl;} + this->heat_map[ y_grid+j ][ x_grid+i ] += filter[i][j]; + this->heat_map[ x_grid+i ][ y_grid+j ] += -filter[i][j]; + } + } + } + + } +}//construct + + +#endif + +}//namespace Gudhi_stat +}//namespace Gudhi diff --git a/src/Gudhi_stat/include/gudhi/Persistence_heat_maps.h b/src/Gudhi_stat/include/gudhi/Persistence_heat_maps.h new file mode 100644 index 00000000..59e58e41 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Persistence_heat_maps.h @@ -0,0 +1,1042 @@ +/* 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 Persistence_heat_maps_H +#define Persistence_heat_maps_H + +//standard include +#include <vector> +#include <sstream> +#include <iostream> +#include <cmath> +#include <limits> +#include <vector> +#include <algorithm> + +//gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> + + + + +namespace Gudhi +{ +namespace Persistence_representations +{ + + +/** + * This is a simple procedure to create n by n (or 2*pixel_radius times 2*pixel_radius cubical approximation of a Gaussian kernel. +**/ +std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius , double sigma ) +{ + bool dbg = false; + //we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a size 2*pixel_radius times 2*pixel_radius. + + double r = 0; + double sigma_sqr = sigma * sigma; + + // sum is for normalization + double sum = 0; + + //initialization of a kernel: + std::vector< std::vector<double> > kernel( 2*pixel_radius +1 ); + for ( size_t i = 0 ; i != kernel.size() ; ++i ) + { + std::vector<double> v( 2*pixel_radius +1 , 0 ); + kernel[i] = v; + } + + if ( dbg ) + { + std::cerr << "Kernel initalize \n"; + std::cerr << "pixel_radius : " << pixel_radius << std::endl; + std::cerr << "kernel.size() : " << kernel.size() << std::endl; + getchar(); + } + + for (int x = -pixel_radius; x <= (int)pixel_radius; x++) + { + for(int y = -pixel_radius; y <= (int)pixel_radius; y++) + { + double real_x = 2*sigma*x/pixel_radius; + double real_y = 2*sigma*y/pixel_radius; + r = sqrt(real_x*real_x + real_y*real_y); + kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r*r)/sigma_sqr))/(3.141592 * sigma_sqr); + sum += kernel[x + pixel_radius][y + pixel_radius]; + } + } + + // normalize the kernel + for( size_t i = 0; i != kernel.size() ; ++i) + { + for( size_t j = 0; j != kernel[i].size() ; ++j) + { + kernel[i][j] /= sum; + } + + } + + if ( dbg ) + { + std::cerr << "Here is the kernel : \n"; + for( size_t i = 0; i != kernel.size() ; ++i) + { + for( size_t j = 0; j != kernel[i].size() ; ++j) + { + std::cerr << kernel[i][j] << " "; + } + std::cerr << std::endl; + } + } + return kernel; +} + + +/* +* There are various options to scale the poits depending on their location. One can for instance: +* (1) do nothing (scale all of them with the weight 1), as in the function constant_function +* (2) Scale them by the distance to the diagonal. This is implemented in function +* (3) Scale them with the square of their distance to diagonal. This is implemented in function +* (4) Scale them with +*/ + + +/** + * This is one of a scaling functions used to weight poits depending on their persistence and/or location in the diagram. + * This particular functiona is a finction which always assign value 1 to a point in the diagram. +**/ +class constant_scaling_function +{ +public: + double operator()( const std::pair< double , double >& point_in_diagram ) + { + return 1; + } +}; + + +/** + * This is one of a scaling functions used to weight poits depending on their persistence and/or location in the diagram. + * The scaling given by this function to a point (b,d) is Euclidean distance of (b,d) from diagonal. +**/ +class distance_from_diagonal_scaling +{ +public: + double operator()( const std::pair< double , double >& point_in_diagram ) + { + //(point_in_diagram.first+point_in_diagram.second)/2.0 + return sqrt( pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2) ); + } +}; + +/** + * This is one of a scaling functions used to weight poits depending on their persistence and/or location in the diagram. + * The scaling given by this function to a point (b,d) is a square of Euclidean distance of (b,d) from diagonal. +**/ +class squared_distance_from_diagonal_scaling +{ +public: + double operator()( const std::pair< double , double >& point_in_diagram ) + { + return pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2); + } +}; + +/** + * This is one of a scaling functions used to weight poits depending on their persistence and/or location in the diagram. + * The scaling given by this function to a point (b,d) is an arctan of a persistence of a point (i.e. arctan( b-d ). +**/ +class arc_tan_of_persistence_of_point +{ +public: + double operator()( const std::pair< double , double >& point_in_diagram ) + { + return atan( point_in_diagram.second - point_in_diagram.first ); + } +}; + +/** + * This is one of a scaling functions used to weight poits depending on their persistence and/or location in the diagram. + * This scaling function do not only depend on a point (p,d) in the diagram, but it depends on the whole diagram. + * The longest persistence pair get a scaling 1. Any other pair get a scaling belong to [0,1], which is proportional + * to the persistence of that pair. +**/ +class weight_by_setting_maximal_interval_to_have_length_one +{ +public: + weight_by_setting_maximal_interval_to_have_length_one( double len ):letngth_of_maximal_interval(len){} + double operator()( const std::pair< double , double >& point_in_diagram ) + { + return (point_in_diagram.second-point_in_diagram.first)/this->letngth_of_maximal_interval; + } +private: + double letngth_of_maximal_interval; +}; + + +/** + * This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product +**/ +template <typename Scalling_of_kernels = constant_scaling_function> +class Persistence_heat_maps +{ +public: + /** + * The default constructor. A scaling function from the diagonal is set up to a constant function. The image is not erased below the diagonal. The gaussian have diameter 5. + **/ + Persistence_heat_maps() + { + Scalling_of_kernels f; + this->f = f; + this->erase_below_diagonal = false; + this->min_ = this->max_ = 0; + this->set_up_parameters_for_basic_classes(); + }; + + /** + * Construction that takes at the input the following parameters: + * (1) A vector of pairs of doubles (representing persistence intervals). All other parameters are optional. They are: + * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), + * (3) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). + * (4) a number of pixels in each direction (set to 1000 by default). + * (5) a min x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data, + * (6) a max x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data. + **/ + Persistence_heat_maps( const std::vector< std::pair< double,double > > & interval , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = std::numeric_limits<double>::max() , double max_ = std::numeric_limits<double>::max() ); + + /** + * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by default), a scaling function (constant by default), a boolean value which determines if the area of image below diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each direction (set to 1000 by default). and min and max values of images (both set to std::numeric_limits<double>::max() by defaulet. If this is the case, the program will pick the right values based on the data). + **/ + /** + * Construction that takes at the input the following parameters: + * (1) A a name of a file with persistence intervals. The file shold be readable by the function read_persistence_intervals_in_one_dimension_from_file. All other parameters are optional. They are: + * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), + * (3) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). + * (4) a number of pixels in each direction (set to 1000 by default). + * (5) a min x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data, + * (6) a max x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data. + **/ + Persistence_heat_maps( const char* filename , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1), bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = std::numeric_limits<double>::max() , double max_ = std::numeric_limits<double>::max() , unsigned dimension = std::numeric_limits<unsigned>::max() ); + + + /** + * Compute a mean value of a collection of heat maps and store it in the current object. Note that all the persistence maps send in a vector to this procedure need to have the same parameters. + * If this is not the case, the program will throw an exception. + **/ + void compute_mean( const std::vector<Persistence_heat_maps*>& maps ); + + /** + * Compute a median value of a collection of heat maps and store it in the current object. Note that all the persistence maps send in a vector to this procedure need to have the same parameters. + * If this is not the case, the program will throw an exception. + **/ + void compute_median( const std::vector<Persistence_heat_maps*>& maps ); + + /** + * Compute a percentage of active (i.e) values above the cutoff of a collection of heat maps. + **/ + void compute_percentage_of_active( const std::vector<Persistence_heat_maps*>& maps , size_t cutoff = 1 ); + + //put to file subroutine + /** + * The function outputs the perssitence image to a text file. The format as follow: + * 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 )const; + + /** + * A function that load a heat map from file to the current object (and arase qhatever was stored in the current object before). + **/ + void load_from_file( const char* filename ); + + + /** + * 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 )const + { + bool dbg = false; + if ( this->heat_map.size() != second.heat_map.size() ) + { + if ( dbg )std::cerr << "this->heat_map.size() : " << this->heat_map.size() << " \n second.heat_map.size() : " << second.heat_map.size() << std::endl; + return false; + } + if ( this->min_ != second.min_ ) + { + if ( dbg )std::cerr << "this->min_ : " << this->min_ << ", second.min_ : " << second.min_ << std::endl; + return false; + } + if ( this->max_ != second.max_ ) + { + if ( dbg )std::cerr << "this->max_ : " << this->max_ << ", second.max_ : " << second.max_ << std::endl; + return false; + } + //in the other case we may assume that the persistence images are defined on the same domain. + return true; + } + + + /** + * Return minimal range value of persistent image. + **/ + inline double get_min()const{return this->min_;} + + /** + * Return maximal range value of persistent image. + **/ + inline double get_max()const{return this->max_;} + + + /** + * Operator == to check if to persistence heat maps are the same. + **/ + bool operator == ( const Persistence_heat_maps& rhs )const + { + bool dbg = false; + if ( !this->check_if_the_same(rhs) ) + { + if ( dbg )std::cerr << "The domains are not the same \n"; + return false;//in this case, the domains are not the same, so the maps cannot be the same. + } + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + if ( !almost_equal(this->heat_map[i][j] , rhs.heat_map[i][j]) ) + { + if ( dbg ) + { + std::cerr << "this->heat_map[" << i << "][" << j << "] = " << this->heat_map[i][j] << std::endl; + std::cerr << "rhs.heat_map[" << i << "][" << j << "] = " << rhs.heat_map[i][j] << std::endl; + } + return false; + } + } + } + return true; + } + + /** + * Operator != to check if to persistence heat maps are different. + **/ + bool operator != ( const Persistence_heat_maps& rhs )const + { + return !( (*this) == rhs ); + } + + + /** + * A function to generate a gnuplot script to vizualize the persistent image. + **/ + void plot( const char* filename )const; + + + template<typename Operation_type> + friend Persistence_heat_maps operation_on_pair_of_heat_maps( const Persistence_heat_maps& first , const Persistence_heat_maps& second , Operation_type operation ) + { + //first check if the heat maps are compatible + if ( !first.check_if_the_same( second ) ) + { + std::cerr << "Sizes of the heat maps are not compatible. The program will now terminate \n"; + throw "Sizes of the heat maps are not compatible. The program will now terminate \n"; + } + Persistence_heat_maps result; + result.min_ = first.min_; + result.max_ = first.max_; + result.heat_map.reserve( first.heat_map.size() ); + for ( size_t i = 0 ; i != first.heat_map.size() ; ++i ) + { + std::vector< double > v; + v.reserve( first.heat_map[i].size() ); + for ( size_t j = 0 ; j != first.heat_map[i].size() ; ++j ) + { + v.push_back( operation( first.heat_map[i][j] , second.heat_map[i][j] ) ); + } + result.heat_map.push_back( v ); + } + return result; + }//operation_on_pair_of_heat_maps + + + /** + * Multiplication of Persistence_heat_maps by scalar (so that all values of the heat map gets multiplied by that scalar). + **/ + Persistence_heat_maps multiply_by_scalar( double scalar )const + { + Persistence_heat_maps result; + result.min_ = this->min_; + result.max_ = this->max_; + result.heat_map.reserve( this->heat_map.size() ); + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + std::vector< double > v; + v.reserve( this->heat_map[i].size() ); + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + v.push_back( this->heat_map[i][j] * scalar ); + } + result.heat_map.push_back( v ); + } + return result; + } + + /** + * This function computes a sum of two objects of a type Persistence_heat_maps. + **/ + friend Persistence_heat_maps operator+( const Persistence_heat_maps& first , const Persistence_heat_maps& second ) + { + return operation_on_pair_of_heat_maps( first , second , std::plus<double>() ); + } + /** + * This function computes a difference of two objects of a type Persistence_heat_maps. + **/ + friend Persistence_heat_maps operator-( const Persistence_heat_maps& first , const Persistence_heat_maps& second ) + { + return operation_on_pair_of_heat_maps( first , second , std::minus<double>() ); + } + /** + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ + friend Persistence_heat_maps operator*( double scalar , const Persistence_heat_maps& A ) + { + return A.multiply_by_scalar( scalar ); + } + /** + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ + friend Persistence_heat_maps operator*( const Persistence_heat_maps& A , double scalar ) + { + return A.multiply_by_scalar( scalar ); + } + /** + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ + Persistence_heat_maps operator*( double scalar ) + { + return this->multiply_by_scalar( scalar ); + } + /** + * += operator for Persistence_heat_maps. + **/ + Persistence_heat_maps operator += ( const Persistence_heat_maps& rhs ) + { + *this = *this + rhs; + return *this; + } + /** + * -= operator for Persistence_heat_maps. + **/ + Persistence_heat_maps operator -= ( const Persistence_heat_maps& rhs ) + { + *this = *this - rhs; + return *this; + } + /** + * *= operator for Persistence_heat_maps. + **/ + Persistence_heat_maps operator *= ( double x ) + { + *this = *this*x; + return *this; + } + /** + * /= operator for Persistence_heat_maps. + **/ + Persistence_heat_maps operator /= ( double x ) + { + if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." ); + *this = *this * (1/x); + return *this; + } + + + //Implementations of functions for various concepts. + + /** + * 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 )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()const + { + return this->number_of_functions_for_vectorization; + } + + /** + * This function is required by the Real_valued_topological_data concept. It returns various projections od the persistence heat map to a real line. + * At the moment this function is not tested, since it is quite likelly to be changed in the future. Given this, when using it, keep in mind that it + * will be most likelly changed in the next versions. + **/ + 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()const + { + return this->number_of_functions_for_projections_to_reals; + } + + /** + * A function to compute distance between persistence heat maps. + * 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. + * For max norm distance, set power to std::numeric_limits<double>::max() + **/ + 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. + * This function is required by Topological_data_with_averages concept. + **/ + void compute_average( const std::vector< Persistence_heat_maps* >& to_average ); + + /** + * A function to compute scalar product of persistence heat maps. + * 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_ )const; + + //end of implementation of functions needed for concepts. + + + /** + * The x-range of the persistence heat map. + **/ + std::pair< double , double > get_x_range()const + { + return std::make_pair( this->min_ , this->max_ ); + } + + /** + * The y-range of the persistence heat map. + **/ + std::pair< double , double > get_y_range()const + { + return this->get_x_range(); + } + + + + +protected: + //private methods + std::vector< std::vector<double> > check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps ); + size_t number_of_functions_for_vectorization; + size_t number_of_functions_for_projections_to_reals; + void construct( const std::vector< std::pair<double,double> >& intervals_ , + std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1), + + bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = std::numeric_limits<double>::max() , double max_ = std::numeric_limits<double>::max() ); + + void set_up_parameters_for_basic_classes() + { + this->number_of_functions_for_vectorization = 1; + this->number_of_functions_for_projections_to_reals = 1; + } + + //data + //double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram ); + Scalling_of_kernels f; + bool erase_below_diagonal; + double min_; + double max_; + std::vector< std::vector< double > > heat_map; +}; + + +//if min_ == max_, then the program is requested to set up the values itself based on persistence intervals +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::construct( const std::vector< std::pair<double,double> >& intervals_ , + std::vector< std::vector<double> > filter, + bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ ) +{ + bool dbg = false; + if ( dbg )std::cerr << "Entering construct procedure \n"; + Scalling_of_kernels f; + this->f = f; + + if ( dbg )std::cerr << "min and max passed to construct() procedure: " << min_ << " " << max_ << std::endl; + + if ( min_ == max_ ) + { + if (dbg)std::cerr << "min and max parameters will be etermined based on intervals \n"; + //in this case, we want the program to set up the min_ and max_ values by itself. + min_ = std::numeric_limits<int>::max(); + max_ = -std::numeric_limits<int>::max(); + + + for ( size_t i = 0 ; i != intervals_.size() ; ++i ) + { + if ( intervals_[i].first < min_ )min_ = intervals_[i].first; + if ( intervals_[i].second > max_ )max_ = intervals_[i].second; + } + //now we have the structure filled in, and moreover we know min_ and max_ values of the interval, so we know the range. + + //add some more space: + min_ -= fabs(max_ - min_)/100; + max_ += fabs(max_ - min_)/100; + } + + if ( dbg ) + { + std::cerr << "min_ : " << min_ << std::endl; + std::cerr << "max_ : " << max_ << std::endl; + std::cerr << "number_of_pixels : " << number_of_pixels << std::endl; + getchar(); + } + + this->min_ = min_; + this->max_ = max_; + + //initialization of the structure heat_map + std::vector< std::vector<double> > heat_map_; + for ( size_t i = 0 ; i != number_of_pixels ; ++i ) + { + std::vector<double> v( number_of_pixels , 0 ); + heat_map_.push_back( v ); + } + this->heat_map = heat_map_; + + if (dbg)std::cerr << "Done creating of the heat map, now we will fill in the structure \n"; + + for ( size_t pt_nr = 0 ; pt_nr != intervals_.size() ; ++pt_nr ) + { + //compute the value of intervals_[pt_nr] in the grid: + int x_grid = (int)((intervals_[pt_nr].first - this->min_)/( this->max_-this->min_ )*number_of_pixels); + int y_grid = (int)((intervals_[pt_nr].second - this->min_)/( this->max_-this->min_ )*number_of_pixels); + + if ( dbg ) + { + std::cerr << "point : " << intervals_[pt_nr].first << " , " << intervals_[pt_nr].second << std::endl; + std::cerr << "x_grid : " << x_grid << std::endl; + std::cerr << "y_grid : " << y_grid << std::endl; + } + + //x_grid and y_grid gives a center of the kernel. We want to have its lower left cordner. To get this, we need to shift x_grid and y_grid by a grid diameter. + x_grid -= filter.size()/2; + y_grid -= filter.size()/2; + //note that the numbers x_grid and y_grid may be negative. + + if ( dbg ) + { + std::cerr << "After shift : \n";; + std::cerr << "x_grid : " << x_grid << std::endl; + std::cerr << "y_grid : " << y_grid << std::endl; + } + + double scaling_value = this->f(intervals_[pt_nr]); + + + for ( size_t i = 0 ; i != filter.size() ; ++i ) + { + for ( size_t j = 0 ; j != filter.size() ; ++j ) + { + //if the point (x_grid+i,y_grid+j) is the correct point in the grid. + if ( + ((x_grid+i)>=0) && (x_grid+i<this->heat_map.size()) + && + ((y_grid+j)>=0) && (y_grid+j<this->heat_map.size()) + ) + { + if ( dbg ){std::cerr << y_grid+j << " " << x_grid+i << std::endl;} + this->heat_map[ y_grid+j ][ x_grid+i ] += scaling_value * filter[i][j]; + if ( dbg ) + { + std::cerr << "Position : (" << x_grid+i << "," << y_grid+j << ") got increased by the value : " << filter[i][j] << std::endl; + } + } + } + } + + } + + //now it remains to cut everything below diagonal if the user wants us to. + if ( erase_below_diagonal ) + { + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = i ; j != this->heat_map.size() ; ++j ) + { + this->heat_map[i][j] = 0; + } + } + } +}//construct + +template <typename Scalling_of_kernels> +Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( const std::vector< std::pair< double,double > > & interval , + std::vector< std::vector<double> > filter, + bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ ) +{ + this->construct( interval , filter , erase_below_diagonal , number_of_pixels , min_ , max_ ); + this->set_up_parameters_for_basic_classes(); +} + + +template <typename Scalling_of_kernels> +Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( const char* filename , + std::vector< std::vector<double> > filter, + bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ , unsigned dimension ) +{ + std::vector< std::pair< double , double > > intervals_; + if ( dimension == std::numeric_limits<unsigned>::max() ) + { + intervals_ = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + intervals_ = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + //std::cerr << "intervals_.size() : " << intervals_.size() << std::endl; + //for ( size_t i = 0 ; i != intervals_.size() ; ++i ) + //{ + // std::cerr << intervals_[i].first << " " << intervals_[i].second << std::endl; + //} + this->construct( intervals_ , filter, erase_below_diagonal , number_of_pixels , min_ , max_ ); + this->set_up_parameters_for_basic_classes(); +} + + +template <typename Scalling_of_kernels> +std::vector< std::vector<double> > Persistence_heat_maps<Scalling_of_kernels>::check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps ) +{ + //checking if all the heat maps are of the same size: + for ( size_t i = 0 ; i != maps.size() ; ++i ) + { + if ( maps[i]->heat_map.size() != maps[0]->heat_map.size() ) + { + std::cerr << "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n"; + throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n"; + } + if ( maps[i]->heat_map[0].size() != maps[0]->heat_map[0].size() ) + { + std::cerr << "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n"; + throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n"; + } + } + std::vector< std::vector<double> > heat_maps( maps[0]->heat_map.size() ); + for ( size_t i = 0 ; i != maps[0]->heat_map.size() ; ++i ) + { + std::vector<double> v( maps[0]->heat_map[0].size() , 0 ); + heat_maps[i] = v; + } + return heat_maps; +} + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::compute_median( const std::vector<Persistence_heat_maps*>& maps ) +{ + std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps ); + + std::vector<double> to_compute_median( maps.size() ); + for ( size_t i = 0 ; i != heat_maps.size() ; ++i ) + { + for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j ) + { + for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no ) + { + to_compute_median[map_no] = maps[map_no]->heat_map[i][j]; + } + std::nth_element(to_compute_median.begin(), to_compute_median.begin() + to_compute_median.size()/2, to_compute_median.end()); + heat_maps[i][j] = to_compute_median[ to_compute_median.size()/2 ]; + } + } + this->heat_map = heat_maps; + this->min_= maps[0]->min_; + this->max_= maps[0]->max_; +} + + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::compute_mean( const std::vector<Persistence_heat_maps*>& maps ) +{ + std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps ); + for ( size_t i = 0 ; i != heat_maps.size() ; ++i ) + { + for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j ) + { + double mean = 0; + for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no ) + { + mean += maps[map_no]->heat_map[i][j]; + } + heat_maps[i][j] = mean/(double)maps.size(); + } + } + this->heat_map = heat_maps; + this->min_ = maps[0]->min_; + this->max_ = maps[0]->max_; +} + + + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::compute_percentage_of_active( const std::vector<Persistence_heat_maps*>& maps , size_t cutoff ) +{ + std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps ); + + for ( size_t i = 0 ; i != heat_maps.size() ; ++i ) + { + for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j ) + { + size_t number_of_active_levels = 0; + for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no ) + { + if ( maps[map_no]->heat_map[i][j] ) number_of_active_levels++; + } + if ( number_of_active_levels > cutoff ) + { + heat_maps[i][j] = number_of_active_levels; + } + else + { + heat_maps[i][j] = 0; + } + } + } + this->heat_map = heat_maps; + this->min_ = maps[0]->min_; + this->max_ = maps[0]->max_; +} + + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::plot( const char* filename )const +{ + std::ofstream out; + std::stringstream ss; + ss << filename << "_GnuplotScript"; + + out.open( ss.str().c_str() ); + out << "plot '-' matrix with image" << std::endl; + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + out << this->heat_map[i][j] << " "; + } + out << std::endl; + } + out.close(); + std::cout << "Gnuplot script have been created. Open gnuplot and type load \'" << ss.str().c_str() << "\' to see the picture." << std::endl; +} + + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::print_to_file( const char* filename )const +{ + + std::ofstream out; + out.open( filename ); + + //First we store this->min_ and this->max_ values: + out << this->min_ << " " << this->max_ << std::endl; + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + out << this->heat_map[i][j] << " "; + } + out << std::endl; + } + out.close(); +} + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::load_from_file( const char* filename ) +{ + bool dbg = false; + + std::ifstream in; + in.open( filename ); + + //checking if the file exist / if it was open. + if ( !in.good() ) + { + 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"; + } + + //now we read the file one by one. + + + + in >> this->min_ >> this->max_; + if ( dbg ) + { + std::cerr << "Reading the following values of min and max : " << this->min_ << " , " << this->max_ << std::endl; + } + + std::string temp; + std::getline(in, temp); + + while (!in.eof()) + { + std::getline(in, temp); + std::stringstream lineSS; + lineSS << temp; + + std::vector<double> line_of_heat_map; + while ( lineSS.good() ) + { + double point; + + lineSS >> point; + line_of_heat_map.push_back( point ); + if ( dbg ) + { + std::cout << point << " "; + } + } + if ( dbg ) + { + std::cout << std::endl; + getchar(); + } + + if ( in.good() )this->heat_map.push_back( line_of_heat_map ); + } + in.close(); + if ( dbg )std::cout << "Done \n"; +} + + +//Concretizations of virtual methods: +template <typename Scalling_of_kernels> +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; + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + size_of_result += this->heat_map[i].size(); + } + + std::vector< double > result; + result.reserve( size_of_result ); + + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + result.push_back( this->heat_map[i][j] ); + } + } + + return result; +} + +template <typename Scalling_of_kernels> +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) ) + { + std::cerr << "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate"; + throw "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate"; + } + + //if we are here, we know that the two persistence iomages are defined on the same domain, so we can start computing their distances: + + double distance = 0; + if ( power < std::numeric_limits<double>::max() ) + { + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + distance += pow( fabs(this->heat_map[i][j] - second.heat_map[i][j]) , power ); + } + } + } + else + { + //in this case, we compute max norm distance + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + if ( distance < fabs(this->heat_map[i][j] - second.heat_map[i][j]) ) + { + distance = fabs(this->heat_map[i][j] - second.heat_map[i][j]); + } + } + } + } + return distance; +} + +template <typename Scalling_of_kernels> +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 ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + result += this->heat_map[i][j]; + } + } + return result; +} + +template <typename Scalling_of_kernels> +void Persistence_heat_maps<Scalling_of_kernels>::compute_average( const std::vector< Persistence_heat_maps* >& to_average ) +{ + this->compute_mean( to_average ); +} + +template <typename Scalling_of_kernels> +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) ) + { + std::cerr << "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate"; + throw "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate"; + } + + //if we are here, we know that the two persistence iomages are defined on the same domain, so we can start computing their scalar product: + double scalar_prod = 0; + for ( size_t i = 0 ; i != this->heat_map.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j ) + { + scalar_prod += this->heat_map[i][j]*second.heat_map[i][j]; + } + } + return scalar_prod; +} + + + + +}//namespace Gudhi_stat +}//namespace Gudhi + + +#endif diff --git a/src/Gudhi_stat/include/gudhi/Persistence_intervals.h b/src/Gudhi_stat/include/gudhi/Persistence_intervals.h new file mode 100644 index 00000000..792c0a28 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Persistence_intervals.h @@ -0,0 +1,700 @@ +/* This file is part of the Gudhi hiLibrary. 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 Persistence_intervals_H_ +#define Persistence_intervals_H_ + +//gudhi include +#include <gudhi/read_persistence_from_file.h> + +//standard include +#include <limits> +#include <iostream> +#include <fstream> +#include <vector> +#include <algorithm> +#include <limits> +#include <cmath> +#include <functional> + +namespace Gudhi +{ +namespace Persistence_representations +{ + +/** + * This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data +**/ +class Persistence_intervals +{ +public: + /** + * This is a constructor of a class Persistence_intervals from a text file. Each line of the input file is supposed to contain two numbers of a type doube (or convertable to double) + * representing the birth and the death of the persistence interval. If the pairs are not sorted so that birth <= death, then the constructor will sort then that way. + * * The second parameter of a constructor is a dimension of intervals to be read from a file. If your file contains only birt-death pairs, use the default value. + **/ + Persistence_intervals( const char* filename , unsigned dimension = std::numeric_limits<unsigned>::max() ); + + /** + * This is a constructor of a class Persistence_intervals from a vector of pairs. Each pair is assumed to represent a persistence interval. We assume that the first elemnets of pairs + * are smaller or equal the second elements of pairs. + **/ + Persistence_intervals( const std::vector< std::pair< double,double > >& intervals ); + + /** + * This procedure returns x-range of a given persistence diagram. + **/ + std::pair< double , double > get_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].second > max_ )max_ = this->intervals[i].second; + } + return std::make_pair( min_ , max_ ); + } + + /** + * This procedure returns y-range of a given persistence diagram. + **/ + std::pair< double , double > get_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 )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 )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 > histogram_of_lengths( size_t number_of_bins = 10 )const; + + /** + * Based on a histogram of intervals lengts computed by the function histogram_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_histogram_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 )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 )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()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 )const; + + /** + * Operator that send the diagram to a stream. + **/ + friend std::ostream& operator << ( std::ostream& out , const Persistence_intervals& intervals ) + { + for ( size_t i = 0 ; i != intervals.intervals.size() ; ++i ) + { + out << intervals.intervals[i].first << " " << intervals.intervals[i].second << std::endl; + } + return out; + } + + /** + * Generating gnuplot script to plot the interval. + **/ + void plot( const char* filename , 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() ) const + { + //this program create a gnuplot script file that allows to plot persistence diagram. + std::ofstream out; + + std::ostringstream nameSS; + nameSS << filename << "_GnuplotScript"; + std::string nameStr = nameSS.str(); + out.open( nameStr ); + + std::pair<double,double> min_max_values = this->get_x_range(); + 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 ) + { + out << this->intervals[i].first << " " << this->intervals[i].second << std::endl; + } + out << "EOF" << std::endl; + out << min_max_values.first - 0.1*(min_max_values.second-min_max_values.first) << " " << min_max_values.first - 0.1*(min_max_values.second-min_max_values.first) << std::endl; + out << min_max_values.second + 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.close(); + + std::cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << std::endl; + } + + + + + /** + * Retun numbr of points in the diagram. + **/ + size_t size()const{return this->intervals.size();} + + /** + * 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)const + { + if ( i >= this->intervals.size() )throw("Index out of range! Operator [], one_d_gaussians class\n"); + return this->intervals[i]; + } + + + + + + + + + + + + + + + + //Implementations of functions for various concepts. + /** + * 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. + * At the moment this function is not tested, since it is quite likelly to be changed in the future. Given this, when using it, keep in mind that it + * will be most likelly changed in the next versions. + **/ + 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()const + { + return this->number_of_functions_for_projections_to_reals; + } + + /** + * 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 )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()const + { + return this->number_of_functions_for_vectorization; + } + + //end of implementation of functions needed for concepts. + //end of implementation of functions needed for concepts. + + + + + + + + + + + + + + //For visualization use output from vectorize and build histograms. + std::vector< std::pair< double,double > > output_for_visualization() + { + return this->intervals; + } + +protected: + + void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals() + { + //warning, this function can be only called after filling in the intervals vector. + this->number_of_functions_for_vectorization = this->intervals.size(); + this->number_of_functions_for_projections_to_reals = 1; + } + + std::vector< std::pair< double,double > > intervals; + size_t number_of_functions_for_vectorization; + size_t number_of_functions_for_projections_to_reals; +}; + + +Persistence_intervals::Persistence_intervals( const char* filename , unsigned dimension ) +{ + //bool dbg = false; + //ifstream in; + //in.open( filename ); + + //if ( !in.good() ) + //{ + // throw("File with the persistence diagram do not exist, the program will now terminate.\n"); + //} + + //while ( true ) + //{ + // double first; + // double second; + // in >> first >> second; + + // if ( first > second ) + // { + // double buf = first; + // first = second; + // second = buf; + // } + + // if ( in.eof() )break; + // this->intervals.push_back( std::make_pair( first,second ) ); + // if ( dbg ) + // { + // std::cerr << "Adding interval [ " << first << " , " << second << " ]\n"; + // getchar(); + // } + //} + //in.close(); + if ( dimension == std::numeric_limits<unsigned>::max() ) + { + this->intervals = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + this->intervals = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); +}//Persistence_intervals + + +Persistence_intervals::Persistence_intervals( const std::vector< std::pair< double , double > >& intervals_ ):intervals(intervals_) +{ + 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 ) + { + result[i] = this->intervals[i].second - this->intervals[i].first; + } + std::sort( result.begin() , result.end() , std::greater<double>() ); + + + result.resize( std::min(where_to_cut,result.size()) ); + return result; +}//length_of_dominant_intervals + + + +bool compare( const std::pair< size_t , double >& first , const std::pair< size_t , double >& second ) +{ + return first.second > second.second; +} + + +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() ); + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + position_length_vector[i] = std::make_pair( i , this->intervals[i].second - this->intervals[i].first ); + } + + std::sort( position_length_vector.begin() , position_length_vector.end() , compare ); + + std::vector< std::pair<double,double> > result; + result.reserve( std::min( where_to_cut , position_length_vector.size() ) ); + + for ( size_t i = 0 ; i != std::min( where_to_cut , position_length_vector.size() ) ; ++i ) + { + result.push_back( this->intervals[ position_length_vector[i].first ] ); + if ( dbg )std::cerr << "Position : " << position_length_vector[i].first << " length : " << position_length_vector[i].second << std::endl; + } + + return result; +}//dominant_intervals + + +std::vector< size_t > Persistence_intervals::histogram_of_lengths( size_t number_of_bins )const +{ + bool dbg = false; + + if ( dbg )std::cerr << "this->intervals.size() : " << this->intervals.size() << std::endl; + //first find the length of the longest interval: + double lengthOfLongest = 0; + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + if ( (this->intervals[i].second - this->intervals[i].first) > lengthOfLongest ) + { + lengthOfLongest = this->intervals[i].second - this->intervals[i].first; + } + } + + if ( dbg ){std::cerr << "lengthOfLongest : " << lengthOfLongest << std::endl;} + + //this is a container we will use to store the resulting histogram + std::vector< size_t > result( number_of_bins + 1 , 0 ); + + //for every persistence interval in our collection. + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + //compute its length relative to the length of the dominant interval: + double relative_length_of_this_interval = (this->intervals[i].second - this->intervals[i].first)/lengthOfLongest; + + //given the relative length (between 0 and 1) compute to which bin should it contribute. + size_t position = (size_t)(relative_length_of_this_interval*number_of_bins); + + + ++result[position]; + + if ( dbg ) + { + std::cerr << "i : " << i << std::endl; + std::cerr << "Interval : [" << this->intervals[i].first << " , " << this->intervals[i].second << " ] \n"; + std::cerr << "relative_length_of_this_interval : " << relative_length_of_this_interval << std::endl; + std::cerr << "position : " << position << std::endl; + getchar(); + } + } + + + if ( dbg ){for ( size_t i = 0 ; i != result.size() ; ++i )std::cerr << result[i] << std::endl;} + return result; +} + + +std::vector< size_t > Persistence_intervals::cumulative_histogram_of_lengths( size_t number_of_bins )const +{ + std::vector< size_t > histogram = this->histogram_of_lengths( number_of_bins ); + std::vector< size_t > result( histogram.size() ); + + size_t sum = 0; + for ( size_t i = 0 ; i != histogram.size() ; ++i ) + { + sum += histogram[i]; + result[i] = sum; + } + return result; +} + + +std::vector< double > Persistence_intervals::characteristic_function_of_diagram( double x_min , double x_max , size_t number_of_bins )const +{ + bool dbg = false; + + std::vector< double > result( number_of_bins ); + std::fill( result.begin() , result.end() , 0 ); + + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + if ( dbg ) + { + std::cerr << "Interval : " << this->intervals[i].first << " , " << this->intervals[i].second << std::endl; + } + + size_t beginIt = 0; + if ( this->intervals[i].first < x_min )beginIt = 0; + if ( this->intervals[i].first >= x_max )beginIt = result.size(); + if ( ( this->intervals[i].first > x_min ) && ( this->intervals[i].first < x_max ) ) + { + beginIt = number_of_bins*(this->intervals[i].first-x_min)/(x_max - x_min); + } + + size_t endIt = 0; + if ( this->intervals[i].second < x_min )endIt = 0; + if ( this->intervals[i].second >= x_max )endIt = result.size(); + if ( ( this->intervals[i].second > x_min ) && ( this->intervals[i].second < x_max ) ) + { + endIt = number_of_bins*( this->intervals[i].second - x_min )/(x_max - x_min); + } + + if ( beginIt > endIt ){beginIt = endIt;} + + if ( dbg ) + { + std::cerr << "beginIt : " << beginIt << std::endl; + std::cerr << "endIt : " << endIt << std::endl; + } + + + for ( size_t pos = beginIt ; pos != endIt ; ++pos ) + { + result[pos] += ( (x_max - x_min)/(double)number_of_bins ) * ( this->intervals[i].second - this->intervals[i].first ); + } + //cerr << "x_max : " << x_max << " x_min : " << x_min << " , number_of_bins : " << number_of_bins << " this->intervals[i].second : " << this->intervals[i].second << " this->intervals[i].first : " << this->intervals[i].first << endl; + if ( dbg ) + { + std::cerr << "Result at this stage \n"; + for ( size_t aa = 0 ; aa != result.size() ; ++aa ) + { + std::cerr << result[aa] << " "; + } + std::cerr << std::endl; + //getchar(); + } + } + return result; +}//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 )const +{ + std::vector< double > intsOfBars = this->characteristic_function_of_diagram( x_min , x_max , number_of_bins ); + std::vector< double > result( intsOfBars.size() ); + double sum = 0; + for ( size_t i = 0 ; i != intsOfBars.size() ; ++i ) + { + sum += intsOfBars[i]; + result[i] = sum; + } + return result; +}//cumulative_characteristic_function_of_diagram + + +template <typename T> +bool compare_first_element_of_pair( const std::pair< T , bool >& f, const std::pair< T , bool >& s ) +{ + return (f.first < s.first); +} + + +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() ); + + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + places_where_pbs_change[2*i] = std::make_pair( this->intervals[i].first , true ); + places_where_pbs_change[2*i+1] = std::make_pair( this->intervals[i].second , false ); + } + + std::sort( places_where_pbs_change.begin() , places_where_pbs_change.end() , compare_first_element_of_pair<double> ); + size_t pbn = 0; + std::vector< std::pair< double , size_t > > pbns( places_where_pbs_change.size() ); + for ( size_t i = 0 ; i != places_where_pbs_change.size() ; ++i ) + { + if ( places_where_pbs_change[i].second == true ) + { + ++pbn; + } + else + { + --pbn; + } + pbns[i] = std::make_pair( places_where_pbs_change[i].first , pbn ); + } + return pbns; +} + + + + + + +inline double compute_euclidean_distance( const std::pair< double,double > & f , const std::pair< double,double > & s ) +{ + return sqrt( (f.first-s.first)*(f.first-s.first) + (f.second-s.second)*(f.second-s.second) ); +} + + +std::vector< double > Persistence_intervals::k_n_n( size_t k , size_t where_to_cut )const +{ + bool dbg = false; + if ( dbg ) + { + std::cerr << "Here are the intervals : \n"; + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + std::cerr << "[ " << this->intervals[i].first << " , " << this->intervals[i].second << "] \n"; + } + getchar(); + } + + std::vector< double > result; + //compute all to all distance between point in the diagram. Also, consider points in the diagonal with the infinite multiplicity. + std::vector< std::vector< double > > distances( this->intervals.size() ); + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + std::vector<double> aa(this->intervals.size()); + std::fill( aa.begin() , aa.end() , 0 ); + distances[i] = aa; + } + std::vector< double > distances_from_diagonal( this->intervals.size() ); + std::fill( distances_from_diagonal.begin() , distances_from_diagonal.end() , 0 ); + + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + std::vector< double > distancesFromI; + for ( size_t j = i+1 ; j != this->intervals.size() ; ++j ) + { + distancesFromI.push_back( compute_euclidean_distance( this->intervals[i] , this->intervals[j] ) ); + } + //distances.push_back( distancesFromI ); + //also add a distance from this guy to daigonal: + double distanceToDiagonal = compute_euclidean_distance( this->intervals[i] , std::make_pair( 0.5*(this->intervals[i].first + this->intervals[i].second) , 0.5*(this->intervals[i].first + this->intervals[i].second) ) ); + distances_from_diagonal[i] = distanceToDiagonal; + + if ( dbg ) + { + std::cerr << "Here are the distances form the point : [" << this->intervals[i].first << " , " << this->intervals[i].second << "] in the diagram \n"; + for ( size_t aa = 0 ; aa != distancesFromI.size() ; ++aa ) + { + std::cerr << "To : " << i+aa << " : " << distancesFromI[aa] << " "; + } + std::cerr << std::endl; + getchar(); + } + + //filling in the distances matrix: + for ( size_t j = i+1 ; j != this->intervals.size() ; ++j ) + { + distances[i][j] = distancesFromI[j-i-1]; + distances[j][i] = distancesFromI[j-i-1]; + } + } + if ( dbg ) + { + std::cerr << "Here is the distance matrix : \n"; + for ( size_t i = 0 ; i != distances.size() ; ++i ) + { + for ( size_t j = 0 ; j != distances.size() ; ++j ) + { + std::cerr << distances[i][j] << " "; + } + std::cerr << std::endl; + } + std::cerr << std::endl << std::endl << "And here are the distances to the diagonal : " << std::endl; + for ( size_t i = 0 ; i != distances_from_diagonal. size() ; ++i ) + { + std::cerr << distances_from_diagonal[i] << " "; + } + std::cerr << std::endl << std::endl; + getchar(); + } + + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + std::vector< double > distancesFromI = distances[i]; + distancesFromI.push_back( distances_from_diagonal[i] ); + + //sort it: + std::sort( distancesFromI.begin() , distancesFromI.end() , std::greater<double>() ); + + if ( k > distancesFromI.size() ) + { + if ( dbg ) + { + std::cerr << "There are not enough neighbors in your set. We set the result to plus infty \n"; + } + result.push_back( std::numeric_limits<double>::max() ); + } + else + { + if ( distances_from_diagonal[i] > distancesFromI[k] ) + { + if ( dbg ) + { + std::cerr << "The k-th n.n. is on a diagonal. Therefore we set up a distance to diagonal \n"; + } + result.push_back( distances_from_diagonal[i] ); + } + else + { + result.push_back( distancesFromI[k] ); + } + } + } + std::sort( result.begin() , result.end() , std::greater<double>() ); + result.resize( std::min( result.size() , where_to_cut ) ); + + return result; +} + + +double Persistence_intervals::project_to_R( int number_of_function )const +{ + double result = 0; + + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + result += ( this->intervals[i].second - this->intervals[i].first )*( this->intervals[i].second - this->intervals[i].first ); + } + + return result; +} + + +}//namespace gudhi stat +}//namespace gudhi + +#endif diff --git a/src/Gudhi_stat/include/gudhi/Persistence_intervals_with_distances.h b/src/Gudhi_stat/include/gudhi/Persistence_intervals_with_distances.h new file mode 100644 index 00000000..7ef711e9 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Persistence_intervals_with_distances.h @@ -0,0 +1,65 @@ +/* This file is part of the Gudhi hiLibrary. 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 Persistence_intervals_WITH_DISTANCES_H_ +#define Persistence_intervals_WITH_DISTANCES_H_ + + +#include <gudhi/Persistence_intervals.h> +#include <gudhi/Bottleneck.h> + +namespace Gudhi +{ +namespace Persistence_representations +{ + +class Persistence_intervals_with_distances : public Persistence_intervals +{ +public: + using Persistence_intervals::Persistence_intervals; + + /** + *Computations of distance from the current persistnce diagram to the persistence diagram given as a parameter of this function. + *The last but one parameter, power, is here in case we would like to compute p=th Wasserstein distance. At the moment, this method only implement Bottleneck distance, + * which is infinity Wasserstein distance. Therefore any power which is not the default std::numeric_limits< double >::max() will be ignored and an + * exception will be thrown. + * The last parameter, tolerance, it is an additiv error of the approimation, set by default to zero. + **/ + double distance( const Persistence_intervals_with_distances& second , double power = std::numeric_limits< double >::max() , double tolerance = 0) const + { + if ( power >= std::numeric_limits< double >::max() ) + { + return Gudhi::persistence_diagram::bottleneck_distance(this->intervals, second.intervals, tolerance); + } + else + { + std::cerr << "At the moment Gudhi do not support Wasserstein distances. We only support Bottleneck distance." << std::endl; + throw "At the moment Gudhi do not support Wasserstein distances. We only support Bottleneck distance."; + } + } +}; + + +}//namespace gudhi stat +}//namespace gudhi + +#endif diff --git a/src/Gudhi_stat/include/gudhi/Persistence_landscape.h b/src/Gudhi_stat/include/gudhi/Persistence_landscape.h new file mode 100644 index 00000000..9a177b60 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Persistence_landscape.h @@ -0,0 +1,1498 @@ +/* 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 Persistence_landscapes_H +#define Persistence_landscapes_H + +//standard include +#include <cmath> +#include <iostream> +#include <vector> +#include <limits> +#include <fstream> +#include <sstream> +#include <algorithm> + + +//gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> + + + + +namespace Gudhi +{ +namespace Persistence_representations +{ + + + +//predeclaration +class Persistence_landscape; +template < typename operation > +Persistence_landscape operation_on_pair_of_landscapes( const Persistence_landscape& land1 , const Persistence_landscape& land2 ); + + + +/** + * A clas implementing persistence landascpes data structures. For theroretical desciritpion, please consult a paper ''Statistical topological data analysis using persistence landscapes'' by Peter Bubenik. + * For details of algorithms, please consult ''A persistence landscapes toolbox for topological statistics'' by Peter Bubenik and Pawel Dlotko. + * Persistence landscapes allow vertorization, computations of distances, computations of projections to Real, computations of averages and scalar products. Therefore they implement suitable interfaces. + * It implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product + * Note that at the moment, due to roundoff errors during the construction of persistence landscapes, elements which are different by 0.000005 are considered the same. If the scale in your persistence diagrams + * is comparable to this value, please rescale them before use this code. +**/ +class Persistence_landscape +{ +public: + /** + * Default constructor. + **/ + Persistence_landscape() + { + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); + } + + /** + * Constructor that takes as an input a vector of birth-death pairs. + **/ + Persistence_landscape( const std::vector< std::pair< double , double > >& p ); + + /** + * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed + * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. + **/ + Persistence_landscape(const char* filename , size_t dimension = std::numeric_limits<unsigned>::max() ); + + + + /** + * This procedure loads a landscape from file. It erase all the data that was previously stored in this landscape. + **/ + void load_landscape_from_file( const char* filename ); + + + /** + * The procedure stores a landscape to a file. The file can be later used by a procedure load_landscape_from_file. + **/ + void print_to_file( const char* filename )const; + + + + /** + * This function compute integral of the landscape (defined formally as sum of integrals on R of all landscape functions) + **/ + double compute_integral_of_landscape()const; + + + /** + * This function compute integral of the 'level'-level of a landscape. + **/ + double compute_integral_of_a_level_of_a_landscape( size_t level )const; + + + /** + * This function compute integral of the landscape p-th power of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions) + **/ + double compute_integral_of_landscape( double p )const;//this function compute integral of p-th power of landscape. + + + /** + * A function that computes the value of a landscape at a given point. The parameters of the function are: unsigned level and double x. + * The procedure will compute the value of the level-landscape at the point x. + **/ + double compute_value_at_a_given_point( unsigned level , double x )const; + + /** + * Writing landscape into a stream. A i-th level landscape starts with a string "lambda_i". Then the discontinuity points of the landscapes follows. + * Shall those points be joined with lines, we will obtain the i-th landscape function. + **/ + friend std::ostream& operator<<(std::ostream& out, Persistence_landscape& land ); + + template < typename operation > + friend Persistence_landscape operation_on_pair_of_landscapes( const Persistence_landscape& land1 , const Persistence_landscape& land2 ); + + + + + /** + *\private A function that compute sum of two landscapes. + **/ + friend Persistence_landscape add_two_landscapes ( const Persistence_landscape& land1 , const Persistence_landscape& land2 ) + { + return operation_on_pair_of_landscapes< std::plus<double> >(land1,land2); + } + + /** + *\private A function that compute difference of two landscapes. + **/ + friend Persistence_landscape subtract_two_landscapes ( const Persistence_landscape& land1 , const Persistence_landscape& land2 ) + { + return operation_on_pair_of_landscapes< std::minus<double> >(land1,land2); + } + + /** + * An operator +, that compute sum of two landscapes. + **/ + friend Persistence_landscape operator+( const Persistence_landscape& first , const Persistence_landscape& second ) + { + return add_two_landscapes( first,second ); + } + + /** + * An operator -, that compute difference of two landscapes. + **/ + friend Persistence_landscape operator-( const Persistence_landscape& first , const Persistence_landscape& second ) + { + return subtract_two_landscapes( first,second ); + } + + /** + * An operator * that allows multipilication of a landscape by a real number. + **/ + friend Persistence_landscape operator*( const Persistence_landscape& first , double con ) + { + return first.multiply_lanscape_by_real_number_not_overwrite(con); + } + + /** + * An operator * that allows multipilication of a landscape by a real number (order of parameters swapped). + **/ + friend Persistence_landscape operator*( double con , const Persistence_landscape& first ) + { + return first.multiply_lanscape_by_real_number_not_overwrite(con); + } + + /** + * Operator +=. The second parameter is persistence landscape. + **/ + Persistence_landscape operator += ( const Persistence_landscape& rhs ) + { + *this = *this + rhs; + return *this; + } + + /** + * Operator -=. The second parameter is a persistence landscape. + **/ + Persistence_landscape operator -= ( const Persistence_landscape& rhs ) + { + *this = *this - rhs; + return *this; + } + + + /** + * Operator *=. The second parameter is a real number by which the y values of all landscape functions are multiplied. The x-values remain unchanged. + **/ + Persistence_landscape operator *= ( double x ) + { + *this = *this*x; + return *this; + } + + /** + * Operator /=. The second parameter is a real number. + **/ + Persistence_landscape operator /= ( double x ) + { + if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." ); + *this = *this * (1/x); + return *this; + } + + /** + * An operator to compare two persistence landscapes. + **/ + bool operator == ( const Persistence_landscape& rhs )const; + + + /** + * An operator to compare two persistence landscapes. + **/ + bool operator != ( const Persistence_landscape& rhs )const + { + return !((*this) == rhs); + } + + + /** + * Computations of maximum (y) value of landscape. + **/ + double compute_maximum()const + { + double maxValue = 0; + if ( this->land.size() ) + { + maxValue = -std::numeric_limits<int>::max(); + for ( size_t i = 0 ; i != this->land[0].size() ; ++i ) + { + if ( this->land[0][i].second > maxValue )maxValue = this->land[0][i].second; + } + } + return maxValue; + } + + + /** + *\private 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; + } + + /** + *\private Computations of a \f$L^i\f$ norm of landscape, where i is the input parameter. + **/ + double compute_norm_of_landscape( double i ) + { + Persistence_landscape l; + if ( i < std::numeric_limits< double >::max() ) + { + return compute_distance_of_landscapes(*this,l,i); + } + else + { + return compute_max_norm_distance_of_landscapes(*this,l); + } + } + + /** + * An operator to compute the value of a landscape in the level 'level' at the argument 'x'. + **/ + double operator()(unsigned level,double x)const{return this->compute_value_at_a_given_point(level,x);} + + /** + *\private Computations of \f$L^{\infty}\f$ distance between two landscapes. + **/ + friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second ); + //friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , unsigned& nrOfLand , double&x , double& y1, double& y2 ); + + + /** + *\private Computations of \f$L^{p}\f$ distance between two landscapes. p is the parameter of the procedure. + **/ + friend double compute_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , double p ); + + + + /** + * Function to compute absolute value of a PL function. The representation of persistence landscapes allow to store general PL-function. When computing distance betwen two landscapes, we compute difference between + * them. In this case, a general PL-function with negative value can appear as a result. Then in order to compute distance, we need to take its absolute value. This is the purpose of this procedure. + **/ + Persistence_landscape abs(); + + /** + * Computes the number of landscape functions. + **/ + size_t size()const{return this->land.size(); } + + /** + * Computate maximal value of lambda-level landscape. + **/ + double find_max( unsigned lambda )const; + + /** + *\private Function to compute inner (scalar) product of two landscapes. + **/ + friend double compute_inner_product( const Persistence_landscape& l1 , const Persistence_landscape& l2 ); + + + + + + + + + + + + + + + + //Implementations of functions for various concepts. + + /** + * 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. + * At the moment this function is not tested, since it is quite likelly to be changed in the future. Given this, when using it, keep in mind that it + * will be most likelly changed in the next versions. + **/ + double project_to_R( int number_of_function )const + { + return this->compute_integral_of_a_level_of_a_landscape( (size_t)number_of_function ); + } + + /** + * 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()const + { + return this->number_of_functions_for_projections_to_reals; + } + + /** + * 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 )const + { + //TODO, think of something smarter over here + std::vector<double> v; + if ( (size_t)number_of_function > this->land.size() ) + { + return v; + } + v.reserve( this->land[number_of_function].size() ); + for ( size_t i = 0 ; i != this->land[number_of_function].size() ; ++i ) + { + v.push_back( this->land[number_of_function][i].second ); + } + return v; + } + /** + * 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()const + { + return this->number_of_functions_for_vectorization; + } + + /** + * A function to compute averaged persistence landscape, based on vector of persistence landscapes. + * This function is required by Topological_data_with_averages concept. + **/ + void compute_average( const std::vector< Persistence_landscape* >& to_average ) + { + bool dbg = false; + + if ( dbg ){std::cerr << "to_average.size() : " << to_average.size() << std::endl;} + + std::vector< Persistence_landscape* > nextLevelMerge( to_average.size() ); + for ( size_t i = 0 ; i != to_average.size() ; ++i ) + { + nextLevelMerge[i] = to_average[i]; + } + bool is_this_first_level = true;//in the loop, we will create dynamically a unmber of intermediate complexes. We have to clean that up, but we cannot erase the initial andscapes we have + //to average. In this case, we simply check if the nextLevelMerge are the input landscapes or the ones created in that loop by usig this extra variable. + + while ( nextLevelMerge.size() != 1 ) + { + if ( dbg ){std::cerr << "nextLevelMerge.size() : " << nextLevelMerge.size() << std::endl;} + std::vector< Persistence_landscape* > nextNextLevelMerge; + nextNextLevelMerge.reserve( to_average.size() ); + for ( size_t i = 0 ; i < nextLevelMerge.size() ; i=i+2 ) + { + if ( dbg ){std::cerr << "i : " << i << std::endl;} + Persistence_landscape* l = new Persistence_landscape; + if ( i+1 != nextLevelMerge.size() ) + { + (*l) = (*nextLevelMerge[i])+(*nextLevelMerge[i+1]); + } + else + { + (*l) = *nextLevelMerge[i]; + } + nextNextLevelMerge.push_back( l ); + } + if ( dbg ){std::cerr << "After this iteration \n";getchar();} + + if ( !is_this_first_level ) + { + //deallocate the memory if the vector nextLevelMerge do not consist of the initial landscapes + for ( size_t i = 0 ; i != nextLevelMerge.size() ; ++i ) + { + delete nextLevelMerge[i]; + } + } + is_this_first_level = false; + nextLevelMerge.swap(nextNextLevelMerge); + } + (*this) = (*nextLevelMerge[0]); + (*this) *= 1/( (double)to_average.size() ); + } + + + /** + * A function to compute distance between persistence landscape. + * The parameter of this functionis a Persistence_landscape. + * This function is required in Topological_data_with_distances concept. + * For max norm distance, set power to std::numeric_limits<double>::max() + **/ + double distance( const Persistence_landscape& second , double power = 1 )const + { + if ( power < std::numeric_limits<double>::max() ) + { + return compute_distance_of_landscapes( *this , second , power ); + } + else + { + return compute_max_norm_distance_of_landscapes( *this , second ); + } + } + + + /** + * A function to compute scalar product of persistence landscapes. + * 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 )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 > get_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 > get_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; + } + + + + //a function used to create a gnuplot script for visualization of landscapes + void plot( const char* filename, double xRangeBegin = std::numeric_limits<double>::max() , double xRangeEnd = std::numeric_limits<double>::max() , + double yRangeBegin = std::numeric_limits<double>::max() , double yRangeEnd = std::numeric_limits<double>::max(), + int from = std::numeric_limits<int>::max(), int to = std::numeric_limits<int>::max() ); + + +protected: + std::vector< std::vector< std::pair<double,double> > > land; + size_t number_of_functions_for_vectorization; + size_t number_of_functions_for_projections_to_reals; + + void construct_persistence_landscape_from_barcode( const std::vector< std::pair< double , double > > & p ); + Persistence_landscape multiply_lanscape_by_real_number_not_overwrite( double x )const; + void multiply_lanscape_by_real_number_overwrite( double x ); + friend double compute_maximal_distance_non_symmetric( const Persistence_landscape& pl1, const Persistence_landscape& pl2 ); + + void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals() + { + //warning, this function can be only called after filling in the intervals vector. + this->number_of_functions_for_vectorization = this->land.size(); + this->number_of_functions_for_projections_to_reals = this->land.size(); + } +}; + + + + + + + +Persistence_landscape::Persistence_landscape(const char* filename , size_t dimension) +{ + std::vector< std::pair< double , double > > barcode; + if ( dimension < std::numeric_limits<double>::max() ) + { + barcode = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + else + { + barcode = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + this->construct_persistence_landscape_from_barcode( barcode ); + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); +} + + +bool operatorEqualDbg = false; +bool Persistence_landscape::operator == ( const Persistence_landscape& rhs )const +{ + if ( this->land.size() != rhs.land.size() ) + { + if (operatorEqualDbg)std::cerr << "1\n"; + return false; + } + for ( size_t level = 0 ; level != this->land.size() ; ++level ) + { + if ( this->land[level].size() != rhs.land[level].size() ) + { + if (operatorEqualDbg)std::cerr << "this->land[level].size() : " << this->land[level].size() << "\n"; + if (operatorEqualDbg)std::cerr << "rhs.land[level].size() : " << rhs.land[level].size() << "\n"; + if (operatorEqualDbg)std::cerr << "2\n"; + return false; + } + for ( size_t i = 0 ; i != this->land[level].size() ; ++i ) + { + if ( !( almost_equal(this->land[level][i].first , rhs.land[level][i].first) && almost_equal(this->land[level][i].second , rhs.land[level][i].second) ) ) + { + //std::cerr<< this->land[level][i].first << " , " << rhs.land[level][i].first << " and " << this->land[level][i].second << " , " << rhs.land[level][i].second << std::endl; + if (operatorEqualDbg)std::cerr << "this->land[level][i] : " << this->land[level][i].first << " " << this->land[level][i].second << "\n"; + if (operatorEqualDbg)std::cerr << "rhs.land[level][i] : " << rhs.land[level][i].first << " " << rhs.land[level][i].second << "\n"; + if (operatorEqualDbg)std::cerr << "3\n"; + return false; + } + } + } + return true; +} + + + + +Persistence_landscape::Persistence_landscape( const std::vector< std::pair< double , double > > & p ) +{ + this->construct_persistence_landscape_from_barcode( p ); + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); +} + + +void Persistence_landscape::construct_persistence_landscape_from_barcode( const std::vector< std::pair< double , double > > & p ) +{ + bool dbg = false; + if ( dbg ){std::cerr << "Persistence_landscape::Persistence_landscape( const std::vector< std::pair< double , double > >& p )" << std::endl;} + + //this is a general algorithm to construct persistence landscapes. + std::vector< std::pair<double,double> > bars; + bars.insert( bars.begin() , p.begin() , p.end() ); + std::sort( bars.begin() , bars.end() , compare_points_sorting ); + + if (dbg) + { + std::cerr << "Bars : \n"; + for ( size_t i = 0 ; i != bars.size() ; ++i ) + { + std::cerr << bars[i].first << " " << bars[i].second << "\n"; + } + getchar(); + } + + std::vector< std::pair<double,double> > characteristicPoints(p.size()); + for ( size_t i = 0 ; i != bars.size() ; ++i ) + { + characteristicPoints[i] = std::make_pair((bars[i].first+bars[i].second)/2.0 , (bars[i].second - bars[i].first)/2.0); + } + std::vector< std::vector< std::pair<double,double> > > Persistence_landscape; + while ( !characteristicPoints.empty() ) + { + if(dbg) + { + for ( size_t i = 0 ; i != characteristicPoints.size() ; ++i ) + { + std::cout << "(" << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n"; + } + std::cin.ignore(); + } + + std::vector< std::pair<double,double> > lambda_n; + lambda_n.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) ); + lambda_n.push_back( std::make_pair(minus_length(characteristicPoints[0]),0) ); + lambda_n.push_back( characteristicPoints[0] ); + + if (dbg) + { + std::cerr << "1 Adding to lambda_n : (" << -std::numeric_limits<int>::max() << " " << 0 << ") , (" << minus_length(characteristicPoints[0]) << " " << 0 << ") , (" << characteristicPoints[0].first << " " << characteristicPoints[0].second << ") \n"; + } + + size_t i = 1; + std::vector< std::pair<double,double> > newCharacteristicPoints; + while ( i < characteristicPoints.size() ) + { + size_t p = 1; + if ( (minus_length(characteristicPoints[i]) >= minus_length(lambda_n[lambda_n.size()-1])) && (birth_plus_deaths(characteristicPoints[i]) > birth_plus_deaths(lambda_n[lambda_n.size()-1])) ) + { + if ( minus_length(characteristicPoints[i]) < birth_plus_deaths(lambda_n[lambda_n.size()-1]) ) + { + std::pair<double,double> point = std::make_pair( (minus_length(characteristicPoints[i])+birth_plus_deaths(lambda_n[lambda_n.size()-1]))/2 , (birth_plus_deaths(lambda_n[lambda_n.size()-1])-minus_length(characteristicPoints[i]))/2 ); + lambda_n.push_back( point ); + if (dbg) + { + std::cerr << "2 Adding to lambda_n : (" << point.first << " " << point.second << ")\n"; + } + + + if ( dbg ) + { + std::cerr << "characteristicPoints[i+p] : " << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << "\n"; + std::cerr << "point : " << point.first << " " << point.second << "\n"; + getchar(); + } + + while ( (i+p < characteristicPoints.size() ) && ( almost_equal(minus_length(point),minus_length(characteristicPoints[i+p])) ) && ( birth_plus_deaths(point) <= birth_plus_deaths(characteristicPoints[i+p]) ) ) + { + newCharacteristicPoints.push_back( characteristicPoints[i+p] ); + if (dbg) + { + std::cerr << "3.5 Adding to newCharacteristicPoints : (" << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << ")\n"; + getchar(); + } + ++p; + } + + + newCharacteristicPoints.push_back( point ); + if (dbg) + { + std::cerr << "4 Adding to newCharacteristicPoints : (" << point.first << " " << point.second << ")\n"; + } + + + while ( (i+p < characteristicPoints.size() ) && ( minus_length(point) <= minus_length(characteristicPoints[i+p]) ) && (birth_plus_deaths(point)>=birth_plus_deaths(characteristicPoints[i+p])) ) + { + newCharacteristicPoints.push_back( characteristicPoints[i+p] ); + if (dbg) + { + std::cerr << "characteristicPoints[i+p] : " << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << "\n"; + std::cerr << "point : " << point.first << " " << point.second << "\n"; + std::cerr << "characteristicPoints[i+p] birth and death : " << minus_length(characteristicPoints[i+p]) << " , " << birth_plus_deaths(characteristicPoints[i+p]) << "\n"; + std::cerr << "point birth and death : " << minus_length(point) << " , " << birth_plus_deaths(point) << "\n"; + + std::cerr << "3 Adding to newCharacteristicPoints : (" << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << ")\n"; + getchar(); + } + ++p; + } + + } + else + { + lambda_n.push_back( std::make_pair( birth_plus_deaths(lambda_n[lambda_n.size()-1]) , 0 ) ); + lambda_n.push_back( std::make_pair( minus_length(characteristicPoints[i]) , 0 ) ); + if (dbg) + { + std::cerr << "5 Adding to lambda_n : (" << birth_plus_deaths(lambda_n[lambda_n.size()-1]) << " " << 0 << ")\n"; + std::cerr << "5 Adding to lambda_n : (" << minus_length(characteristicPoints[i]) << " " << 0 << ")\n"; + } + } + lambda_n.push_back( characteristicPoints[i] ); + if (dbg) + { + std::cerr << "6 Adding to lambda_n : (" << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n"; + } + } + else + { + newCharacteristicPoints.push_back( characteristicPoints[i] ); + if (dbg) + { + std::cerr << "7 Adding to newCharacteristicPoints : (" << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n"; + } + } + i = i+p; + } + lambda_n.push_back( std::make_pair(birth_plus_deaths(lambda_n[lambda_n.size()-1]),0) ); + lambda_n.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) ); + + characteristicPoints = newCharacteristicPoints; + + lambda_n.erase(std::unique(lambda_n.begin(), lambda_n.end()), lambda_n.end()); + this->land.push_back( lambda_n ); + } +} + + + +//this function find maximum of lambda_n +double Persistence_landscape::find_max( unsigned lambda )const +{ + if ( this->land.size() < lambda )return 0; + double maximum = -std::numeric_limits<int>::max(); + for ( size_t i = 0 ; i != this->land[lambda].size() ; ++i ) + { + if ( this->land[lambda][i].second > maximum )maximum = this->land[lambda][i].second; + } + return maximum; +} + + +double Persistence_landscape::compute_integral_of_landscape()const +{ + double result = 0; + for ( size_t i = 0 ; i != this->land.size() ; ++i ) + { + for ( size_t nr = 2 ; nr != this->land[i].size()-1 ; ++nr ) + { + //it suffices to compute every planar integral and then sum them ap for each lambda_n + result += 0.5*( this->land[i][nr].first - this->land[i][nr-1].first )*(this->land[i][nr].second + this->land[i][nr-1].second); + } + } + return result; +} + +double Persistence_landscape::compute_integral_of_a_level_of_a_landscape( size_t level )const +{ + double result = 0; + if ( level >= this->land.size() ) + { + //this landscape function is constantly equal 0, so is the intergral. + return result; + } + //also negative landscapes are assumed to be zero. + if ( level < 0 )return 0; + + for ( size_t nr = 2 ; nr != this->land[ level ].size()-1 ; ++nr ) + { + //it suffices to compute every planar integral and then sum them ap for each lambda_n + result += 0.5*( this->land[ level ][nr].first - this->land[ level ][nr-1].first )*(this->land[ level ][nr].second + this->land[ level ][nr-1].second); + } + + return result; +} + + +double Persistence_landscape::compute_integral_of_landscape( double p )const +{ + bool dbg = false; + double result = 0; + for ( size_t i = 0 ; i != this->land.size() ; ++i ) + { + for ( size_t nr = 2 ; nr != this->land[i].size()-1 ; ++nr ) + { + if (dbg)std::cout << "nr : " << nr << "\n"; + //In this interval, the landscape has a form f(x) = ax+b. We want to compute integral of (ax+b)^p = 1/a * (ax+b)^{p+1}/(p+1) + std::pair<double,double> coef = compute_parameters_of_a_line( this->land[i][nr] , this->land[i][nr-1] ); + double a = coef.first; + double b = coef.second; + + if (dbg)std::cout << "(" << this->land[i][nr].first << "," << this->land[i][nr].second << ") , " << this->land[i][nr-1].first << "," << this->land[i][nr].second << ")" << std::endl; + if ( this->land[i][nr].first == this->land[i][nr-1].first )continue; + if ( a != 0 ) + { + result += 1/(a*(p+1)) * ( pow((a*this->land[i][nr].first+b),p+1) - pow((a*this->land[i][nr-1].first+b),p+1)); + } + else + { + result += ( this->land[i][nr].first - this->land[i][nr-1].first )*( pow(this->land[i][nr].second,p) ); + } + if ( dbg ) + { + std::cout << "a : " <<a << " , b : " << b << std::endl; + std::cout << "result : " << result << std::endl; + } + } + //if (compute_integral_of_landscapeDbg) std::cin.ignore(); + } + return result; +} + + +//this is O(log(n)) algorithm, where n is number of points in this->land. +double Persistence_landscape::compute_value_at_a_given_point( unsigned level , double x )const +{ + bool compute_value_at_a_given_pointDbg = false; + //in such a case lambda_level = 0. + if ( level > this->land.size() ) return 0; + + //we know that the points in this->land[level] are ordered according to x coordinate. Therefore, we can find the point by using bisection: + unsigned coordBegin = 1; + unsigned coordEnd = this->land[level].size()-2; + + if ( compute_value_at_a_given_pointDbg ) + { + std::cerr << "Here \n"; + std::cerr << "x : " << x << "\n"; + std::cerr << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n"; + std::cerr << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n"; + } + + //in this case x is outside the support of the landscape, therefore the value of the landscape is 0. + if ( x <= this->land[level][coordBegin].first )return 0; + if ( x >= this->land[level][coordEnd].first )return 0; + + if (compute_value_at_a_given_pointDbg)std::cerr << "Entering to the while loop \n"; + + while ( coordBegin+1 != coordEnd ) + { + if (compute_value_at_a_given_pointDbg) + { + std::cerr << "coordBegin : " << coordBegin << "\n"; + std::cerr << "coordEnd : " << coordEnd << "\n"; + std::cerr << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n"; + std::cerr << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n"; + } + + + unsigned newCord = (unsigned)floor((coordEnd+coordBegin)/2.0); + + if (compute_value_at_a_given_pointDbg) + { + std::cerr << "newCord : " << newCord << "\n"; + std::cerr << "this->land[level][newCord].first : " << this->land[level][newCord].first << "\n"; + std::cin.ignore(); + } + + if ( this->land[level][newCord].first <= x ) + { + coordBegin = newCord; + if ( this->land[level][newCord].first == x )return this->land[level][newCord].second; + } + else + { + coordEnd = newCord; + } + } + + if (compute_value_at_a_given_pointDbg) + { + std::cout << "x : " << x << " is between : " << this->land[level][coordBegin].first << " a " << this->land[level][coordEnd].first << "\n"; + std::cout << "the y coords are : " << this->land[level][coordBegin].second << " a " << this->land[level][coordEnd].second << "\n"; + std::cerr << "coordBegin : " << coordBegin << "\n"; + std::cerr << "coordEnd : " << coordEnd << "\n"; + std::cin.ignore(); + } + return function_value( this->land[level][coordBegin] , this->land[level][coordEnd] , x ); +} + +std::ostream& operator<<(std::ostream& out, Persistence_landscape& land ) +{ + for ( size_t level = 0 ; level != land.land.size() ; ++level ) + { + out << "Lambda_" << level << ":" << std::endl; + for ( size_t i = 0 ; i != land.land[level].size() ; ++i ) + { + if ( land.land[level][i].first == -std::numeric_limits<int>::max() ) + { + out << "-inf"; + } + else + { + if ( land.land[level][i].first == std::numeric_limits<int>::max() ) + { + out << "+inf"; + } + else + { + out << land.land[level][i].first; + } + } + out << " , " << land.land[level][i].second << std::endl; + } + } + return out; +} + + + + +void Persistence_landscape::multiply_lanscape_by_real_number_overwrite( double x ) +{ + for ( size_t dim = 0 ; dim != this->land.size() ; ++dim ) + { + for ( size_t i = 0 ; i != this->land[dim].size() ; ++i ) + { + this->land[dim][i].second *= x; + } + } +} + +bool AbsDbg = false; +Persistence_landscape Persistence_landscape::abs() +{ + Persistence_landscape result; + for ( size_t level = 0 ; level != this->land.size() ; ++level ) + { + if ( AbsDbg ){ std::cout << "level: " << level << std::endl; } + std::vector< std::pair<double,double> > lambda_n; + lambda_n.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) ); + for ( size_t i = 1 ; i != this->land[level].size() ; ++i ) + { + if ( AbsDbg ){std::cout << "this->land[" << level << "][" << i << "] : " << this->land[level][i].first << " " << this->land[level][i].second << std::endl;} + //if a line segment between this->land[level][i-1] and this->land[level][i] crosses the x-axis, then we have to add one landscape point t oresult + if ( (this->land[level][i-1].second)*(this->land[level][i].second) < 0 ) + { + double zero = find_zero_of_a_line_segment_between_those_two_points( this->land[level][i-1] , this->land[level][i] ); + + lambda_n.push_back( std::make_pair(zero , 0) ); + lambda_n.push_back( std::make_pair(this->land[level][i].first , fabs(this->land[level][i].second)) ); + if ( AbsDbg ) + { + std::cout << "Adding pair : (" << zero << ",0)" << std::endl; + std::cout << "In the same step adding pair : (" << this->land[level][i].first << "," << fabs(this->land[level][i].second) << ") " << std::endl; + std::cin.ignore(); + } + } + else + { + lambda_n.push_back( std::make_pair(this->land[level][i].first , fabs(this->land[level][i].second)) ); + if ( AbsDbg ) + { + std::cout << "Adding pair : (" << this->land[level][i].first << "," << fabs(this->land[level][i].second) << ") " << std::endl; + std::cin.ignore(); + } + } + } + result.land.push_back( lambda_n ); + } + return result; +} + + +Persistence_landscape Persistence_landscape::multiply_lanscape_by_real_number_not_overwrite( double x )const +{ + std::vector< std::vector< std::pair<double,double> > > result(this->land.size()); + for ( size_t dim = 0 ; dim != this->land.size() ; ++dim ) + { + std::vector< std::pair<double,double> > lambda_dim( this->land[dim].size() ); + for ( size_t i = 0 ; i != this->land[dim].size() ; ++i ) + { + lambda_dim[i] = std::make_pair( this->land[dim][i].first , x*this->land[dim][i].second ); + } + result[dim] = lambda_dim; + } + Persistence_landscape res; + //CHANGE + //res.land = result; + res.land.swap(result); + return res; +}//multiply_lanscape_by_real_number_overwrite + + +void Persistence_landscape::print_to_file( const char* filename )const +{ + std::ofstream write; + write.open(filename); + for ( size_t dim = 0 ; dim != this->land.size() ; ++dim ) + { + write << "#lambda_" << dim << std::endl; + for ( size_t i = 1 ; i != this->land[dim].size()-1 ; ++i ) + { + write << this->land[dim][i].first << " " << this->land[dim][i].second << std::endl; + } + } + write.close(); +} + +void Persistence_landscape::load_landscape_from_file( const char* filename ) +{ + bool dbg = false; + //removing the current content of the persistence landscape. + this->land.clear(); + + + //this constructor reads persistence landscape form a file. This file have to be created by this software beforehead + std::ifstream in; + in.open( filename ); + if ( !in.good() ) + { + 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::string line; + std::vector< std::pair<double,double> > landscapeAtThisLevel; + + bool isThisAFirsLine = true; + while ( !in.eof() ) + { + getline(in,line); + if ( !(line.length() == 0 || line[0] == '#') ) + { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + landscapeAtThisLevel.push_back( std::make_pair( beginn , endd ) ); + if (dbg){std::cerr << "Reading a pont : " << beginn << " , " << endd << std::endl;} + } + else + { + if (dbg) + { + std::cout << "IGNORE LINE\n"; + getchar(); + } + if ( !isThisAFirsLine ) + { + landscapeAtThisLevel.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) ); + this->land.push_back(landscapeAtThisLevel); + std::vector< std::pair<double,double> > newLevelOdLandscape; + landscapeAtThisLevel.swap(newLevelOdLandscape); + } + landscapeAtThisLevel.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) ); + isThisAFirsLine = false; + } + } + if ( landscapeAtThisLevel.size() > 1 ) + { + //seems that the last line of the file is not finished with the newline sign. We need to put what we have in landscapeAtThisLevel to the constructed landscape. + landscapeAtThisLevel.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) ); + this->land.push_back(landscapeAtThisLevel); + } + + in.close(); +} + + +template < typename T > +Persistence_landscape operation_on_pair_of_landscapes ( const Persistence_landscape& land1 , const Persistence_landscape& land2 ) +{ + bool operation_on_pair_of_landscapesDBG = false; + if ( operation_on_pair_of_landscapesDBG ){std::cout << "operation_on_pair_of_landscapes\n";std::cin.ignore();} + Persistence_landscape result; + std::vector< std::vector< std::pair<double,double> > > land( std::max( land1.land.size() , land2.land.size() ) ); + result.land = land; + T oper; + + if ( operation_on_pair_of_landscapesDBG ) + { + for ( size_t i = 0 ; i != std::min( land1.land.size() , land2.land.size() ) ; ++i ) + { + std::cerr << "land1.land[" << i << "].size() : " << land1.land[i].size() << std::endl; + std::cerr << "land2.land[" << i << "].size() : " << land2.land[i].size() << std::endl; + } + getchar(); + } + + for ( size_t i = 0 ; i != std::min( land1.land.size() , land2.land.size() ) ; ++i ) + { + std::vector< std::pair<double,double> > lambda_n; + size_t p = 0; + size_t q = 0; + while ( (p+1 < land1.land[i].size()) && (q+1 < land2.land[i].size()) ) + { + if ( operation_on_pair_of_landscapesDBG ) + { + std::cerr << "p : " << p << "\n"; + std::cerr << "q : " << q << "\n"; + std::cerr << "land1.land.size() : " << land1.land.size() << std::endl; + std::cerr << "land2.land.size() : " << land2.land.size() << std::endl; + std::cerr << "land1.land[" << i << "].size() : " << land1.land[i].size() << std::endl; + std::cerr << "land2.land[" << i << "].size() : " << land2.land[i].size() << std::endl; + std::cout << "land1.land[i][p].first : " << land1.land[i][p].first << "\n"; + std::cout << "land2.land[i][q].first : " << land2.land[i][q].first << "\n"; + //getchar(); + } + + if ( land1.land[i][p].first < land2.land[i][q].first ) + { + if ( operation_on_pair_of_landscapesDBG ) + { + std::cout << "first \n"; + std::cout << " function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) : "<< function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) << "\n"; + //std::cout << "oper( " << land1.land[i][p].second <<"," << function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) << " : " << oper( land1.land[i][p].second , function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) ) << "\n"; + } + lambda_n.push_back( + std::make_pair( + land1.land[i][p].first , + oper( (double)land1.land[i][p].second , function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) ) + ) + ); + ++p; + continue; + } + if ( land1.land[i][p].first > land2.land[i][q].first ) + { + if ( operation_on_pair_of_landscapesDBG ) + { + std::cout << "Second \n"; + std::cout << "function_value("<< land1.land[i][p-1].first << " " << land1.land[i][p-1].second <<" ,"<< land1.land[i][p].first << " " << land1.land[i][p].second <<", " << land2.land[i][q].first<<" ) : " << function_value( land1.land[i][p-1] , land1.land[i][p-1] ,land2.land[i][q].first ) << "\n"; + std::cout << "oper( " << function_value( land1.land[i][p] , land1.land[i][p-1] ,land2.land[i][q].first ) <<"," << land2.land[i][q].second <<" : " << oper( land2.land[i][q].second , function_value( land1.land[i][p] , land1.land[i][p-1] ,land2.land[i][q].first ) ) << "\n"; + } + lambda_n.push_back( std::make_pair( land2.land[i][q].first , oper( function_value( land1.land[i][p] , land1.land[i][p-1] ,land2.land[i][q].first ) , land2.land[i][q].second ) ) ); + ++q; + continue; + } + if ( land1.land[i][p].first == land2.land[i][q].first ) + { + if (operation_on_pair_of_landscapesDBG)std::cout << "Third \n"; + lambda_n.push_back( std::make_pair( land2.land[i][q].first , oper( land1.land[i][p].second , land2.land[i][q].second ) ) ); + ++p;++q; + } + if (operation_on_pair_of_landscapesDBG){std::cout << "Next iteration \n";} + } + while ( (p+1 < land1.land[i].size())&&(q+1 >= land2.land[i].size()) ) + { + if (operation_on_pair_of_landscapesDBG) + { + std::cout << "New point : " << land1.land[i][p].first << " oper(land1.land[i][p].second,0) : " << oper(land1.land[i][p].second,0) << std::endl; + } + lambda_n.push_back( std::make_pair(land1.land[i][p].first , oper(land1.land[i][p].second,0) ) ); + ++p; + } + while ( (p+1 >= land1.land[i].size())&&(q+1 < land2.land[i].size()) ) + { + if (operation_on_pair_of_landscapesDBG) + { + std::cout << "New point : " << land2.land[i][q].first << " oper(0,land2.land[i][q].second) : " << oper(0,land2.land[i][q].second) << std::endl; + } + lambda_n.push_back( std::make_pair(land2.land[i][q].first , oper(0,land2.land[i][q].second) ) ); + ++q; + } + lambda_n.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) ); + //CHANGE + //result.land[i] = lambda_n; + result.land[i].swap(lambda_n); + } + if ( land1.land.size() > std::min( land1.land.size() , land2.land.size() ) ) + { + if (operation_on_pair_of_landscapesDBG){std::cout << "land1.land.size() > std::min( land1.land.size() , land2.land.size() )" << std::endl;} + for ( size_t i = std::min( land1.land.size() , land2.land.size() ) ; i != std::max( land1.land.size() , land2.land.size() ) ; ++i ) + { + std::vector< std::pair<double,double> > lambda_n( land1.land[i] ); + for ( size_t nr = 0 ; nr != land1.land[i].size() ; ++nr ) + { + lambda_n[nr] = std::make_pair( land1.land[i][nr].first , oper( land1.land[i][nr].second , 0 ) ); + } + //CHANGE + //result.land[i] = lambda_n; + result.land[i].swap(lambda_n); + } + } + if ( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) ) + { + if (operation_on_pair_of_landscapesDBG){std::cout << "( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) ) " << std::endl;} + for ( size_t i = std::min( land1.land.size() , land2.land.size() ) ; i != std::max( land1.land.size() , land2.land.size() ) ; ++i ) + { + std::vector< std::pair<double,double> > lambda_n( land2.land[i] ); + for ( size_t nr = 0 ; nr != land2.land[i].size() ; ++nr ) + { + lambda_n[nr] = std::make_pair( land2.land[i][nr].first , oper( 0 , land2.land[i][nr].second ) ); + } + //CHANGE + //result.land[i] = lambda_n; + result.land[i].swap(lambda_n); + } + } + if ( operation_on_pair_of_landscapesDBG ){std::cout << "operation_on_pair_of_landscapes END\n";std::cin.ignore();} + return result; +}//operation_on_pair_of_landscapes + + + +double compute_maximal_distance_non_symmetric( const Persistence_landscape& pl1, const Persistence_landscape& pl2 ) +{ + bool dbg = false; + if (dbg)std::cerr << " compute_maximal_distance_non_symmetric \n"; + //this distance is not symmetric. It compute ONLY distance between inflection points of pl1 and pl2. + double maxDist = 0; + size_t minimalNumberOfLevels = std::min( pl1.land.size() , pl2.land.size() ); + for ( size_t level = 0 ; level != minimalNumberOfLevels ; ++ level ) + { + if (dbg) + { + std::cerr << "Level : " << level << std::endl; + std::cerr << "PL1 : \n"; + for ( size_t i = 0 ; i != pl1.land[level].size() ; ++i ) + { + std::cerr << "(" <<pl1.land[level][i].first << "," << pl1.land[level][i].second << ") \n"; + } + std::cerr << "PL2 : \n"; + for ( size_t i = 0 ; i != pl2.land[level].size() ; ++i ) + { + std::cerr << "(" <<pl2.land[level][i].first << "," << pl2.land[level][i].second << ") \n"; + } + std::cin.ignore(); + } + + int p2Count = 0; + for ( size_t i = 1 ; i != pl1.land[level].size()-1 ; ++i ) //w tym przypadku nie rozwarzam punktow w nieskocznosci + { + while ( true ) + { + if ( (pl1.land[level][i].first>=pl2.land[level][p2Count].first) && (pl1.land[level][i].first<=pl2.land[level][p2Count+1].first) )break; + p2Count++; + } + double val = fabs( function_value( pl2.land[level][p2Count] , pl2.land[level][p2Count+1] , pl1.land[level][i].first ) - pl1.land[level][i].second); + if ( maxDist <= val )maxDist = val; + + if (dbg) + { + std::cerr << pl1.land[level][i].first <<"in [" << pl2.land[level][p2Count].first << "," << pl2.land[level][p2Count+1].first <<"] \n"; + std::cerr << "pl1[level][i].second : " << pl1.land[level][i].second << std::endl; + std::cerr << "function_value( pl2[level][p2Count] , pl2[level][p2Count+1] , pl1[level][i].first ) : " << function_value( pl2.land[level][p2Count] , pl2.land[level][p2Count+1] , pl1.land[level][i].first ) << std::endl; + std::cerr << "val : " << val << std::endl; + std::cin.ignore(); + } + } + } + + if (dbg)std::cerr << "minimalNumberOfLevels : " << minimalNumberOfLevels << std::endl; + + if ( minimalNumberOfLevels < pl1.land.size() ) + { + for ( size_t level = minimalNumberOfLevels ; level != pl1.land.size() ; ++ level ) + { + for ( size_t i = 0 ; i != pl1.land[level].size() ; ++i ) + { + if (dbg)std::cerr << "pl1[level][i].second : " << pl1.land[level][i].second << std::endl; + if ( maxDist < pl1.land[level][i].second )maxDist = pl1.land[level][i].second; + } + } + } + return maxDist; +} + + + + +double compute_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , double p ) +{ + bool dbg = false; + //This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a time: + + //first-second : + Persistence_landscape lan = first-second; + + //| first-second |: + lan = lan.abs(); + + if ( dbg ){std::cerr << "Abs of difference ; " << lan << std::endl;getchar();} + + if ( p < std::numeric_limits<double>::max() ) + { + //\int_{- \infty}^{+\infty}| first-second |^p + double result; + if ( p != 1 ) + { + if ( dbg )std::cerr << "Power != 1, compute integral to the power p\n"; + result = lan.compute_integral_of_landscape( (double)p ); + } + else + { + if ( dbg )std::cerr << "Power = 1, compute integral \n"; + result = lan.compute_integral_of_landscape(); + } + //(\int_{- \infty}^{+\infty}| first-second |^p)^(1/p) + return pow( result , 1/(double)p ); + } + else + { + //p == infty + if ( dbg )std::cerr << "Power = infty, compute maximum \n"; + return lan.compute_maximum(); + } +} + +double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second ) +{ + return std::max( compute_maximal_distance_non_symmetric(first,second) , compute_maximal_distance_non_symmetric(second,first) ); +} + + +bool comparePairsForMerging( std::pair< double , unsigned > first , std::pair< double , unsigned > second ) +{ + return (first.first < second.first); +} + + + + +double compute_inner_product( const Persistence_landscape& l1 , const Persistence_landscape& l2 ) +{ + bool dbg = false; + double result = 0; + + for ( size_t level = 0 ; level != std::min( l1.size() , l2.size() ) ; ++level ) + { + if ( dbg ){std::cerr << "Computing inner product for a level : " << level << std::endl;getchar();} + if ( l1.land[level].size() * l2.land[level].size() == 0 )continue; + + //endpoints of the interval on which we will compute the inner product of two locally linear functions: + double x1 = -std::numeric_limits<int>::max(); + double x2; + if ( l1.land[level][1].first < l2.land[level][1].first ) + { + x2 = l1.land[level][1].first; + } + else + { + x2 = l2.land[level][1].first; + } + + //iterators for the landscapes l1 and l2 + size_t l1It = 0; + size_t l2It = 0; + + while ( (l1It < l1.land[level].size()-1) && (l2It < l2.land[level].size()-1) ) + { + //compute the value of a inner product on a interval [x1,x2] + + double a,b,c,d; + + if ( l1.land[level][l1It+1].first != l1.land[level][l1It].first ) + { + a = (l1.land[level][l1It+1].second - l1.land[level][l1It].second)/(l1.land[level][l1It+1].first - l1.land[level][l1It].first); + } + else + { + a = 0; + } + b = l1.land[level][l1It].second - a*l1.land[level][l1It].first; + if ( l2.land[level][l2It+1].first != l2.land[level][l2It].first ) + { + c = (l2.land[level][l2It+1].second - l2.land[level][l2It].second)/(l2.land[level][l2It+1].first - l2.land[level][l2It].first); + } + else + { + c = 0; + } + d = l2.land[level][l2It].second - c*l2.land[level][l2It].first; + + double contributionFromThisPart + = + (a*c*x2*x2*x2/3 + (a*d+b*c)*x2*x2/2 + b*d*x2) - (a*c*x1*x1*x1/3 + (a*d+b*c)*x1*x1/2 + b*d*x1); + + result += contributionFromThisPart; + + if ( dbg ) + { + std::cerr << "[l1.land[level][l1It].first,l1.land[level][l1It+1].first] : " << l1.land[level][l1It].first << " , " << l1.land[level][l1It+1].first << std::endl; + std::cerr << "[l2.land[level][l2It].first,l2.land[level][l2It+1].first] : " << l2.land[level][l2It].first << " , " << l2.land[level][l2It+1].first << std::endl; + std::cerr << "a : " << a << ", b : " << b << " , c: " << c << ", d : " << d << std::endl; + std::cerr << "x1 : " << x1 << " , x2 : " << x2 << std::endl; + std::cerr << "contributionFromThisPart : " << contributionFromThisPart << std::endl; + std::cerr << "result : " << result << std::endl; + getchar(); + } + + //we have two intervals in which functions are constant: + //[l1.land[level][l1It].first , l1.land[level][l1It+1].first] + //and + //[l2.land[level][l2It].first , l2.land[level][l2It+1].first] + //We also have an interval [x1,x2]. Since the intervals in the landscapes cover the whole R, then it is clear that x2 + //is either l1.land[level][l1It+1].first of l2.land[level][l2It+1].first or both. Lets test it. + if ( x2 == l1.land[level][l1It+1].first ) + { + if ( x2 == l2.land[level][l2It+1].first ) + { + //in this case, we increment both: + ++l2It; + if ( dbg ){std::cerr << "Incrementing both \n";} + } + else + { + if ( dbg ){std::cerr << "Incrementing first \n";} + } + ++l1It; + } + else + { + //in this case we increment l2It + ++l2It; + if ( dbg ){std::cerr << "Incrementing second \n";} + } + //Now, we shift x1 and x2: + x1 = x2; + if ( l1.land[level][l1It+1].first < l2.land[level][l2It+1].first ) + { + x2 = l1.land[level][l1It+1].first; + } + else + { + x2 = l2.land[level][l2It+1].first; + } + + } + + } + return result; +} + + +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. + std::ofstream out; + + std::ostringstream nameSS; + nameSS << filename << "_GnuplotScript"; + std::string nameStr = nameSS.str(); + out.open( nameStr ); + + if ( (xRangeBegin != std::numeric_limits<double>::max()) || (xRangeEnd != std::numeric_limits<double>::max()) || (yRangeBegin != std::numeric_limits<double>::max()) || (yRangeEnd != std::numeric_limits<double>::max()) ) + { + out << "set xrange [" << xRangeBegin << " : " << xRangeEnd << "]" << std::endl; + out << "set yrange [" << yRangeBegin << " : " << yRangeEnd << "]" << std::endl; + } + + if ( from == std::numeric_limits<int>::max() ){from = 0;} + if ( to == std::numeric_limits<int>::max() ){to = this->land.size();} + + out << "plot "; + for ( size_t lambda= std::min((size_t)from,this->land.size()) ; lambda != std::min((size_t)to,this->land.size()) ; ++lambda ) + { + //out << " '-' using 1:2 title 'l" << lambda << "' with lp"; + out << " '-' using 1:2 notitle with lp"; + if ( lambda+1 != std::min((size_t)to,this->land.size()) ) + { + out << ", \\"; + } + out << std::endl; + } + + for ( size_t lambda= std::min((size_t)from,this->land.size()) ; lambda != std::min((size_t)to,this->land.size()) ; ++lambda ) + { + for ( size_t i = 1 ; i != this->land[lambda].size()-1 ; ++i ) + { + out << this->land[lambda][i].first << " " << this->land[lambda][i].second << std::endl; + } + out << "EOF" << std::endl; + } + std::cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << std::endl; +} + + + + +}//namespace gudhi stat +}//namespace gudhi + + +#endif diff --git a/src/Gudhi_stat/include/gudhi/Persistence_landscape_on_grid.h b/src/Gudhi_stat/include/gudhi/Persistence_landscape_on_grid.h new file mode 100644 index 00000000..5703163a --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Persistence_landscape_on_grid.h @@ -0,0 +1,1563 @@ +/** 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 Persistence_landscape_on_grid_H_ +#define Persistence_landscape_on_grid_H_ + + +//standard include +#include <iostream> +#include <vector> +#include <limits> +#include <fstream> +#include <sstream> +#include <algorithm> +#include <cmath> +#include <limits> +#include <functional> +#include <functional> + + +//gudhi include + +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> + + + +namespace Gudhi +{ +namespace Persistence_representations +{ + +//predeclaration +class Persistence_landscape_on_grid; +template < typename operation > +Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid( const Persistence_landscape_on_grid& land1 , const Persistence_landscape_on_grid& land2 ); + +/** + * A clas implementing persistence landascpes by approximating them on a collection of grid points. * Persistence landscapes on grid allow vertorization, computations of distances, computations + * of projections to Real, computations of averages and scalar products. Therefore they implement suitable interfaces. + * It implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product + * Note that at the moment, due to roundoff errors during the construction of persistence landscapes on a grid, elements which are different by 0.000005 are considered the same. If the scale in your persistence diagrams + * is comparable to this value, please rescale them before use this code. +**/ + +//this class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product +class Persistence_landscape_on_grid +{ +public: + /** + * Default constructor. + **/ + Persistence_landscape_on_grid() + { + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); + this->grid_min = this->grid_max = 0; + } + + /** + * Constructor that takes as an input a vector of birth-death pairs. + **/ + Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ ); + + + /** + * Constructor that takes as an input a vector of birth-death pairs, parameters of the grid and number of landscape function to be created. + **/ + Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ , unsigned number_of_levels_of_landscape ); + + /** + * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed + * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. The additional parameters of this procedure are: ranges of grid, resoltion of a grid + * number of landscape functions to be created and the dimension of intervals that are need to be read from a file (in case of Gudhi format files). + **/ + Persistence_landscape_on_grid(const char* filename , double grid_min_, double grid_max_ , size_t number_of_points_ , unsigned number_of_levels_of_landscape , unsigned short dimension_ = std::numeric_limits<unsigned short>::max() ); + + /** + * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed + * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. The additional parameters of this procedure are: ranges of grid, resoltion of a grid + * and the dimension of intervals that are need to be read from a file (in case of Gudhi format files). + **/ + Persistence_landscape_on_grid(const char* filename , double grid_min_, double grid_max_ , size_t number_of_points_ , unsigned short dimension_ = std::numeric_limits<unsigned short>::max() ); + + + /** + * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed + * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. The additional parameter is the resoution of a grid and the number of landscape + * functions to be created. The remaning parameters are calculated based on data. + **/ + Persistence_landscape_on_grid(const char* filename , size_t number_of_points , unsigned number_of_levels_of_landscape , unsigned short dimension = std::numeric_limits<unsigned short>::max() ); + + /** + * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed + * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. The additional parameter is the resoution of a grid. The last parameter is the dimension + * of a persistence to read from the file. If your file contains only persistence pair in a single dimension, please set it up to std::numeric_limits<unsigned>::max(). + * The remaning parameters are calculated based on data. + **/ + Persistence_landscape_on_grid(const char* filename , size_t number_of_points , unsigned short dimension = std::numeric_limits<unsigned short>::max() ); + + + /** + * This procedure loads a landscape from file. It erase all the data that was previously stored in this landscape. + **/ + void load_landscape_from_file( const char* filename ); + + + /** + * The procedure stores a landscape to a file. The file can be later used by a procedure load_landscape_from_file. + **/ + void print_to_file( const char* filename )const; + + + + /** + * This function compute integral of the landscape (defined formally as sum of integrals on R of all landscape functions) + **/ + double compute_integral_of_landscape()const + { + size_t maximal_level = this->number_of_nonzero_levels(); + double result = 0; + for ( size_t i = 0 ; i != maximal_level ; ++i ) + { + result += this->compute_integral_of_landscape(i); + } + return result; + } + + + /** + * This function compute integral of the 'level'-level of a landscape. + **/ + double compute_integral_of_landscape( size_t level )const + { + bool dbg = false; + double result = 0; + double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1); + + if ( dbg ) + { + std::cerr << "this->grid_max : " << this->grid_max << std::endl; + std::cerr << "this->grid_min : " << this->grid_min << std::endl; + std::cerr << "this->values_of_landscapes.size() : " << this->values_of_landscapes.size() << std::endl; + getchar(); + } + + + double previous_x = this->grid_min-dx; + double previous_y = 0; + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + double current_x = previous_x + dx; + double current_y = 0; + if ( this->values_of_landscapes[i].size() > level )current_y = this->values_of_landscapes[i][level]; + + if ( dbg ) + { + std::cerr << "this->values_of_landscapes[i].size() : " << this->values_of_landscapes[i].size() << " , level : " << level << std::endl; + if ( this->values_of_landscapes[i].size() > level )std::cerr << "this->values_of_landscapes[i][level] : " << this->values_of_landscapes[i][level] << std::endl; + std::cerr << "previous_y : " << previous_y << std::endl; + std::cerr << "current_y : " << current_y << std::endl; + std::cerr << "dx : " << dx << std::endl; + std::cerr << "0.5*dx*( previous_y + current_y ); " << 0.5*dx*( previous_y + current_y ) << std::endl; + } + + result += 0.5*dx*( previous_y + current_y ); + previous_x = current_x; + previous_y = current_y; + } + return result; + } + + /** + * This function compute integral of the landscape p-th power of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions) + **/ + double compute_integral_of_landscape( double p )const + { + size_t maximal_level = this->number_of_nonzero_levels(); + double result = 0; + for ( size_t i = 0 ; i != maximal_level ; ++i ) + { + result += this->compute_integral_of_landscape(p,i); + } + return result; + } + + /** + * This function compute integral of the landscape p-th power of a level of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions) + **/ + double compute_integral_of_landscape( double p , size_t level )const + { + bool dbg = false; + + double result = 0; + double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1); + double previous_x = this->grid_min; + double previous_y = 0; + if ( this->values_of_landscapes[0].size() > level )previous_y = this->values_of_landscapes[0][level]; + + if ( dbg ) + { + std::cerr << "dx : " << dx << std::endl; + std::cerr << "previous_x : " << previous_x << std::endl; + std::cerr << "previous_y : " << previous_y << std::endl; + std::cerr << "power : " << p << std::endl; + getchar(); + } + + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + double current_x = previous_x + dx; + double current_y = 0; + if ( this->values_of_landscapes[i].size() > level )current_y = this->values_of_landscapes[i][level]; + + if ( dbg )std::cerr << "current_y : " << current_y << std::endl; + + if ( current_y == previous_y )continue; + + std::pair<double,double> coef = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y ) , std::make_pair( current_x , current_y ) ); + double a = coef.first; + double b = coef.second; + + if ( dbg ) + { + std::cerr << "A line passing through points : (" << previous_x << "," << previous_y << ") and (" << current_x << "," << current_y << ") is : " << a << "x+" << b << std::endl; + } + + //In this interval, the landscape has a form f(x) = ax+b. We want to compute integral of (ax+b)^p = 1/a * (ax+b)^{p+1}/(p+1) + double value_to_add = 0; + if ( a != 0 ) + { + value_to_add = 1/(a*(p+1)) * ( pow((a*current_x+b),p+1) - pow((a*previous_x+b),p+1)); + } + else + { + value_to_add = ( current_x - previous_x )*( pow(b,p) ); + } + result += value_to_add; + if ( dbg ) + { + std::cerr << "Increasing result by : " << value_to_add << std::endl; + std::cerr << "restult : " << result << std::endl; + getchar(); + } + previous_x = current_x; + previous_y = current_y; + } + if ( dbg )std::cerr << "The total result is : " << result << std::endl; + return result; + } + + /** + * Writing landscape into a stream. A i-th level landscape starts with a string "lambda_i". Then the discontinuity points of the landscapes follows. + * Shall those points be joined with lines, we will obtain the i-th landscape function. + **/ + friend std::ostream& operator<<(std::ostream& out, const Persistence_landscape_on_grid& land ) + { + double dx = (land.grid_max - land.grid_min)/(double)(land.values_of_landscapes.size()-1); + double x = land.grid_min; + for ( size_t i = 0 ; i != land.values_of_landscapes.size() ; ++i ) + { + out << x << " : "; + for ( size_t j = 0 ; j != land.values_of_landscapes[i].size() ; ++j ) + { + out << land.values_of_landscapes[i][j] << " "; + } + out << std::endl; + x += dx; + } + return out; + } + + template < typename oper > + friend Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid( const Persistence_landscape_on_grid& land1 , const Persistence_landscape_on_grid& land2 ); + + + /** + * A function that computes the value of a landscape at a given point. The parameters of the function are: unsigned level and double x. + * The procedure will compute the value of the level-landscape at the point x. + **/ + double compute_value_at_a_given_point( unsigned level , double x )const + { + bool dbg = false; + if ( (x < this->grid_min) || (x > this->grid_max) )return 0; + + //find a position of a vector closest to x: + double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1); + size_t position = size_t((x-this->grid_min)/dx); + + if ( dbg ) + { + std::cerr << "This is a procedure compute_value_at_a_given_point \n"; + std::cerr << "level : " << level << std::endl; + std::cerr << "x : " << x << std::endl; + std::cerr << "psoition : " << position << std::endl; + } + //check if we are not exacly in the grid point: + if ( almost_equal( position*dx+ this->grid_min , x) ) + { + if ( this->values_of_landscapes[position].size() < level ) + { + return this->values_of_landscapes[position][level]; + } + else + { + return 0; + } + } + //in the other case, approximate with a line: + std::pair<double,double> line; + if ( (this->values_of_landscapes[position].size() > level) && (this->values_of_landscapes[position+1].size() > level) ) + { + line = compute_parameters_of_a_line( std::make_pair( position*dx+ this->grid_min , this->values_of_landscapes[position][level] ) , std::make_pair( (position+1)*dx+ this->grid_min , this->values_of_landscapes[position+1][level] ) ); + } + else + { + if ( (this->values_of_landscapes[position].size() > level) || (this->values_of_landscapes[position+1].size() > level) ) + { + if ( (this->values_of_landscapes[position].size() > level) ) + { + line = compute_parameters_of_a_line( std::make_pair( position*dx+ this->grid_min , this->values_of_landscapes[position][level] ) , std::make_pair( (position+1)*dx+ this->grid_min , 0 ) ); + } + else + { + //(this->values_of_landscapes[position+1].size() > level) + line = compute_parameters_of_a_line( std::make_pair( position*dx+ this->grid_min , 0 ) , std::make_pair( (position+1)*dx+ this->grid_min , this->values_of_landscapes[position+1][level] ) ); + } + } + else + { + return 0; + } + } + //compute the value of the linear function parametrized by line on a point x: + return line.first*x+line.second; + } + + + public: + /** + *\private A function that compute sum of two landscapes. + **/ + friend Persistence_landscape_on_grid add_two_landscapes ( const Persistence_landscape_on_grid& land1 , const Persistence_landscape_on_grid& land2 ) + { + return operation_on_pair_of_landscapes_on_grid< std::plus<double> >(land1,land2); + } + + /** + *\private A function that compute difference of two landscapes. + **/ + friend Persistence_landscape_on_grid subtract_two_landscapes ( const Persistence_landscape_on_grid& land1 , const Persistence_landscape_on_grid& land2 ) + { + return operation_on_pair_of_landscapes_on_grid< std::minus<double> >(land1,land2); + } + + /** + * An operator +, that compute sum of two landscapes. + **/ + friend Persistence_landscape_on_grid operator+( const Persistence_landscape_on_grid& first , const Persistence_landscape_on_grid& second ) + { + return add_two_landscapes( first,second ); + } + + /** + * An operator -, that compute difference of two landscapes. + **/ + friend Persistence_landscape_on_grid operator-( const Persistence_landscape_on_grid& first , const Persistence_landscape_on_grid& second ) + { + return subtract_two_landscapes( first,second ); + } + + /** + * An operator * that allows multipilication of a landscape by a real number. + **/ + friend Persistence_landscape_on_grid operator*( const Persistence_landscape_on_grid& first , double con ) + { + return first.multiply_lanscape_by_real_number_not_overwrite(con); + } + + /** + * An operator * that allows multipilication of a landscape by a real number (order of parameters swapped). + **/ + friend Persistence_landscape_on_grid operator*( double con , const Persistence_landscape_on_grid& first ) + { + return first.multiply_lanscape_by_real_number_not_overwrite(con); + } + + friend bool check_if_defined_on_the_same_domain( const Persistence_landscape_on_grid& land1, const Persistence_landscape_on_grid& land2 ) + { + if ( land1.values_of_landscapes.size() != land2.values_of_landscapes.size() )return false; + if ( land1.grid_min != land2.grid_min )return false; + if ( land1.grid_max != land2.grid_max )return false; + return true; + } + + /** + * Operator +=. The second parameter is persistnece landwscape. + **/ + Persistence_landscape_on_grid operator += ( const Persistence_landscape_on_grid& rhs ) + { + *this = *this + rhs; + return *this; + } + + /** + * Operator -=. The second parameter is persistnece landwscape. + **/ + Persistence_landscape_on_grid operator -= ( const Persistence_landscape_on_grid& rhs ) + { + *this = *this - rhs; + return *this; + } + + + /** + * Operator *=. The second parameter is a real number by which the y values of all landscape functions are multiplied. The x-values remain unchanged. + **/ + Persistence_landscape_on_grid operator *= ( double x ) + { + *this = *this*x; + return *this; + } + + /** + * Operator /=. The second parameter is a real number. + **/ + Persistence_landscape_on_grid operator /= ( double x ) + { + if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." ); + *this = *this * (1/x); + return *this; + } + + /** + * An operator to compare two persistence landscapes. + **/ + bool operator == ( const Persistence_landscape_on_grid& rhs )const + { + bool dbg = true; + if ( this->values_of_landscapes.size() != rhs.values_of_landscapes.size() ) + { + if (dbg) std::cerr << "values_of_landscapes of incompatable sizes\n"; + return false; + } + if ( !almost_equal( this->grid_min , rhs.grid_min ) ) + { + if (dbg) std::cerr << "grid_min not equal\n"; + return false; + } + if ( !almost_equal(this->grid_max,rhs.grid_max ) ) + { + if (dbg) std::cerr << "grid_max not equal\n"; + return false; + } + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + for ( size_t aa = 0 ; aa != this->values_of_landscapes[i].size() ; ++aa ) + { + if ( !almost_equal( this->values_of_landscapes[i][aa] , rhs.values_of_landscapes[i][aa] ) ) + { + if (dbg) + { + std::cerr << "Problem in the position : " << i << " of values_of_landscapes. \n"; + std::cerr << this->values_of_landscapes[i][aa] << " " << rhs.values_of_landscapes[i][aa] << std::endl; + } + return false; + } + } + } + return true; + } + + + /** + * An operator to compare two persistence landscapes. + **/ + bool operator != ( const Persistence_landscape_on_grid& rhs )const + { + return !((*this) == rhs); + } + + + /** + * Computations of maximum (y) value of landscape. + **/ + double compute_maximum()const + { + //since the function can only be entirely positive or negative, the maximal value will be an extremal value in the arrays: + double max_value = -std::numeric_limits<double>::max(); + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + if ( this->values_of_landscapes[i].size() ) + { + if ( this->values_of_landscapes[i][0] > max_value )max_value = this->values_of_landscapes[i][0]; + if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] > max_value )max_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ]; + } + } + return max_value; + } + + /** + * Computations of minimum and maximum value of landscape. + **/ + std::pair<double,double> compute_minimum_maximum()const + { + //since the function can only be entirely positive or negative, the maximal value will be an extremal value in the arrays: + double max_value = -std::numeric_limits<double>::max(); + double min_value = 0; + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + if ( this->values_of_landscapes[i].size() ) + { + if ( this->values_of_landscapes[i][0] > max_value )max_value = this->values_of_landscapes[i][0]; + if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] > max_value )max_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ]; + + if ( this->values_of_landscapes[i][0] < min_value )min_value = this->values_of_landscapes[i][0]; + if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] < min_value )min_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ]; + } + } + return std::make_pair(min_value , max_value); + } + + /** + * 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 > get_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 > get_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 + { + size_t result = 0; + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + if ( this->values_of_landscapes[i].size() > result )result = this->values_of_landscapes[i].size(); + } + return result; + } + + /** + * Computations of a \f$L^i\f$ norm of landscape, where i is the input parameter. + **/ + 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); + + if ( i < std::numeric_limits<double>::max() ) + { + return compute_distance_of_landscapes_on_grid(*this,l,i); + } + else + { + return compute_max_norm_distance_of_landscapes(*this,l); + } + } + + /** + * An operator to compute the value of a landscape in the level 'level' at the argument 'x'. + **/ + double operator()(unsigned level,double x)const{return this->compute_value_at_a_given_point(level,x);} + + /** + * Computations of \f$L^{\infty}\f$ distance between two landscapes. + **/ + friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second ); + //friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second , unsigned& nrOfLand , double&x , double& y1, double& y2 ); + + + + + /** + * Function to compute absolute value of a PL function. The representation of persistence landscapes allow to store general PL-function. When computing distance betwen two landscapes, we compute difference between + * them. In this case, a general PL-function with negative value can appear as a result. Then in order to compute distance, we need to take its absolute value. This is the purpose of this procedure. + **/ + void abs() + { + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->values_of_landscapes[i].size() ; ++j ) + { + this->values_of_landscapes[i][j] = std::abs( this->values_of_landscapes[i][j] ); + } + } + } + + /** + * Computes the number of landscape functions. + **/ + size_t size()const{return this->number_of_nonzero_levels(); } + + /** + * Computate maximal value of lambda-level landscape. + **/ + double find_max( unsigned lambda )const + { + double max_value = -std::numeric_limits<double>::max(); + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + if ( this->values_of_landscapes[i].size() > lambda ) + { + if ( this->values_of_landscapes[i][lambda] > max_value )max_value = this->values_of_landscapes[i][lambda]; + } + } + return max_value; + } + + /** + * Function to compute inner (scalar) product of two landscapes. + **/ + friend double compute_inner_product( const Persistence_landscape_on_grid& l1 , const Persistence_landscape_on_grid& l2 ) + { + if ( !check_if_defined_on_the_same_domain(l1,l2) )throw "Landscapes are not defined on the same grid, the program will now terminate"; + size_t maximal_level = l1.number_of_nonzero_levels(); + double result = 0; + for ( size_t i = 0 ; i != maximal_level ; ++i ) + { + result += compute_inner_product(l1,l2,i); + } + return result; + } + + + + /** + * Function to compute inner (scalar) product of given levels of two landscapes. + **/ + friend double compute_inner_product( const Persistence_landscape_on_grid& l1 , const Persistence_landscape_on_grid& l2 , size_t level ) + { + bool dbg = false; + + if ( !check_if_defined_on_the_same_domain(l1,l2) )throw "Landscapes are not defined on the same grid, the program will now terminate"; + double result = 0; + + double dx = (l1.grid_max - l1.grid_min)/(double)(l1.values_of_landscapes.size()-1); + + double previous_x = l1.grid_min-dx; + double previous_y_l1 = 0; + double previous_y_l2 = 0; + for ( size_t i = 0 ; i != l1.values_of_landscapes.size() ; ++i ) + { + if ( dbg )std::cerr << "i : " << i << std::endl; + + double current_x = previous_x + dx; + double current_y_l1 = 0; + if ( l1.values_of_landscapes[i].size() > level )current_y_l1 = l1.values_of_landscapes[i][level]; + + double current_y_l2 = 0; + if ( l2.values_of_landscapes[i].size() > level )current_y_l2 = l2.values_of_landscapes[i][level]; + + if ( dbg ) + { + std::cerr << "previous_x : " << previous_x << std::endl; + std::cerr << "previous_y_l1 : " << previous_y_l1 << std::endl; + std::cerr << "current_y_l1 : " << current_y_l1 << std::endl; + std::cerr << "previous_y_l2 : " << previous_y_l2 << std::endl; + std::cerr << "current_y_l2 : " << current_y_l2 << std::endl; + } + + std::pair<double,double> l1_coords = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y_l1 ) , std::make_pair( current_x , current_y_l1 ) ); + std::pair<double,double> l2_coords = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y_l2 ) , std::make_pair( current_x , current_y_l2 ) ); + + //let us assume that the first line is of a form y = ax+b, and the second one is of a form y = cx + d. Then here are a,b,c,d: + double a = l1_coords.first; + double b = l1_coords.second; + + double c = l2_coords.first; + double d = l2_coords.second; + + if ( dbg ) + { + std::cerr << "Here are the formulas for a line: \n"; + std::cerr << "a : " << a << std::endl; + std::cerr << "b : " << b << std::endl; + std::cerr << "c : " << c << std::endl; + std::cerr << "d : " << d << std::endl; + } + + //now, to compute the inner product in this interval we need to compute the integral of (ax+b)(cx+d) = acx^2 + (ad+bc)x + bd in the interval from previous_x to current_x: + //The integal is ac/3*x^3 + (ac+bd)/2*x^2 + bd*x + + double added_value = (a*c/3*current_x*current_x*current_x + (a*d+b*c)/2*current_x*current_x + b*d*current_x)- + (a*c/3*previous_x*previous_x*previous_x + (a*d+b*c)/2*previous_x*previous_x + b*d*previous_x); + + if ( dbg ) + { + std::cerr << "Value of the integral on the left end ie : " << previous_x << " is : " << a*c/3*previous_x*previous_x*previous_x + (a*d+b*c)/2*previous_x*previous_x + b*d*previous_x << std::endl; + std::cerr << "Value of the integral on the right end i.e. : " << current_x << " is " << a*c/3*current_x*current_x*current_x + (a*d+b*c)/2*current_x*current_x + b*d*current_x << std::endl; + } + + result += added_value; + + if ( dbg ) + { + std::cerr << "added_value : " << added_value << std::endl; + std::cerr << "result : " << result << std::endl; + getchar(); + } + + + previous_x = current_x; + previous_y_l1 = current_y_l1; + previous_y_l2 = current_y_l2; + + } + return result; + } + + + /** + * Computations of \f$L^{p}\f$ distance between two landscapes on a grid. p is the parameter of the procedure. + * FIXME: Note that, due to the grid representation, the method below may give non--accurate results in case when the landscape P and Q the difference of which we want to compute + * are interxsecting. This is a consequence of a general way they are computed. In the future, an integral of absolute value of a difference of P and Q will be given as a separated + * function to fix that inaccuracy. + **/ + friend double compute_distance_of_landscapes_on_grid( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second , double p ) + { + bool dbg = false; + //This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a time: + + if ( dbg ) + { + std::cerr << "first : " << first << std::endl; + std::cerr << "second : " << second << std::endl; + getchar(); + } + + //first-second : + Persistence_landscape_on_grid lan = first-second; + + if ( dbg ) + { + std::cerr << "Difference : " << lan << std::endl; + } + + //| first-second |: + lan.abs(); + + if ( dbg ) + { + std::cerr << "Abs : " << lan << std::endl; + } + + if ( p < std::numeric_limits< double >::max() ) + { + //\int_{- \infty}^{+\infty}| first-second |^p + double result; + if ( p != 1 ) + { + if (dbg){std::cerr << "p : " << p << std::endl; getchar();} + result = lan.compute_integral_of_landscape( (double)p ); + if (dbg){std::cerr << "integral : " << result << std::endl;getchar();} + } + else + { + result = lan.compute_integral_of_landscape(); + if (dbg){std::cerr << "integral, wihtout power : " << result << std::endl;getchar();} + } + //(\int_{- \infty}^{+\infty}| first-second |^p)^(1/p) + return pow( result , 1/(double)p ); + } + else + { + //p == infty + return lan.compute_maximum(); + } + } + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + //Functions that are needed for that class to implement the concept. + + /** + * 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. + * At the moment this function is not tested, since it is quite likelly to be changed in the future. Given this, when using it, keep in mind that it + * will be most likelly changed in the next versions. + **/ + double project_to_R( int number_of_function )const + { + return this->compute_integral_of_landscape( (size_t)number_of_function ); + } + + /** + * 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()const + { + return number_of_functions_for_projections_to_reals; + } + + + + + /** + * 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 )const + { + //TODO, think of something smarter over here + if ( ( number_of_function < 0 ) || ( (size_t)number_of_function >= this->values_of_landscapes.size() ) ) + { + throw "Wrong number of function\n"; + } + std::vector<double> v( this->values_of_landscapes.size() ); + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + v[i] = 0; + if ( this->values_of_landscapes[i].size() > (size_t)number_of_function ) + { + v[ i ] = this->values_of_landscapes[i][number_of_function]; + } + } + return v; + } + + /** + * 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()const + { + return number_of_functions_for_vectorization; + } + + + + + + /** + * A function to compute averaged persistence landscape on a grid, based on vector of persistence landscapes on grid. + * This function is required by Topological_data_with_averages concept. + **/ + void compute_average( const std::vector< Persistence_landscape_on_grid* >& to_average ) + { + + bool dbg = false; + //After execution of this procedure, the average is supposed to be in the current object. To make sure that this is the case, we need to do some cleaning first. + this->values_of_landscapes .clear(); + this->grid_min = this->grid_max = 0; + + //if there is nothing to averate, then the average is a zero landscape. + if ( to_average.size() == 0 )return; + + //now we need to check if the grids in all objects of to_average are the same: + for ( size_t i = 0 ; i != to_average.size() ; ++i ) + { + if ( !check_if_defined_on_the_same_domain(*(to_average[0]),*(to_average[i])) )throw "Two grids are not compatible"; + } + + this->values_of_landscapes = std::vector< std::vector<double> >( (to_average[0])->values_of_landscapes.size() ); + this->grid_min = (to_average[0])->grid_min; + this->grid_max = (to_average[0])->grid_max; + + if ( dbg ) + { + std::cerr << "Computations of average. The data from the current landscape have been cleared. We are ready to do the computations. \n"; + } + + //for every point in the grid: + for ( size_t grid_point = 0 ; grid_point != (to_average[0])->values_of_landscapes.size() ; ++grid_point ) + { + + //set up a vector of the correct size: + size_t maximal_size_of_vector = 0; + for ( size_t land_no = 0 ; land_no != to_average.size() ; ++land_no ) + { + if ( (to_average[land_no])->values_of_landscapes[grid_point].size() > maximal_size_of_vector ) + maximal_size_of_vector = (to_average[land_no])->values_of_landscapes[grid_point].size(); + } + this->values_of_landscapes[grid_point] = std::vector<double>( maximal_size_of_vector ); + + if ( dbg ) + { + std::cerr << "We are considering the point : " << grid_point << " of the grid. In this point, there are at most : " << maximal_size_of_vector << " nonzero landscape functions \n"; + } + + //and compute an arythmetic average: + for ( size_t land_no = 0 ; land_no != to_average.size() ; ++land_no ) + { + //summing: + for ( size_t i = 0 ; i != (to_average[land_no])->values_of_landscapes[grid_point].size() ; ++i ) + { + //compute the average in a smarter way. + this->values_of_landscapes[grid_point][i] += (to_average[land_no])->values_of_landscapes[grid_point][i]; + } + } + //normalizing: + for ( size_t i = 0 ; i != this->values_of_landscapes[grid_point].size() ; ++i ) + { + this->values_of_landscapes[grid_point][i] /= (double)to_average.size(); + } + } + }//compute_average + + + /** + * A function to compute distance between persistence landscape on a grid. + * The parameter of this functionis a Persistence_landscape_on_grid. + * This function is required in Topological_data_with_distances concept. + * For max norm distance, set power to std::numeric_limits<double>::max() + **/ + double distance( const Persistence_landscape_on_grid& second , double power = 1 )const + { + if ( power < std::numeric_limits<double>::max() ) + { + return compute_distance_of_landscapes_on_grid( *this , second , power ); + } + else + { + return compute_max_norm_distance_of_landscapes( *this , second ); + } + } + + /** + * A function to compute scalar product of persistence landscape on a grid. + * The parameter of this functionis a Persistence_landscape_on_grid. + * This function is required in Topological_data_with_scalar_product concept. + **/ + double compute_scalar_product( const Persistence_landscape_on_grid& second ) + { + return compute_inner_product( (*this) , second ); + } + + //end of implementation of functions needed for concepts. + + + + + + + + + + + + + + + + + + + + + /** + * A function that returns values of landsapes. It can be used for vizualization + **/ + 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. Over here we need to specify which landscapes do we want to plot. + * In addition, the user may specify the range (min and max) where landscape is plot. The fefault values for min and max are std::numeric_limits<double>::max(). If the procedure detect those + * values, it will determine the range so that the whole landscape is supported there. If at least one min or max value is different from std::numeric_limits<double>::max(), then the values + * provided by the user will be used. + **/ + void plot( const char* filename , size_t from_ , size_t to_ )const + { + this->plot( filename , std::numeric_limits<double>::max() , std::numeric_limits<double>::max(), std::numeric_limits<double>::max() , std::numeric_limits<double>::max() , 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 = 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() , size_t from_ = std::numeric_limits<size_t>::max(), size_t to_= std::numeric_limits<size_t>::max() )const; + + +protected: + double grid_min; + double grid_max; + std::vector< std::vector< double > > values_of_landscapes; + size_t number_of_functions_for_vectorization; + size_t number_of_functions_for_projections_to_reals; + + void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals() + { + //warning, this function can be only called after filling in the values_of_landscapes vector. + this->number_of_functions_for_vectorization = this->values_of_landscapes.size(); + this->number_of_functions_for_projections_to_reals = this->values_of_landscapes.size(); + } + void set_up_values_of_landscapes( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ , unsigned number_of_levels = std::numeric_limits<unsigned>::max() ); + Persistence_landscape_on_grid multiply_lanscape_by_real_number_not_overwrite( double x )const; +}; + + +void Persistence_landscape_on_grid::set_up_values_of_landscapes( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_, unsigned number_of_levels ) +{ + bool dbg = false; + if ( dbg ) + { + std::cerr << "Here is the procedure : set_up_values_of_landscapes. The parameters are : grid_min_ : " << grid_min_ << ", grid_max_ : " << grid_max_ << ", number_of_points_ : " << number_of_points_ << ", number_of_levels: " << number_of_levels << std::endl; + std::cerr << "Here are the intervals at our disposal : \n"; + for ( size_t i = 0 ; i != p.size() ; ++i ) + { + std::cerr << p[i].first << " , " << p[i].second << std::endl; + } + } + + if ( (grid_min_ == std::numeric_limits<double>::max()) || (grid_max_ == std::numeric_limits<double>::max()) ) + { + //in this case, we need to find grid_min_ and grid_min_ based on the data. + double min = std::numeric_limits<double>::max(); + double max = std::numeric_limits<double>::min(); + for ( size_t i = 0 ; i != p.size() ; ++i ) + { + if ( p[i].first < min )min = p[i].first; + if ( p[i].second > max )max = p[i].second; + } + if ( grid_min_ == std::numeric_limits<double>::max() ) + { + grid_min_ = min; + } + else + { + //in this case grid_max_ == std::numeric_limits<double>::max() + grid_max_ = max; + } + } + + //if number_of_levels == std::numeric_limits<size_t>::max(), then we will have all the nonzero values of landscapes, and will store them in a vector + //if number_of_levels != std::numeric_limits<size_t>::max(), then we will use those vectors as heaps. + this->values_of_landscapes = std::vector< std::vector< double > >( number_of_points_+1 ); + + + this->grid_min = grid_min_; + this->grid_max = grid_max_; + + if ( grid_max_ <= grid_min_ ) + { + throw "Wrong parameters of grid_min and grid_max given to the procedure. THe grid have negative, or zero size. The program will now terminate.\n"; + } + + double dx = ( grid_max_ - grid_min_ )/(double)(number_of_points_); + //for every interval in the diagram: + for ( size_t int_no = 0 ; int_no != p.size() ; ++int_no ) + { + size_t grid_interval_begin = (p[int_no].first-grid_min_)/dx; + size_t grid_interval_end = (p[int_no].second-grid_min_)/dx; + size_t grid_interval_midpoint = (size_t)(0.5*(grid_interval_begin+grid_interval_end)); + + if ( dbg ) + { + std::cerr << "Considering an interval : " << p[int_no].first << "," << p[int_no].second << std::endl; + + std::cerr << "grid_interval_begin : " << grid_interval_begin << std::endl; + std::cerr << "grid_interval_end : " << grid_interval_end << std::endl; + std::cerr << "grid_interval_midpoint : " << grid_interval_midpoint << std::endl; + } + + double landscape_value = dx; + for ( size_t i = grid_interval_begin+1 ; i < grid_interval_midpoint ; ++i ) + { + if ( dbg ) + { + std::cerr << "Adding landscape value (going up) for a point : " << i << " equal : " << landscape_value << std::endl; + } + if ( number_of_levels != std::numeric_limits<unsigned>::max() ) + { + //we have a heap of no more that number_of_levels values. + //Note that if we are using heaps, we want to know the shortest distance in the heap. + //This is achieved by putting -distance to the heap. + if ( this->values_of_landscapes[i].size() >= number_of_levels ) + { + //in this case, the full heap is build, and we need to check if the landscape_value is not larger than the smallest element in the heap. + if ( -landscape_value < this->values_of_landscapes[i].front() ) + { + //if it is, we remove the largest value in the heap, and move on. + std::pop_heap (this->values_of_landscapes[i].begin(),this->values_of_landscapes[i].end()); + this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] = -landscape_value; + std::push_heap (this->values_of_landscapes[i].begin(),this->values_of_landscapes[i].end()); + } + } + else + { + + //in this case we are still filling in the array. + this->values_of_landscapes[i].push_back( -landscape_value ); + if ( this->values_of_landscapes[i].size() == number_of_levels-1 ) + { + //this->values_of_landscapes[i].size() == number_of_levels + //in this case we need to create the heep. + std::make_heap (this->values_of_landscapes[i].begin(),this->values_of_landscapes[i].end()); + } + } + } + else + { + //we have vector of all values + this->values_of_landscapes[i].push_back( landscape_value ); + } + landscape_value += dx; + } + for ( size_t i = grid_interval_midpoint ; i <= grid_interval_end ; ++i ) + { + if ( landscape_value > 0 ) + { + if ( number_of_levels != std::numeric_limits< unsigned >::max() ) + { + //we have a heap of no more that number_of_levels values + if ( this->values_of_landscapes[i].size() >= number_of_levels ) + { + //in this case, the full heap is build, and we need to check if the landscape_value is not larger than the smallest element in the heap. + if ( -landscape_value < this->values_of_landscapes[i].front() ) + { + //if it is, we remove the largest value in the heap, and move on. + std::pop_heap (this->values_of_landscapes[i].begin(),this->values_of_landscapes[i].end()); + this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] = -landscape_value; + std::push_heap (this->values_of_landscapes[i].begin(),this->values_of_landscapes[i].end()); + } + } + else + { + + //in this case we are still filling in the array. + this->values_of_landscapes[i].push_back( -landscape_value ); + if ( this->values_of_landscapes[i].size() == number_of_levels-1 ) + { + //this->values_of_landscapes[i].size() == number_of_levels + //in this case we need to create the heep. + std::make_heap (this->values_of_landscapes[i].begin(),this->values_of_landscapes[i].end()); + } + } + } + else + { + this->values_of_landscapes[i].push_back( landscape_value ); + } + + + if ( dbg ) + { + std::cerr << "AAdding landscape value (going down) for a point : " << i << " equal : " << landscape_value << std::endl; + } + } + landscape_value -= dx; + } + } + + if ( number_of_levels != std::numeric_limits< unsigned >::max() ) + { + //in this case, vectors are used as heaps. And, since we want to have the smallest element at the top of + //each heap, we store mminus distances. To get if right at the end, we need to multiply each value + //in the heap by -1 to get real vector of distances. + for ( size_t pt = 0 ; pt != this->values_of_landscapes.size() ; ++pt ) + { + //std::cerr << this->values_of_landscapes[pt].size() <<std::endl; + for ( size_t j = 0 ; j != this->values_of_landscapes[pt].size() ; ++j ) + { + this->values_of_landscapes[pt][j] *= -1; + } + } + } + + //and now we need to sort the values: + for ( size_t pt = 0 ; pt != this->values_of_landscapes.size() ; ++pt ) + { + std::sort( this->values_of_landscapes[pt].begin() , this->values_of_landscapes[pt].end() , std::greater<double>() ); + } +}//set_up_values_of_landscapes + + + +Persistence_landscape_on_grid::Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ ) +{ + this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ ); +}//Persistence_landscape_on_grid + + +Persistence_landscape_on_grid::Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ , unsigned number_of_levels_of_landscape ) +{ + this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ , number_of_levels_of_landscape ); +} + + +Persistence_landscape_on_grid::Persistence_landscape_on_grid(const char* filename , double grid_min_, double grid_max_ , size_t number_of_points_ , unsigned short dimension ) +{ + std::vector< std::pair< double , double > > p; + if ( dimension == std::numeric_limits<unsigned short>::max() ) + { + p = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + p = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ ); +} + +Persistence_landscape_on_grid::Persistence_landscape_on_grid(const char* filename , double grid_min_, double grid_max_ , size_t number_of_points_, unsigned number_of_levels_of_landscape, unsigned short dimension ) +{ + std::vector< std::pair< double , double > > p; + if ( dimension == std::numeric_limits<unsigned short>::max() ) + { + p = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + p = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ , number_of_levels_of_landscape ); +} + +Persistence_landscape_on_grid::Persistence_landscape_on_grid(const char* filename , size_t number_of_points_ , unsigned short dimension ) +{ + std::vector< std::pair< double , double > > p; + if ( dimension == std::numeric_limits<unsigned short>::max() ) + { + p = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + p = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + double grid_min_ = std::numeric_limits<double>::max(); + double grid_max_ = -std::numeric_limits<double>::max(); + for ( size_t i = 0 ; i != p.size() ; ++i ) + { + if ( p[i].first < grid_min_ )grid_min_ = p[i].first; + if ( p[i].second > grid_max_ )grid_max_ = p[i].second; + } + this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ ); +} + +Persistence_landscape_on_grid::Persistence_landscape_on_grid(const char* filename , size_t number_of_points_ , unsigned number_of_levels_of_landscape , unsigned short dimension ) +{ + std::vector< std::pair< double , double > > p; + if ( dimension == std::numeric_limits<unsigned short>::max() ) + { + p = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + p = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + double grid_min_ = std::numeric_limits<double>::max(); + double grid_max_ = -std::numeric_limits<double>::max(); + for ( size_t i = 0 ; i != p.size() ; ++i ) + { + if ( p[i].first < grid_min_ )grid_min_ = p[i].first; + if ( p[i].second > grid_max_ )grid_max_ = p[i].second; + } + this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ , number_of_levels_of_landscape ); +} + + +void Persistence_landscape_on_grid::load_landscape_from_file( const char* filename ) +{ + std::ifstream in; + in.open( filename ); + //check if the file exist. + if ( !in.good() ) + { + 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"; + } + + size_t number_of_points_in_the_grid = 0; + in >> this->grid_min >> this->grid_max >> number_of_points_in_the_grid; + + std::vector< std::vector< double > > v(number_of_points_in_the_grid); + std::string line; + std::getline(in, line); + double number; + for ( size_t i = 0 ; i != number_of_points_in_the_grid ; ++i ) + { + //read a line of a file and convert it to a vector. + std::vector< double > vv; + std::getline(in, line); + //std::cerr << "Reading line : " << line << std::endl;getchar(); + std::istringstream stream(line); + while (stream >> number) + { + vv.push_back(number); + } + v[i] = vv; + } + this->values_of_landscapes = v; + in.close(); +} + +void Persistence_landscape_on_grid::print_to_file( const char* filename )const +{ + std::ofstream out; + out.open( filename ); + + //first we store the parameters of the grid: + out << grid_min << std::endl << grid_max << std::endl << this->values_of_landscapes.size() << std::endl; + + //and now in the following lines, the values of this->values_of_landscapes for the following arguments: + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + for ( size_t j = 0 ; j != this->values_of_landscapes[i].size() ; ++j ) + { + out << this->values_of_landscapes[i][j] << " "; + } + out << std::endl; + } + + out.close(); +} + +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. + std::ofstream out; + + std::ostringstream nameSS; + nameSS << filename << "_GnuplotScript"; + std::string nameStr = nameSS.str(); + out.open( nameStr ); + + if ( min_x == max_x ) + { + std::pair<double,double> min_max = compute_minimum_maximum(); + out << "set xrange [" << this->grid_min << " : " << this->grid_max << "]" << std::endl; + out << "set yrange [" << min_max.first << " : " << min_max.second << "]" << std::endl; + } + else + { + out << "set xrange [" << min_x << " : " << max_x << "]" << std::endl; + out << "set yrange [" << min_y << " : " << max_y << "]" << std::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); + + + size_t from = 0; + if ( from_ != std::numeric_limits<size_t>::max() ) + { + if ( from_ < number_of_nonzero_levels ) + { + from = from_; + } + else + { + return; + } + } + size_t to = number_of_nonzero_levels; + if ( to_ != std::numeric_limits<size_t>::max() ) + { + if ( to_ < number_of_nonzero_levels ) + { + to = to_; + } + } + + + out << "plot "; + for ( size_t lambda= from ; lambda != to ; ++lambda ) + { + //out << " '-' using 1:2 title 'l" << lambda << "' with lp"; + out << " '-' using 1:2 notitle with lp"; + if ( lambda+1 != to ) + { + out << ", \\"; + } + out << std::endl; + } + + for ( size_t lambda = from ; lambda != to ; ++lambda ) + { + double point = this->grid_min; + for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i ) + { + double value = 0; + if ( this->values_of_landscapes[i].size() > lambda ) + { + value = this->values_of_landscapes[i][lambda]; + } + out << point << " " << value << std::endl; + point += dx; + } + out << "EOF" << std::endl; + } + std::cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << std::endl; +} + +template < typename T > +Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid ( const Persistence_landscape_on_grid& land1 , const Persistence_landscape_on_grid& land2 ) +{ + //first we need to check if the domains are the same: + if ( !check_if_defined_on_the_same_domain(land1,land2) )throw "Two grids are not compatible"; + + T oper; + Persistence_landscape_on_grid result; + result.values_of_landscapes = std::vector< std::vector< double > >( land1.values_of_landscapes.size() ); + result.grid_min = land1.grid_min; + result.grid_max = land1.grid_max; + + //now we perorm the operations: + for ( size_t grid_point = 0 ; grid_point != land1.values_of_landscapes.size() ; ++grid_point ) + { + result.values_of_landscapes[grid_point] = std::vector< double >( std::max( land1.values_of_landscapes[grid_point].size() , land2.values_of_landscapes[grid_point].size() ) ); + for ( size_t lambda = 0 ; lambda != std::max( land1.values_of_landscapes[grid_point].size() , land2.values_of_landscapes[grid_point].size() ) ; ++lambda ) + { + double value1 = 0; + double value2 = 0; + if ( lambda < land1.values_of_landscapes[grid_point].size() )value1 = land1.values_of_landscapes[grid_point][lambda]; + if ( lambda < land2.values_of_landscapes[grid_point].size() )value2 = land2.values_of_landscapes[grid_point][lambda]; + result.values_of_landscapes[grid_point][lambda] = oper( value1 , value2 ); + } + } + + return result; +} + +Persistence_landscape_on_grid Persistence_landscape_on_grid::multiply_lanscape_by_real_number_not_overwrite( double x )const +{ + Persistence_landscape_on_grid result; + result.values_of_landscapes = std::vector< std::vector< double > >( this->values_of_landscapes.size() ); + result.grid_min = this->grid_min; + result.grid_max = this->grid_max; + + for ( size_t grid_point = 0 ; grid_point != this->values_of_landscapes.size() ; ++grid_point ) + { + result.values_of_landscapes[grid_point] = std::vector< double >( this->values_of_landscapes[grid_point].size() ); + for ( size_t i = 0 ; i != this->values_of_landscapes[grid_point].size() ; ++i ) + { + result.values_of_landscapes[grid_point][i] = x*this->values_of_landscapes[grid_point][i]; + } + } + + return result; +} + +double compute_max_norm_distance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second ) +{ + double result = 0; + + //first we need to check if first and second is defined on the same domain" + if ( !check_if_defined_on_the_same_domain(first, second) )throw "Two grids are not compatible"; + + for ( size_t i = 0 ; i != first.values_of_landscapes.size() ; ++i ) + { + for ( size_t j = 0 ; j != std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) ; ++j ) + { + if ( result < abs( first.values_of_landscapes[i][j] - second.values_of_landscapes[i][j] ) ) + { + result = abs( first.values_of_landscapes[i][j] - second.values_of_landscapes[i][j] ); + } + } + if ( first.values_of_landscapes[i].size() == std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) ) + { + for ( size_t j = first.values_of_landscapes[i].size() ; j != second.values_of_landscapes[i].size() ; ++j ) + { + if ( result < second.values_of_landscapes[i][j] )result = second.values_of_landscapes[i][j]; + } + } + if ( second.values_of_landscapes[i].size() == std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) ) + { + for ( size_t j = second.values_of_landscapes[i].size() ; j != first.values_of_landscapes[i].size() ; ++j ) + { + if ( result < first.values_of_landscapes[i][j] )result = first.values_of_landscapes[i][j]; + } + } + } + return result; +} + + + +}//namespace Gudhi_stat +}//namespace Gudhi + +#endif diff --git a/src/Gudhi_stat/include/gudhi/Persistence_vectors.h b/src/Gudhi_stat/include/gudhi/Persistence_vectors.h new file mode 100644 index 00000000..7fde3413 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/Persistence_vectors.h @@ -0,0 +1,778 @@ +/* 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 Vector_distances_in_diagram_H +#define Vector_distances_in_diagram_H + +#include <fstream> +#include <cmath> +#include <algorithm> +#include <iostream> +#include <limits> +#include <functional> + +//gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/distance_functions.h> + + +namespace Gudhi +{ +namespace Persistence_representations +{ + +/* +template <typename T> +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 +{ + double operator() ( const std::pair< T,T >& f , const std::pair<T,T>& s ) + { + return std::max( fabs( f.first - s.first ) , fabs( f.second - s.second ) ); + } +}; + + + + +/** +* This is an implementation of idea presented in the paper 'Stable Topological Signatures for Points on 3D Shapes' by +* M. Carriere, S. Y. Oudot and M. Ovsjanikov published in Computer Graphics Forum (proc. SGP 2015). +* The parameter of the class is the class that computes distance used to construct the vectors. The typical function is +* either Eucludean of maximum (Manhattan) distance. +* This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product +* +**/ + +template <typename F> +class Vector_distances_in_diagram +{ +public: + /** + * The default constructor. + **/ + Vector_distances_in_diagram(){}; + + /** + * The constructor that takes as an input a multiset of persistence intervals (given as vector of birth-death pairs). The second parameter is the desiered length of the output vectors. + **/ + Vector_distances_in_diagram( const std::vector< std::pair< double , double > >& intervals , size_t where_to_cut ); + + /** + * The constructor taking as an input a file with birth-death pairs. The second parameter is the desiered length of the output vectors. + **/ + Vector_distances_in_diagram( const char* filename , size_t where_to_cut , unsigned dimension = std::numeric_limits<unsigned>::max() ); + + + /** + * Writing to a stream. + **/ + template <typename K> + friend std::ostream& operator << ( std::ostream& out , const Vector_distances_in_diagram<K>& d ) + { + for ( size_t i = 0 ; i != std::min( d.sorted_vector_of_distances.size() , d.where_to_cut) ; ++i ) + { + out << d.sorted_vector_of_distances[i] << " "; + } + return out; + } + + /** + * This procedure gives the value of a vector on a given 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]; + } + + /** + * Return a size of a vector. + **/ + inline size_t size()const{return this->sorted_vector_of_distances.size();} + + /** + * Write a vector to a file. + **/ + void write_to_file( const char* filename )const; + + /** + * Write a vector to a file. + **/ + void print_to_file( const char* filename )const + { + this->write_to_file( filename ); + } + + /** + * Loading a vector to a file. + **/ + void load_from_file( const char* filename ); + + /** + * Comparision operators: + **/ + 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 ) + { + if ( !almost_equal(this->sorted_vector_of_distances[i] , second.sorted_vector_of_distances[i]) )return false; + } + return true; + } + + bool operator != ( const Vector_distances_in_diagram& second )const + { + return !( *this == second ); + } + + + + + //Implementations of functions for various concepts. + /** + * Compute projection to real numbers of persistence vector. This function is required by the Real_valued_topological_data concept + * At the moment this function is not tested, since it is quite likelly to be changed in the future. Given this, when using it, keep in mind that it + * will be most likelly changed in the next versions. + **/ + 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()const + { + return this->number_of_functions_for_projections_to_reals; + } + + /** + * Compute a vectorization of a persistent vectors. It is required in a concept Vectorized_topological_data. + **/ + 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()const + { + return this->number_of_functions_for_vectorization; + } + + /** + * Compute a average of two persistent vectors. This function is required by Topological_data_with_averages concept. + **/ + void compute_average( const std::vector< Vector_distances_in_diagram* >& to_average ); + + /** + * Compute a distance of two persistent vectors. This function is required in Topological_data_with_distances concept. + * For max norm distance, set power to std::numeric_limits<double>::max() + **/ + 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 )const; + //end of implementation of functions needed for concepts. + + + + + + + + + + + + + + + + + + + /** + * For visualization use output from vectorize and build histograms. + **/ + std::vector< double > output_for_visualization()const + { + return this->sorted_vector_of_distances; + } + + + /** + * Create a gnuplot script to vizualize the data structure. + **/ + void plot( const char* filename )const + { + std::stringstream gnuplot_script; + gnuplot_script << filename << "_GnuplotScript"; + std::ofstream out; + out.open( gnuplot_script.str().c_str() ); + out << "set style data histogram" << std::endl; + out << "set style histogram cluster gap 1" << std::endl; + out << "set style fill solid border -1" << std::endl; + out << "plot '-' notitle" << std::endl; + for ( size_t i = 0 ; i != this->sorted_vector_of_distances.size() ; ++i ) + { + out << this->sorted_vector_of_distances[i] << std::endl; + } + out <<std::endl; + 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 > get_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 > get_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); + } + + //arythmetic operations: + template < typename Operation_type > + friend Vector_distances_in_diagram operation_on_pair_of_vectors( const Vector_distances_in_diagram& first , const Vector_distances_in_diagram& second , Operation_type opertion ) + { + Vector_distances_in_diagram result; + //Operation_type operation; + result.sorted_vector_of_distances.reserve(std::max( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ); + for ( size_t i = 0 ; i != std::min( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ; ++i ) + { + result.sorted_vector_of_distances.push_back( opertion( first.sorted_vector_of_distances[i] , second.sorted_vector_of_distances[i]) ); + } + if ( first.sorted_vector_of_distances.size() == std::min( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ) + { + for ( size_t i = std::min( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ; + i != std::max( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ; ++i ) + { + result.sorted_vector_of_distances.push_back( opertion(0,second.sorted_vector_of_distances[i]) ); + } + } + else + { + for ( size_t i = std::min( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ; + i != std::max( first.sorted_vector_of_distances.size() , second.sorted_vector_of_distances.size() ) ; ++i ) + { + result.sorted_vector_of_distances.push_back( opertion(first.sorted_vector_of_distances[i],0) ); + } + } + return result; + }//operation_on_pair_of_vectors + + /** + * This function implements an operation of multiplying Vector_distances_in_diagram by a scalar. + **/ + Vector_distances_in_diagram multiply_by_scalar( double scalar )const + { + Vector_distances_in_diagram result; + result.sorted_vector_of_distances.reserve( this->sorted_vector_of_distances.size() ); + for ( size_t i = 0 ; i != this->sorted_vector_of_distances.size() ; ++i ) + { + result.sorted_vector_of_distances.push_back( scalar * this->sorted_vector_of_distances[i] ); + } + return result; + }//multiply_by_scalar + + + + /** + * This function computes a sum of two objects of a type Vector_distances_in_diagram. + **/ + friend Vector_distances_in_diagram operator+( const Vector_distances_in_diagram& first , const Vector_distances_in_diagram& second ) + { + return operation_on_pair_of_vectors( first , second , std::plus<double>() ); + } + /** + * This function computes a difference of two objects of a type Vector_distances_in_diagram. + **/ + friend Vector_distances_in_diagram operator-( const Vector_distances_in_diagram& first , const Vector_distances_in_diagram& second ) + { + return operation_on_pair_of_vectors( first , second , std::minus<double>() ); + } + /** + * This function computes a product of an object of a type Vector_distances_in_diagram with real number. + **/ + friend Vector_distances_in_diagram operator*( double scalar , const Vector_distances_in_diagram& A ) + { + return A.multiply_by_scalar( scalar ); + } + /** + * This function computes a product of an object of a type Vector_distances_in_diagram with real number. + **/ + friend Vector_distances_in_diagram operator*( const Vector_distances_in_diagram& A , double scalar ) + { + return A.multiply_by_scalar( scalar ); + } + /** + * This function computes a product of an object of a type Vector_distances_in_diagram with real number. + **/ + Vector_distances_in_diagram operator*( double scalar ) + { + return this->multiply_by_scalar( scalar ); + } + /** + * += operator for Vector_distances_in_diagram. + **/ + Vector_distances_in_diagram operator += ( const Vector_distances_in_diagram& rhs ) + { + *this = *this + rhs; + return *this; + } + /** + * -= operator for Vector_distances_in_diagram. + **/ + Vector_distances_in_diagram operator -= ( const Vector_distances_in_diagram& rhs ) + { + *this = *this - rhs; + return *this; + } + /** + * *= operator for Vector_distances_in_diagram. + **/ + Vector_distances_in_diagram operator *= ( double x ) + { + *this = *this*x; + return *this; + } + /** + * /= operator for Vector_distances_in_diagram. + **/ + Vector_distances_in_diagram operator /= ( double x ) + { + if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." ); + *this = *this * (1/x); + return *this; + } + + +private: + std::vector< std::pair< double , double > > intervals; + std::vector< double > sorted_vector_of_distances; + size_t number_of_functions_for_vectorization; + size_t number_of_functions_for_projections_to_reals; + size_t where_to_cut; + + void compute_sorted_vector_of_distances_via_heap( size_t where_to_cut ); + void compute_sorted_vector_of_distances_via_vector_sorting( size_t where_to_cut ); + + Vector_distances_in_diagram( const std::vector< double >& sorted_vector_of_distances_ ):sorted_vector_of_distances(sorted_vector_of_distances_) + { + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); + } + + void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals() + { + //warning, this function can be only called after filling in the intervals vector. + this->number_of_functions_for_vectorization = this->sorted_vector_of_distances.size(); + this->number_of_functions_for_projections_to_reals = this->sorted_vector_of_distances.size(); + } +}; + + +template <typename F> +Vector_distances_in_diagram<F>::Vector_distances_in_diagram( const std::vector< std::pair< double,double > >& intervals_ , size_t where_to_cut_ ):where_to_cut(where_to_cut_) +{ + std::vector< std::pair< double,double > > i( intervals_ ); + this->intervals = i; + //this->compute_sorted_vector_of_distances_via_heap( where_to_cut ); + this->compute_sorted_vector_of_distances_via_vector_sorting(where_to_cut); + this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); +} + +template <typename F> +Vector_distances_in_diagram<F>::Vector_distances_in_diagram( const char* filename , size_t where_to_cut , unsigned dimension ):where_to_cut(where_to_cut) +{ + std::vector< std::pair< double , double > > intervals; + if ( dimension == std::numeric_limits<unsigned>::max() ) + { + intervals = read_persistence_intervals_in_one_dimension_from_file( filename ); + } + else + { + intervals = read_persistence_intervals_in_one_dimension_from_file( filename , dimension ); + } + this->intervals = intervals; + this->compute_sorted_vector_of_distances_via_heap( where_to_cut ); + //this->compute_sorted_vector_of_distances_via_vector_sorting( where_to_cut ); + set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); +} + +template < typename F> +void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distances_via_heap( size_t where_to_cut ) +{ + + bool dbg = false; + if ( dbg ) + { + std::cerr << "Here are the intervals : \n"; + for ( size_t i = 0 ; i != this->intervals.size() ; ++i ) + { + std::cerr << this->intervals[i].first << " , " << this->intervals[i].second <<std::endl; + } + } + where_to_cut = std::min(where_to_cut , (size_t)(0.5 * this->intervals.size() * ( this->intervals.size() - 1 ) + this->intervals.size())); + + std::vector< double > heap( where_to_cut , std::numeric_limits<int>::max() ); + std::make_heap (heap.begin(),heap.end()); + F f; + + //for every pair of points in the diagram, compute the minimum of their distance, and distance of those points from diagonal + for ( size_t i = 0 ; i < this->intervals.size() ; ++i ) + { + for ( size_t j = i+1 ; j < this->intervals.size() ; ++j ) + { + double value = std::min( + f( this->intervals[i] , this->intervals[j] ), + std::min( + f( this->intervals[i] , std::make_pair( 0.5*(this->intervals[i].first+this->intervals[i].second) , 0.5*(this->intervals[i].first+this->intervals[i].second) ) ), + f( this->intervals[j] , std::make_pair( 0.5*(this->intervals[j].first+this->intervals[j].second) , 0.5*(this->intervals[j].first+this->intervals[j].second) ) ) + ) + ); + + if ( dbg ) + { + std::cerr << "Value : " << value <<std::endl; + std::cerr << "heap.front() : " << heap.front() <<std::endl; + getchar(); + } + + if ( -value < heap.front() ) + { + if ( dbg ){std::cerr << "Replacing : " << heap.front() << " with : " << -value <<std::endl;getchar();} + //remove the first element from the heap + std::pop_heap (heap.begin(),heap.end()); + //heap.pop_back(); + //and put value there instead: + //heap.push_back(-value); + heap[ where_to_cut-1 ] = -value; + std::push_heap (heap.begin(),heap.end()); + } + } + } + + //now add distances of all points from diagonal + for ( size_t i = 0 ; i < this->intervals.size() ; ++i ) + { + double value = f( this->intervals[i] , std::make_pair( 0.5*(this->intervals[i].first+this->intervals[i].second) , 0.5*(this->intervals[i].first+this->intervals[i].second) ) ); + if ( -value < heap.front() ) + { + //std::cerr << "Replacing : " << heap.front() << " with : " << -value <<std::endl;getchar(); + //remove the first element from the heap + std::pop_heap (heap.begin(),heap.end()); + //heap.pop_back(); + //and put value there instead: + //heap.push_back(-value); + heap[ where_to_cut-1 ] = -value; + std::push_heap (heap.begin(),heap.end()); + } + } + + + std::sort_heap (heap.begin(),heap.end()); + for ( size_t i = 0 ; i != heap.size() ; ++i ) + { + if ( heap[i] == std::numeric_limits<int>::max() ) + { + heap[i] = 0; + } + else + { + heap[i] *= -1; + } + } + + if ( dbg ) + { + std::cerr << "This is the heap after all the operations :\n"; + for ( size_t i = 0 ; i != heap.size() ; ++i ) + { + std::cout << heap[i] << " "; + } + std::cout <<std::endl; + } + + this->sorted_vector_of_distances = heap; +} + + + + +template < typename F> +void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distances_via_vector_sorting( size_t where_to_cut ) +{ + std::vector< double > distances; + distances.reserve( (size_t)(0.5 * this->intervals.size() * ( this->intervals.size() - 1 ) + this->intervals.size()) ); + F f; + + //for every pair of points in the diagram, compute the minimum of their distance, and distance of those points from diagonal + for ( size_t i = 0 ; i < this->intervals.size() ; ++i ) + { + //add distance of i-th point in the diagram from the diagonal to the distances vector + distances.push_back( f( this->intervals[i] , std::make_pair( 0.5*(this->intervals[i].first+this->intervals[i].second) , 0.5*(this->intervals[i].first+this->intervals[i].second) ) ) ); + for ( size_t j = i+1 ; j < this->intervals.size() ; ++j ) + { + double value = std::min( + f( this->intervals[i] , this->intervals[j] ), + std::min( + f( this->intervals[i] , std::make_pair( 0.5*(this->intervals[i].first+this->intervals[i].second) , 0.5*(this->intervals[i].first+this->intervals[i].second) ) ), + f( this->intervals[j] , std::make_pair( 0.5*(this->intervals[j].first+this->intervals[j].second) , 0.5*(this->intervals[j].first+this->intervals[j].second) ) ) + ) + ); + distances.push_back( value ); + + } + } + std::sort( distances.begin() , distances.end() , std::greater<double>() ); + if ( distances.size() > where_to_cut )distances.resize( where_to_cut ); + + this->sorted_vector_of_distances = distances; +} + + + +//Implementations of functions for various concepts. +template <typename F> +double Vector_distances_in_diagram<F>::project_to_R( int number_of_function )const +{ + if ( (size_t)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"; + + double result = 0; + for ( size_t i = 0 ; i != (size_t)number_of_function ; ++i ) + { + result += sorted_vector_of_distances[i]; + } + return result; +} + +template <typename F> +void Vector_distances_in_diagram<F>::compute_average( const std::vector< Vector_distances_in_diagram* >& to_average ) +{ + + if ( to_average.size() == 0 ) + { + (*this) = Vector_distances_in_diagram<F>(); + return; + } + + size_t maximal_length_of_vector = 0; + for ( size_t i = 0 ; i != to_average.size() ; ++i ) + { + if ( to_average[i]->sorted_vector_of_distances.size() > maximal_length_of_vector ) + { + maximal_length_of_vector = to_average[i]->sorted_vector_of_distances.size(); + } + } + + std::vector< double > av( maximal_length_of_vector , 0 ); + for ( size_t i = 0 ; i != to_average.size() ; ++i ) + { + for ( size_t j = 0 ; j != to_average[i]->sorted_vector_of_distances.size() ; ++j ) + { + av[j] += to_average[i]->sorted_vector_of_distances[j]; + } + } + + for ( size_t i = 0 ; i != maximal_length_of_vector ; ++i ) + { + av[i] /= (double)to_average.size(); + } + this->sorted_vector_of_distances = av; + this->where_to_cut = av.size(); +} + +template <typename F> +double Vector_distances_in_diagram<F>::distance( const Vector_distances_in_diagram& second_ , double power )const +{ + bool dbg = false; + + if ( dbg ) + { + std::cerr << "Entering double Vector_distances_in_diagram<F>::distance( const Abs_Topological_data_with_distances* second , double power ) procedure \n"; + std::cerr << "Power : " << power << std::endl; + std::cerr << "This : " << *this << std::endl; + std::cerr << "second : " << second_ << std::endl; + } + + + double result = 0; + for ( size_t i = 0 ; i != std::min(this->sorted_vector_of_distances.size(), second_.sorted_vector_of_distances.size()) ; ++i ) + { + if ( power == 1 ) + { + if ( dbg ) + { + std::cerr << "|" << this->sorted_vector_of_distances[i] << " - " << second_.sorted_vector_of_distances[i] << " | : " << fabs( this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i] ) <<std::endl; + } + result += fabs( this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i] ); + } + else + { + if ( power < std::numeric_limits<double>::max() ) + { + result += std::pow( fabs( this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i] ) , power ); + } + else + { + //nax morm + if ( result < fabs( this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i] ) )result = fabs( this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i] ); + } + if ( dbg ) + { + std::cerr << "| " << this->sorted_vector_of_distances[i] << " - " << second_.sorted_vector_of_distances[i] << " : " << fabs( this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i] ) <<std::endl; + } + } + } + if ( this->sorted_vector_of_distances.size() != second_.sorted_vector_of_distances.size() ) + { + if ( this->sorted_vector_of_distances.size() > second_.sorted_vector_of_distances.size() ) + { + for ( size_t i = second_.sorted_vector_of_distances.size() ; i != this->sorted_vector_of_distances.size() ; ++i ) + { + result += fabs( this->sorted_vector_of_distances[i] ); + } + } + else + { + //this->sorted_vector_of_distances.size() < second_.sorted_vector_of_distances.size() + for ( size_t i = this->sorted_vector_of_distances.size() ; i != second_.sorted_vector_of_distances.size() ; ++i ) + { + result += fabs( second_.sorted_vector_of_distances[i] ); + } + } + } + + + if ( power != 1 ) + { + result = std::pow( result , (1.0/power) ); + } + return result; +} + +template < typename F> +std::vector<double> Vector_distances_in_diagram<F>::vectorize( int number_of_function )const +{ + if ( (size_t)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"; + + std::vector< double > result( std::min( (size_t)number_of_function , this->sorted_vector_of_distances.size() ) ); + for ( size_t i = 0 ; i != std::min( (size_t)number_of_function , this->sorted_vector_of_distances.size() ) ; ++i ) + { + result[i] = this->sorted_vector_of_distances[i]; + } + return result; +} + + +template < typename F> +void Vector_distances_in_diagram<F>::write_to_file( const char* filename )const +{ + std::ofstream out; + out.open( filename ); + + for ( size_t i = 0 ; i != this->sorted_vector_of_distances.size() ; ++i ) + { + out << this->sorted_vector_of_distances[i] << " "; + } + + out.close(); +} + +template < typename F> +void Vector_distances_in_diagram<F>::load_from_file( const char* filename ) +{ + std::ifstream in; + in.open( filename ); + //check if the file exist. + if ( !in.good() ) + { + 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"; + } + + double number; + while ( true ) + { + in >> number; + if ( in.eof() )break; + this->sorted_vector_of_distances.push_back(number); + } + in.close(); +} + +template < typename F> +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 ) + { + result += this->sorted_vector_of_distances[i] * second_vector.sorted_vector_of_distances[i]; + } + return result; +} + + + + + +}//namespace Gudhi_stat +}//namespace Gudhi + + +#endif // Vector_distances_in_diagram_H diff --git a/src/Gudhi_stat/include/gudhi/common_persistence_representations.h b/src/Gudhi_stat/include/gudhi/common_persistence_representations.h new file mode 100644 index 00000000..f223079a --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/common_persistence_representations.h @@ -0,0 +1,154 @@ +/* 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 common_gudhi_stat_H +#define common_gudhi_stat_H + + +namespace Gudhi +{ +namespace Persistence_representations +{ + //this file contain an implementation of some common procedures used in Gudhi_stat. + + //double epsi = std::numeric_limits<double>::epsilon(); +double epsi = 0.000005; + + + + + +/** + * A procedure used to compare doubles. Typically gien two doubles A and B, comparing A == B is not good idea. In this case, we use the procedure almostEqual with the epsi defined at + * the top of the file. Setting up the epsi give the user a tolerance on what should be consider equal. +**/ +inline bool almost_equal( double a , double b ) +{ + if ( fabs(a-b) < epsi ) + return true; + return false; +} + + + + +//landscapes +/** + * Extra functions needed in construction of barcodes. +**/ +double minus_length( std::pair<double,double> a ) +{ + return a.first-a.second; +} +double birth_plus_deaths( std::pair<double,double> a ) +{ + return a.first+a.second; +} + + +//landscapes +/** + * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two points. +**/ +std::pair<double,double> compute_parameters_of_a_line( std::pair<double,double> p1 , std::pair<double,double> p2 ) +{ + double a = (p2.second-p1.second)/( p2.first - p1.first ); + double b = p1.second - a*p1.first; + return std::make_pair(a,b); +} + +//landscapes +/** + * This procedure given two points which lies on the opposide sides of x axis, compute x for which the line connecting those two points crosses x axis. +**/ +double find_zero_of_a_line_segment_between_those_two_points ( std::pair<double,double> p1, std::pair<double,double> p2 ) +{ + if ( p1.first == p2.first )return p1.first; + if ( p1.second*p2.second > 0 ) + { + std::ostringstream errMessage; + errMessage <<"In function find_zero_of_a_line_segment_between_those_two_points the agguments are: (" << p1.first << "," << p1.second << ") and (" << p2.first << "," << p2.second << "). There is no zero in line between those two points. Program terminated."; + std::string errMessageStr = errMessage.str(); + const char* err = errMessageStr.c_str(); + throw(err); + } + //we assume here, that x \in [ p1.first, p2.first ] and p1 and p2 are points between which we will put the line segment + double a = (p2.second - p1.second)/(p2.first - p1.first); + double b = p1.second - a*p1.first; + //cerr << "Line crossing points : (" << p1.first << "," << p1.second << ") oraz (" << p2.first << "," << p2.second << ") : \n"; + //cerr << "a : " << a << " , b : " << b << " , x : " << x << endl; + return -b/a; +} + + + +//landscapes +/** + * This method provides a comparision of points that is used in construction of persistence landscapes. The orderign is lexicographical for the first coordinate, and reverse-lexicographical for the + * second coordinate. +**/ +bool compare_points_sorting( std::pair<double,double> f, std::pair<double,double> s ) +{ + if ( f.first < s.first ) + { + return true; + } + else + {//f.first >= s.first + if ( f.first > s.first ) + { + return false; + } + else + {//f.first == s.first + if ( f.second > s.second ) + { + return true; + } + else + { + return false; + } + } + } +} + + +//landscapes +/** + * This procedure takes two points in R^2 and a double value x. It conputes the line pasing through those two points and return the value of that linear function at x. +**/ +double function_value ( std::pair<double,double> p1, std::pair<double,double> p2 , double x ) +{ + //we assume here, that x \in [ p1.first, p2.first ] and p1 and p2 are points between which we will put the line segment + double a = (p2.second - p1.second)/(p2.first - p1.first); + double b = p1.second - a*p1.first; + return (a*x+b); +} + + + + +}//namespace Gudhi_stat +}//namespace Gudhi + +#endif diff --git a/src/Gudhi_stat/include/gudhi/read_persistence_from_file.h b/src/Gudhi_stat/include/gudhi/read_persistence_from_file.h new file mode 100644 index 00000000..058c77a4 --- /dev/null +++ b/src/Gudhi_stat/include/gudhi/read_persistence_from_file.h @@ -0,0 +1,300 @@ +/* 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 Read_Persitence_From_File_H +#define Read_Persitence_From_File_H + +#include <iostream> +#include <fstream> +#include <sstream> +#include <vector> +#include <algorithm> + + + +namespace Gudhi +{ +namespace Persistence_representations +{ + + + +/** + * This procedure reads names of files which are stored in a file. +**/ +std::vector< std::string > readFileNames( const char* filenameWithFilenames ) +{ + bool dbg = false; + + std::ifstream in(filenameWithFilenames); + if ( !in.good() ) + { + std::cerr << "The file : " << filenameWithFilenames << " do not exist. The program will now terminate \n"; + throw "The file from which you are trying to read do not exist. The program will now terminate \n"; + } + + std::vector< std::string > result; + std::string line; + while (!in.eof()) + { + getline(in,line); + line.erase( std::remove_if( line.begin(), line.end(), ::isspace) , line.end() ); + + if (dbg){std::cerr << "line : " << line << std::endl;} + + if ( (line.length() == 0) || (line[0] == '#') ) + { + //in this case we have a file name. First we should remove all the white spaces. + if ( dbg ){std::cerr << "This is a line with comment, it will be ignored n";} + } + else + { + result.push_back( line ); + if (dbg){std::cerr << "Line after removing white spaces : " << line << std::endl;} + } + } + in.close(); + + return result; +}//readFileNames + + + +std::vector< std::vector< double > > read_numbers_from_file_line_by_line( const char* filename ) +{ + bool dbg = false; + std::ifstream in(filename); + if ( !in.good() ) + { + 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::string line; + while ( in.good() ) + { + std::getline(in,line); + std::stringstream 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 + + +/** + * Universal procedure to read files with persistence. It ignores the lines starting from # (treat them as comments). + * It reads the fist line which is not a comment and assume that there are some numerical entries over there. The program assume + * that each other line in the file, which is not a comment, have the same number of numerical entries (2, 3 or 4). + * If there are two numerical entries per line, then the function assume that they are birth/death coordinates. + * If there are three numerical entries per line, then the function assume that they are: dimension and birth/death coordinates. + * If there are four numerical entries per line, then the function assume that they are: thc characteristic of a filed over which + * persistence was computed, dimension and birth/death coordinates. + * The 'inf' string can appear only as a last element of a line. + * The procedure returns vector of persistence pairs. +**/ +std::vector<std::pair<double,double> > read_persistence_intervals_in_one_dimension_from_file(std::string const& filename, int dimension=-1 , double what_to_substitute_for_infinite_bar = -1 ) +{ + bool dbg = false; + std::ifstream in; + in.open( filename ); + //checking if the file exist: + if ( !in.good() ) + { + 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::string line; + std::vector< std::pair<double,double> > barcode; + + int number_of_entries_per_line = -1; + + while (!in.eof()) + { + getline(in,line); + if ( dbg )std::cerr << "Reading line : " << line << std::endl; + if ( !(line.length() == 0 || line[0] == '#') ) + { + //If we do not know how many entries per line we have, we check it in below. + if ( number_of_entries_per_line == -1 ) + { + number_of_entries_per_line = 0; + std::string line_copy(line); + if ( line_copy.find("inf") != std::string::npos ) + { + size_t np = line_copy.find("inf"); + //replace symbols 'inf' in line_copy with whilespaces: + line_copy[np] = ' '; + line_copy[np+1] = ' '; + line_copy[np+2] = ' '; + number_of_entries_per_line = 1; + } + //check how many entries we have in the line. + std::stringstream ss( line_copy ); + double number; + std::vector<double> this_line; + while ( ss >> number ) + { + this_line.push_back( number ); + } + number_of_entries_per_line += (int)this_line.size(); + if ( dbg ) + { + std::cerr << "number_of_entries_per_line : " << number_of_entries_per_line << ". This number was obtained by analyzing this line : " << line << std::endl; + } + if ( (number_of_entries_per_line < 2) || ( number_of_entries_per_line > 4 ) ) + { + std::cerr << "The input file you have provided have wrong number of numerical entries per line. The program will now terminate. \n"; + throw "The input file you have provided have wrong number of numerical entries per line. The program will now terminate. \n"; + } + } + //In case there is an 'inf' string in this line, we are dealing with this situation in below. + if ( line.find("inf") != std::string::npos ) + { + if ( dbg ) + { + std::cerr << "This line: " << line << " contains infinite interval. \n"; + } + //first we substitute inf by whitespaces: + size_t np = line.find("inf"); + line[np] = ' '; + line[np+1] = ' '; + line[np+2] = ' '; + if ( what_to_substitute_for_infinite_bar != -1 ) + { + double beginn, field, dim; + std::stringstream lineSS(line); + if ( number_of_entries_per_line == 4 )lineSS >> field; + if ( number_of_entries_per_line >= 3 ) + { + lineSS >> dim; + } + else + { + dim = dimension; + } + lineSS >> beginn; + if ( dim == dimension ) + { + if ( beginn > what_to_substitute_for_infinite_bar ) + { + barcode.push_back( std::make_pair( what_to_substitute_for_infinite_bar , beginn ) ); + } + else + { + barcode.push_back( std::make_pair( beginn , what_to_substitute_for_infinite_bar ) ); + } + if (dbg) + { + std::cerr << "this is the line that is going to the output : " << beginn << " , " << what_to_substitute_for_infinite_bar << std::endl; + } + } + } + else + { + //this is a line with infinity. Since the variable what_to_substitute_for_infinite_bar have not been set up, it means that this line will be skipped. + if ( dbg ) + { + std::cerr << "We will skip it \n"; + } + } + continue; + } + else + { + //Then, we read the content of the line. We know that it do not contain 'inf' substring. + std::stringstream lineSS(line); + double beginn, endd, field, dim; + if ( number_of_entries_per_line == 4 )lineSS >> field; + if ( number_of_entries_per_line >= 3 ) + { + lineSS >> dim; + } + else + { + dim = dimension; + } + lineSS >> beginn; + lineSS >> endd; + if ( beginn > endd ) + { + std::swap(beginn,endd); + } + if ( dim == dimension ) + { + barcode.push_back( std::make_pair( beginn , endd ) ); + if (dbg) + { + std::cerr << "This is a line that is going to the output : " << beginn << " , " << endd << std::endl; + } + } + else + { + if ( (number_of_entries_per_line==3) && (dimension == -1) ) + { + barcode.push_back( std::make_pair( beginn , endd ) ); + } + } + } + } + else + { + if ( dbg ) + { + std::cerr << "This is a comment line \n"; + } + } + } + in.close(); + if ( dbg )std::cerr << "End of reading \n"; + + return barcode; +}//read_persistence_intervals_in_one_dimension_from_file + +}//namespace Gudhi_stat +}//namespace Gudhi + + + + +#endif + |