diff options
Diffstat (limited to 'src/Persistence_representations/include')
10 files changed, 5703 insertions, 0 deletions
diff --git a/src/Persistence_representations/include/gudhi/PSSK.h b/src/Persistence_representations/include/gudhi/PSSK.h new file mode 100644 index 00000000..fc90d0f4 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/PSSK.h @@ -0,0 +1,156 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PSSK_H_ +#define PSSK_H_ + +// gudhi include +#include <gudhi/Persistence_heat_maps.h> + +#include <limits> +#include <utility> +#include <vector> + +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 = + static_cast<int>((intervals_[pt_nr].first - this->min_) / (this->max_ - this->min_) * number_of_pixels); + int y_grid = + static_cast<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 corner. 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 + +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // PSSK_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h new file mode 100644 index 00000000..b1af3503 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h @@ -0,0 +1,993 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko and Mathieu Carriere + * + * Modifications: + * - 2018/04 MC: Add discrete/non-discrete mechanism and non-discrete version + * + * Copyright (C) 2019 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_HEAT_MAPS_H_ +#define PERSISTENCE_HEAT_MAPS_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> + +// standard include +#include <vector> +#include <sstream> +#include <iostream> +#include <cmath> +#include <limits> +#include <algorithm> +#include <utility> +#include <string> +#include <functional> + +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 initialize \n"; + std::cerr << "pixel_radius : " << pixel_radius << std::endl; + std::cerr << "kernel.size() : " << kernel.size() << std::endl; + getchar(); + } + + for (int x = -pixel_radius; x <= static_cast<int>(pixel_radius); x++) { + for (int y = -pixel_radius; y <= static_cast<int>(pixel_radius); y++) { + double real_x = 2 * sigma * x / pixel_radius; + double real_y = 2 * sigma * y / pixel_radius; + r = std::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 points 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 points depending on their persistence and/or location in the + *diagram. + * This particular functionality is a function 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 points 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 std::sqrt(std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + + std::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 points 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 std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + + std::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 points 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 points 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) : length_of_maximal_interval(len) {} + double operator()(const std::pair<double, double>& point_in_diagram) { + return (point_in_diagram.second - point_in_diagram.first) / this->length_of_maximal_interval; + } + + private: + double length_of_maximal_interval; +}; + +/** + * \class Persistence_heat_maps Persistence_heat_maps.h gudhi/Persistence_heat_maps.h + * \brief A class implementing persistence heat maps. + * + * \ingroup Persistence_representations + **/ + +// 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 Gaussian filter generated by create_Gaussian_filter filter (the default value of this variable 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 default. 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 name of a file with persistence intervals. The file should be readable by the function + *read_persistence_intervals_in_one_dimension_from_file. All other parameters are optional. They are: + * (2) a Gaussian filter generated by create_Gaussian_filter filter (the default value of this variable 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()); + + /** + * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for + * x and y axes) + **/ + Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x = 0, double max_x = 1, + double min_y = 0, double max_y = 1); + + /** + * Construction that takes only the diagram as input (weight and 2D kernel are template parameters) + **/ + Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel); + + /** + * 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 persistence 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 erase whatever 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 visualize 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 persistence 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 on the + *persistence heat map to a real line. + * At the moment this function is not tested, since it is quite likely to be changed in the future. Given this, when + *using it, keep in mind that it + * will be most likely 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 function is 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; + } + + // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false) + bool discrete = true; + std::function<double(std::pair<double, double>, std::pair<double, double>)> kernel; + std::vector<std::pair<double, double> > interval; + std::vector<double> weights; + + // data + Scalling_of_kernels f; + bool erase_below_diagonal; + double min_; + double max_; + std::vector<std::vector<double> > heat_map; +}; + +template <typename Scalling_of_kernels> +Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( + const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x, double max_x, double min_y, double max_y) { + this->discrete = true; + this->min_ = min_x; + this->max_ = max_x; + this->heat_map.resize(number_of_y_pixels); + double step_x = (max_x - min_x) / (number_of_x_pixels - 1); + double step_y = (max_y - min_y) / (number_of_y_pixels - 1); + + int num_pts = interval.size(); + for (size_t i = 0; i < number_of_y_pixels; i++) { + double y = min_y + i * step_y; + this->heat_map[i].reserve(number_of_x_pixels); + for (size_t j = 0; j < number_of_x_pixels; j++) { + double x = min_x + j * step_x; + std::pair<double, double> grid_point(x, y); + double pixel_value = 0; + for (int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point); + this->heat_map[i].push_back(pixel_value); + } + } + this->set_up_parameters_for_basic_classes(); +} + +template <typename Scalling_of_kernels> +Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( + const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel) { + this->discrete = false; + this->interval = interval; + this->kernel = kernel; + int num_pts = this->interval.size(); + for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i])); + this->set_up_parameters_for_basic_classes(); +} + +// 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 determined 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 = + static_cast<int>((intervals_[pt_nr].first - this->min_) / (this->max_ - this->min_) * number_of_pixels); + int y_grid = + static_cast<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 corner. 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); + } + 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 / static_cast<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 gnuplot_script; + gnuplot_script << filename << "_GnuplotScript"; + + out.open(gnuplot_script.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 << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" + << gnuplot_script.str().c_str() << "\'\"" << 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 persistence landscape file 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.good()) { + 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 { + std::vector<double> result; + if (!discrete) { + std::cout << "No vectorize method in case of infinite dimensional vectorization" << std::endl; + return result; + } + + // 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(); + } + + 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 { + if (this->discrete) { + // 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 non compatible sizes. We cannot therefore compute distance between " + "them. The program will now terminate"; + throw "The persistence images are of non compatible sizes. The program will now terminate"; + } + + // if we are here, we know that the two persistence images 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 += std::pow(std::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 < std::fabs(this->heat_map[i][j] - second.heat_map[i][j])) { + distance = std::fabs(this->heat_map[i][j] - second.heat_map[i][j]); + } + } + } + } + return distance; + } else { + return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) - + 2 * this->compute_scalar_product(second)); + } +} + +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 { + if (discrete) { + // 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 non compatible sizes. We cannot therefore compute distance between " + "them. The program will now terminate"; + throw "The persistence images are of non compatible sizes. The program will now terminate"; + } + + // if we are here, we know that the two persistence images 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; + } else { + int num_pts1 = this->interval.size(); + int num_pts2 = second.interval.size(); + double kernel_val = 0; + for (int i = 0; i < num_pts1; i++) { + std::pair<double, double> pi = this->interval[i]; + for (int j = 0; j < num_pts2; j++) { + std::pair<double, double> pj = second.interval[j]; + kernel_val += this->weights[i] * second.weights[j] * this->kernel(pi, pj); + } + } + return kernel_val; + } +} + +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_HEAT_MAPS_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_intervals.h b/src/Persistence_representations/include/gudhi/Persistence_intervals.h new file mode 100644 index 00000000..e2db4572 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_intervals.h @@ -0,0 +1,558 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#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 <cmath> +#include <functional> +#include <utility> +#include <string> + +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 double (or convertible 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 birth-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 elements 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 interval'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 lengths 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 sum 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 function 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 parameter 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::stringstream gnuplot_script; + gnuplot_script << filename << "_GnuplotScript"; + + out.open(gnuplot_script.str().c_str()); + + 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 << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" + << gnuplot_script.str().c_str() << "\'\"" << std::endl; + } + + /** +* Return number 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 projecting the persistence intervals to a real number. The function we use here is a sum + *of squared lengths 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 likely to be changed in the future. Given this, when + *using it, keep in mind that it + * will be most likely 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 family of vectors obtained from the persistence diagram. The i-th vector consist of the length 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 persistence 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. + + // 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) { + 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) / static_cast<double>(number_of_bins)) * + (this->intervals[i].second - this->intervals[i].first); + } + 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; + } + } + 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])); + } + // also add a distance from this guy to diagonal: + 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 Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_INTERVALS_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_intervals_with_distances.h b/src/Persistence_representations/include/gudhi/Persistence_intervals_with_distances.h new file mode 100644 index 00000000..98543f2f --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_intervals_with_distances.h @@ -0,0 +1,51 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_INTERVALS_WITH_DISTANCES_H_ +#define PERSISTENCE_INTERVALS_WITH_DISTANCES_H_ + +#include <gudhi/Persistence_intervals.h> +#include <gudhi/Bottleneck.h> + +#include <limits> + +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 = (std::numeric_limits<double>::min)()) 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 Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_INTERVALS_WITH_DISTANCES_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_landscape.h b/src/Persistence_representations/include/gudhi/Persistence_landscape.h new file mode 100644 index 00000000..b819ccb6 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_landscape.h @@ -0,0 +1,1371 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_LANDSCAPE_H_ +#define PERSISTENCE_LANDSCAPE_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> + +// standard include +#include <cmath> +#include <iostream> +#include <vector> +#include <limits> +#include <fstream> +#include <sstream> +#include <algorithm> +#include <string> +#include <utility> +#include <functional> + +namespace Gudhi { +namespace Persistence_representations { + +// pre declaration +class Persistence_landscape; +template <typename operation> +Persistence_landscape operation_on_pair_of_landscapes(const Persistence_landscape& land1, + const Persistence_landscape& land2); + +/** + * \class Persistence_landscape Persistence_landscape.h gudhi/Persistence_landscape.h + * \brief A class implementing persistence landscapes data structures. + * + * \ingroup Persistence_representations + * + * \details + * For theoretical description, please consult <i>Statistical topological data analysis using persistence + * landscapes</i>\cite bubenik_landscapes_2015 , and for details of algorithms, + * <i>A persistence landscapes toolbox for topological statistics</i>\cite bubenik_dlotko_landscapes_2016. + * + * Persistence landscapes allow vectorization, 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 rounding 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, + size_t number_of_levels = std::numeric_limits<size_t>::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. + **/ + Persistence_landscape(const char* filename, size_t dimension = std::numeric_limits<unsigned>::max(), + size_t number_of_levels = std::numeric_limits<size_t>::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 multiplication 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 multiplication 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); + + /** + *\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 between 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(); + + Persistence_landscape* new_abs(); + + /** + * Computes the number of landscape functions. + **/ + size_t size() const { return this->land.size(); } + + /** + * Compute 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 likely to be changed in the future. Given this, when + *using it, keep in mind that it + * will be most likely 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(PD) 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 landscape. 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 number of intermediate complexes. We + // have to clean that up, but we cannot erase the initial landscapes 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 using 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 / static_cast<double>(to_average.size()); + } + + /** + * A function to compute distance between persistence landscape. + * The parameter of this function is 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 function is 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 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, + size_t number_of_levels = std::numeric_limits<size_t>::max()); + 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, size_t number_of_levels) { + 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, number_of_levels); + 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))) { + 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, + size_t number_of_levels) { + this->construct_persistence_landscape_from_barcode(p, number_of_levels); + 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, size_t number_of_levels) { + 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; + size_t number_of_levels_in_the_landscape = 0; + 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); + + ++number_of_levels_in_the_landscape; + if (number_of_levels == number_of_levels_in_the_landscape) { + break; + } + } +} + +// 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 up 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 integral. + 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 up 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; + } + } + } + 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 o result + 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::new_abs() { + Persistence_landscape* result = new Persistence_landscape(*this); + 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 o result + 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 before head + 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 persistence landscape file do not exist. The program will now terminate \n"; + } + + std::string line; + std::vector<std::pair<double, double> > landscapeAtThisLevel; + + bool isThisAFirsLine = true; + while (in.good()) { + 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 point : " << 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"; + } + + 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"; + } + lambda_n.push_back( + std::make_pair(land1.land[i][p].first, + oper(static_cast<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; + // In this case, I consider points at the infinity + for (size_t i = 1; i != pl1.land[level].size() - 1; ++i) { + 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(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.0 / 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(); + } + auto&& l1_land_level = l1.land[level]; + auto&& l2_land_level = l2.land[level]; + + 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"; + } + } + + if ( l1It + 1 >= l1_land_level.size() )break; + if ( l2It + 1 >= l2_land_level.size() )break; + + // 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 gnuplot_script; + gnuplot_script << filename << "_GnuplotScript"; + out.open(gnuplot_script.str().c_str()); + + 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 << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" + << gnuplot_script.str().c_str() << "\'\"" << std::endl; +} + +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_LANDSCAPE_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h new file mode 100644 index 00000000..68bce336 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_landscape_on_grid.h @@ -0,0 +1,1335 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_LANDSCAPE_ON_GRID_H_ +#define PERSISTENCE_LANDSCAPE_ON_GRID_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> + +// standard include +#include <iostream> +#include <vector> +#include <limits> +#include <fstream> +#include <sstream> +#include <algorithm> +#include <cmath> +#include <functional> +#include <utility> +#include <string> +#include <cstdint> + +namespace Gudhi { +namespace Persistence_representations { + +// pre declaration +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); + +/** + * \class Persistence_landscape_on_grid Persistence_landscape_on_grid.h gudhi/Persistence_landscape_on_grid.h + * \brief A class implementing persistence landscapes by approximating them on a collection of grid points. + * + * \ingroup Persistence_representations + * + * \details + * Persistence landscapes on grid allows vectorization, 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 rounding 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, resolution 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, + uint16_t dimension_ = std::numeric_limits<uint16_t>::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, resolution 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_, + uint16_t dimension_ = std::numeric_limits<uint16_t>::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 resolution of a grid and the number of landscape + * functions to be created. The remaining parameters are calculated based on data. + **/ + Persistence_landscape_on_grid(const char* filename, size_t number_of_points, unsigned number_of_levels_of_landscape, + uint16_t dimension = std::numeric_limits<uint16_t>::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 resolution 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 remaining parameters are calculated based on data. + **/ + Persistence_landscape_on_grid(const char* filename, size_t number_of_points, + uint16_t dimension = std::numeric_limits<uint16_t>::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) / static_cast<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) / static_cast<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 << "result : " << 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) / static_cast<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) / static_cast<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 << "position : " << position << std::endl; + } + // check if we are not exactly 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 { + 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 multiplication 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 multiplication 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 persistence landscape. + **/ + Persistence_landscape_on_grid operator+=(const Persistence_landscape_on_grid& rhs) { + *this = *this + rhs; + return *this; + } + + /** + * Operator -=. The second parameter is persistence landscape. + **/ + 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 incompatible 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); + } + + /** + * 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(); } + + /** + * 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); + + /** + * Function to compute absolute value of a PL function. The representation of persistence landscapes allow to store + *general PL-function. When computing distance between 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(); } + + /** + * Compute 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) / static_cast<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 integral 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 i.e. : " << 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 intersecting. 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(p); + if (dbg) { + std::cerr << "integral : " << result << std::endl; + getchar(); + } + } else { + result = lan.compute_integral_of_landscape(); + if (dbg) { + std::cerr << "integral, without power : " << result << std::endl; + getchar(); + } + } + // (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p) + return pow(result, 1.0 / 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 likely to be changed in the future. Given this, when + *using it, keep in mind that it + * will be most likely 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(PD) 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 landscape. 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 average, 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 arithmetic 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] /= static_cast<double>(to_average.size()); + } + } + } // compute_average + + /** + * A function to compute distance between persistence landscape on a grid. + * The parameter of this function is 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 function is 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 landscapes. It can be used for visualization + **/ + std::vector<std::vector<double> > output_for_visualization() const { return this->values_of_landscapes; } + + /** + * function used to create a gnuplot script for visualization of landscapes. 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 default 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 program will now terminate.\n"; + } + + double dx = (grid_max_ - grid_min_) / static_cast<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 heap. + 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 heap. + 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 << "Adding 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 minus 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) { + 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_, uint16_t dimension) { + std::vector<std::pair<double, double> > p; + if (dimension == std::numeric_limits<uint16_t>::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, + uint16_t dimension) { + std::vector<std::pair<double, double> > p; + if (dimension == std::numeric_limits<uint16_t>::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_, + uint16_t dimension) { + std::vector<std::pair<double, double> > p; + if (dimension == std::numeric_limits<uint16_t>::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, + uint16_t dimension) { + std::vector<std::pair<double, double> > p; + if (dimension == std::numeric_limits<uint16_t>::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 persistence landscape file 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::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 gnuplot_script; + gnuplot_script << filename << "_GnuplotScript"; + out.open(gnuplot_script.str().c_str()); + + 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) / static_cast<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 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 << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" + << gnuplot_script.str().c_str() << "\'\"" << 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 perform 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 Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_LANDSCAPE_ON_GRID_H_ diff --git a/src/Persistence_representations/include/gudhi/Persistence_vectors.h b/src/Persistence_representations/include/gudhi/Persistence_vectors.h new file mode 100644 index 00000000..6776f4a3 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_vectors.h @@ -0,0 +1,628 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_VECTORS_H_ +#define PERSISTENCE_VECTORS_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/distance_functions.h> + +#include <fstream> +#include <cmath> +#include <algorithm> +#include <iostream> +#include <limits> +#include <functional> +#include <utility> +#include <vector> + +namespace Gudhi { +namespace Persistence_representations { + +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)); + } +}; + +/** + * \class Vector_distances_in_diagram Persistence_vectors.h gudhi/Persistence_vectors.h + * \brief A class implementing persistence vectors. + * + * \ingroup Persistence_representations + * + * \details + * This is an implementation of idea presented in the paper <i>Stable Topological Signatures for Points on 3D + * Shapes</i> \cite Carriere_Oudot_Ovsjanikov_top_signatures_3d .<br> + * The parameter of the class is the class that computes distance used to construct the vectors. The typical function + * is either Euclidean 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 desired 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 desired 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); + + /** + * Comparison 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 likely to be changed in the future. Given this, when + *using it, keep in mind that it + * will be most likely 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 persistence 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 visualize 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 visualize, install gnuplot and type the command: gnuplot -persist -e \"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); + } + + // arithmetic 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()) { + // 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] /= static_cast<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 { + // max norm + 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 persistence landscape file do not exist. The program will now terminate \n"; + } + + double number; + while (in >> number) { + 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 Persistence_representations +} // namespace Gudhi + +#endif // PERSISTENCE_VECTORS_H_ diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h new file mode 100644 index 00000000..e3ed2f6a --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -0,0 +1,380 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SLICED_WASSERSTEIN_H_ +#define SLICED_WASSERSTEIN_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Debug_utils.h> + +#include <vector> // for std::vector<> +#include <utility> // for std::pair<>, std::move +#include <algorithm> // for std::sort, std::max, std::merge +#include <cmath> // for std::abs, std::sqrt +#include <stdexcept> // for std::invalid_argument +#include <random> // for std::random_device + +namespace Gudhi { +namespace Persistence_representations { + +/** + * \class Sliced_Wasserstein gudhi/Sliced_Wasserstein.h + * \brief A class implementing the Sliced Wasserstein kernel. + * + * \ingroup Persistence_representations + * + * \details + * In this class, we compute infinite-dimensional representations of persistence diagrams by using the + * Sliced Wasserstein kernel (see \ref sec_persistence_kernels for more details on kernels). We recall that + * infinite-dimensional representations are defined implicitly, so only scalar products and distances are available for + * the representations defined in this class. + * The Sliced Wasserstein kernel is defined as a Gaussian-like kernel between persistence diagrams, where the distance + * used for comparison is the Sliced Wasserstein distance \f$SW\f$ between persistence diagrams, defined as the + * integral of the 1-norm between the sorted projections of the diagrams onto all lines passing through the origin: + * + * \f$ SW(D_1,D_2)=\int_{\theta\in\mathbb{S}}\,\|\pi_\theta(D_1\cup\pi_\Delta(D_2))-\pi_\theta(D_2\cup\pi_\Delta(D_1))\ + * |_1{\rm d}\theta\f$, + * + * where \f$\pi_\theta\f$ is the projection onto the line defined with angle \f$\theta\f$ in the unit circle + * \f$\mathbb{S}\f$, and \f$\pi_\Delta\f$ is the projection onto the diagonal. + * Assuming that the diagrams are in general position (i.e. there is no collinear triple), the integral can be computed + * exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points in the diagrams. We provide two + * approximations of the integral: one in which we slightly perturb the diagram points so that they are in general + * position, and another in which we approximate the integral by sampling \f$N\f$ lines in the circle in + * \f$O(Nn{\rm log}(n))\f$ time. The Sliced Wasserstein Kernel is then computed as: + * + * \f$ k(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$ + * + * The first method is usually much more accurate but also + * much slower. For more details, please see \cite pmlr-v70-carriere17a . + * + **/ + +class Sliced_Wasserstein { + protected: + Persistence_diagram diagram; + int approx; + double sigma; + std::vector<std::vector<double> > projections, projections_diagonal; + + // ********************************** + // Utils. + // ********************************** + + void build_rep() { + if (approx > 0) { + double step = pi / this->approx; + int n = diagram.size(); + + for (int i = 0; i < this->approx; i++) { + std::vector<double> l, l_diag; + for (int j = 0; j < n; j++) { + double px = diagram[j].first; + double py = diagram[j].second; + double proj_diag = (px + py) / 2; + + l.push_back(px * cos(-pi / 2 + i * step) + py * sin(-pi / 2 + i * step)); + l_diag.push_back(proj_diag * cos(-pi / 2 + i * step) + proj_diag * sin(-pi / 2 + i * step)); + } + + std::sort(l.begin(), l.end()); + std::sort(l_diag.begin(), l_diag.end()); + projections.push_back(std::move(l)); + projections_diagonal.push_back(std::move(l_diag)); + } + + diagram.clear(); + } + } + + // Compute the angle formed by two points of a PD + double compute_angle(const Persistence_diagram& diag, int i, int j) const { + if (diag[i].second == diag[j].second) + return pi / 2; + else + return atan((diag[j].first - diag[i].first) / (diag[i].second - diag[j].second)); + } + + // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in + // [0,pi] + double compute_int_cos(double alpha, double beta) const { + double res = 0; + if (alpha >= 0 && alpha <= pi) { + if (cos(alpha) >= 0) { + if (pi / 2 <= beta) { + res = 2 - sin(alpha) - sin(beta); + } else { + res = sin(beta) - sin(alpha); + } + } else { + if (1.5 * pi <= beta) { + res = 2 + sin(alpha) + sin(beta); + } else { + res = sin(alpha) - sin(beta); + } + } + } + if (alpha >= -pi && alpha <= 0) { + if (cos(alpha) <= 0) { + if (-pi / 2 <= beta) { + res = 2 + sin(alpha) + sin(beta); + } else { + res = sin(alpha) - sin(beta); + } + } else { + if (pi / 2 <= beta) { + res = 2 - sin(alpha) - sin(beta); + } else { + res = sin(beta) - sin(alpha); + } + } + } + return res; + } + + double compute_int(double theta1, double theta2, int p, int q, const Persistence_diagram& diag1, + const Persistence_diagram& diag2) const { + double norm = std::sqrt((diag1[p].first - diag2[q].first) * (diag1[p].first - diag2[q].first) + + (diag1[p].second - diag2[q].second) * (diag1[p].second - diag2[q].second)); + double angle1; + if (diag1[p].first == diag2[q].first) + angle1 = theta1 - pi / 2; + else + angle1 = theta1 - atan((diag1[p].second - diag2[q].second) / (diag1[p].first - diag2[q].first)); + double angle2 = angle1 + theta2 - theta1; + double integral = compute_int_cos(angle1, angle2); + return norm * integral; + } + + // Evaluation of the Sliced Wasserstein Distance between a pair of diagrams. + double compute_sliced_wasserstein_distance(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->approx == second.approx, + std::invalid_argument("Error: different approx values for representations")); + + Persistence_diagram diagram1 = this->diagram; + Persistence_diagram diagram2 = second.diagram; + double sw = 0; + + if (this->approx == -1) { + // Add projections onto diagonal. + int n1, n2; + n1 = diagram1.size(); + n2 = diagram2.size(); + double min_ordinate = std::numeric_limits<double>::max(); + double min_abscissa = std::numeric_limits<double>::max(); + double max_ordinate = std::numeric_limits<double>::lowest(); + double max_abscissa = std::numeric_limits<double>::lowest(); + for (int i = 0; i < n2; i++) { + min_ordinate = std::min(min_ordinate, diagram2[i].second); + min_abscissa = std::min(min_abscissa, diagram2[i].first); + max_ordinate = std::max(max_ordinate, diagram2[i].second); + max_abscissa = std::max(max_abscissa, diagram2[i].first); + diagram1.emplace_back((diagram2[i].first + diagram2[i].second) / 2, + (diagram2[i].first + diagram2[i].second) / 2); + } + for (int i = 0; i < n1; i++) { + min_ordinate = std::min(min_ordinate, diagram1[i].second); + min_abscissa = std::min(min_abscissa, diagram1[i].first); + max_ordinate = std::max(max_ordinate, diagram1[i].second); + max_abscissa = std::max(max_abscissa, diagram1[i].first); + diagram2.emplace_back((diagram1[i].first + diagram1[i].second) / 2, + (diagram1[i].first + diagram1[i].second) / 2); + } + int num_pts_dgm = diagram1.size(); + + // Slightly perturb the points so that the PDs are in generic positions. + double epsilon = 0.0001; + double thresh_y = (max_ordinate - min_ordinate) * epsilon; + double thresh_x = (max_abscissa - min_abscissa) * epsilon; + std::random_device rd; + std::default_random_engine re(rd()); + std::uniform_real_distribution<double> uni(-1, 1); + for (int i = 0; i < num_pts_dgm; i++) { + double u = uni(re); + diagram1[i].first += u * thresh_x; + diagram1[i].second += u * thresh_y; + diagram2[i].first += u * thresh_x; + diagram2[i].second += u * thresh_y; + } + + // Compute all angles in both PDs. + std::vector<std::pair<double, std::pair<int, int> > > angles1, angles2; + for (int i = 0; i < num_pts_dgm; i++) { + for (int j = i + 1; j < num_pts_dgm; j++) { + double theta1 = compute_angle(diagram1, i, j); + double theta2 = compute_angle(diagram2, i, j); + angles1.emplace_back(theta1, std::pair<int, int>(i, j)); + angles2.emplace_back(theta2, std::pair<int, int>(i, j)); + } + } + + // Sort angles. + std::sort(angles1.begin(), angles1.end(), + [](const std::pair<double, std::pair<int, int> >& p1, + const std::pair<double, std::pair<int, int> >& p2) { return (p1.first < p2.first); }); + std::sort(angles2.begin(), angles2.end(), + [](const std::pair<double, std::pair<int, int> >& p1, + const std::pair<double, std::pair<int, int> >& p2) { return (p1.first < p2.first); }); + + // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2). + std::vector<int> orderp1, orderp2; + for (int i = 0; i < num_pts_dgm; i++) { + orderp1.push_back(i); + orderp2.push_back(i); + } + std::sort(orderp1.begin(), orderp1.end(), [&](int i, int j) { + if (diagram1[i].second != diagram1[j].second) + return (diagram1[i].second < diagram1[j].second); + else + return (diagram1[i].first > diagram1[j].first); + }); + std::sort(orderp2.begin(), orderp2.end(), [&](int i, int j) { + if (diagram2[i].second != diagram2[j].second) + return (diagram2[i].second < diagram2[j].second); + else + return (diagram2[i].first > diagram2[j].first); + }); + + // Find the inverses of the orders. + std::vector<int> order1(num_pts_dgm); + std::vector<int> order2(num_pts_dgm); + for (int i = 0; i < num_pts_dgm; i++) { + order1[orderp1[i]] = i; + order2[orderp2[i]] = i; + } + + // Record all inversions of points in the orders as theta varies along the positive half-disk. + std::vector<std::vector<std::pair<int, double> > > anglePerm1(num_pts_dgm); + std::vector<std::vector<std::pair<int, double> > > anglePerm2(num_pts_dgm); + + int m1 = angles1.size(); + for (int i = 0; i < m1; i++) { + double theta = angles1[i].first; + int p = angles1[i].second.first; + int q = angles1[i].second.second; + anglePerm1[order1[p]].emplace_back(p, theta); + anglePerm1[order1[q]].emplace_back(q, theta); + int a = order1[p]; + int b = order1[q]; + order1[p] = b; + order1[q] = a; + } + + int m2 = angles2.size(); + for (int i = 0; i < m2; i++) { + double theta = angles2[i].first; + int p = angles2[i].second.first; + int q = angles2[i].second.second; + anglePerm2[order2[p]].emplace_back(p, theta); + anglePerm2[order2[q]].emplace_back(q, theta); + int a = order2[p]; + int b = order2[q]; + order2[p] = b; + order2[q] = a; + } + + for (int i = 0; i < num_pts_dgm; i++) { + anglePerm1[order1[i]].emplace_back(i, pi / 2); + anglePerm2[order2[i]].emplace_back(i, pi / 2); + } + + // Compute the SW distance with the list of inversions. + for (int i = 0; i < num_pts_dgm; i++) { + std::vector<std::pair<int, double> > u, v; + u = anglePerm1[i]; + v = anglePerm2[i]; + double theta1, theta2; + theta1 = -pi / 2; + unsigned int ku, kv; + ku = 0; + kv = 0; + theta2 = std::min(u[ku].second, v[kv].second); + while (theta1 != pi / 2) { + if (diagram1[u[ku].first].first != diagram2[v[kv].first].first || + diagram1[u[ku].first].second != diagram2[v[kv].first].second) + if (theta1 != theta2) sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2); + theta1 = theta2; + if ((theta2 == u[ku].second) && ku < u.size() - 1) ku++; + if ((theta2 == v[kv].second) && kv < v.size() - 1) kv++; + theta2 = std::min(u[ku].second, v[kv].second); + } + } + } else { + double step = pi / this->approx; + std::vector<double> v1, v2; + for (int i = 0; i < this->approx; i++) { + v1.clear(); + v2.clear(); + std::merge(this->projections[i].begin(), this->projections[i].end(), second.projections_diagonal[i].begin(), + second.projections_diagonal[i].end(), std::back_inserter(v1)); + std::merge(second.projections[i].begin(), second.projections[i].end(), this->projections_diagonal[i].begin(), + this->projections_diagonal[i].end(), std::back_inserter(v2)); + + int n = v1.size(); + double f = 0; + for (int j = 0; j < n; j++) f += std::abs(v1[j] - v2[j]); + sw += f * step; + } + } + + return sw / pi; + } + + public: + /** \brief Sliced Wasserstein kernel constructor. + * \implements Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product + * \ingroup Sliced_Wasserstein + * + * @param[in] _diagram persistence diagram. + * @param[in] _sigma bandwidth parameter. + * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set + * to -1 for random perturbation. If positive, then projections of the diagram points on all + * directions are stored in memory to reduce computation time. + * + */ + Sliced_Wasserstein(const Persistence_diagram& _diagram, double _sigma = 1.0, int _approx = 10) + : diagram(_diagram), approx(_approx), sigma(_sigma) { + build_rep(); + } + + /** \brief Evaluation of the kernel on a pair of diagrams. + * \ingroup Sliced_Wasserstein + * + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. + * + */ + double compute_scalar_product(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->sigma == second.sigma, + std::invalid_argument("Error: different sigma values for representations")); + return std::exp(-compute_sliced_wasserstein_distance(second) / (2 * this->sigma * this->sigma)); + } + + /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. + * \ingroup Sliced_Wasserstein + * + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. + * + */ + double distance(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->sigma == second.sigma, + std::invalid_argument("Error: different sigma values for representations")); + return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) - + 2 * this->compute_scalar_product(second)); + } + +}; // class Sliced_Wasserstein +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // SLICED_WASSERSTEIN_H_ diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h new file mode 100644 index 00000000..5eff0192 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h @@ -0,0 +1,123 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef COMMON_PERSISTENCE_REPRESENTATIONS_H_ +#define COMMON_PERSISTENCE_REPRESENTATIONS_H_ + +#include <utility> +#include <string> +#include <cmath> +#include <boost/math/constants/constants.hpp> + + + +namespace Gudhi { +namespace Persistence_representations { +// this file contain an implementation of some common procedures used in Persistence_representations. + +static constexpr double pi = boost::math::constants::pi<double>(); + + +/** + * In this module, we use the name Persistence_diagram for the representation of a diagram in a vector of pairs of two double. + */ +using Persistence_diagram = std::vector<std::pair<double, double> >; + +// double epsi = std::numeric_limits<double>::epsilon(); +double epsi = 0.000005; + +/** + * A procedure used to compare doubles. Typically given 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 gives the user a tolerance on what should be consider equal. +**/ +inline bool almost_equal(double a, double b) { + if (std::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 opposite 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 arguments 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; + return -b / a; +} + +// landscapes +/** + * This method provides a comparison of points that is used in construction of persistence landscapes. The ordering 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 computes the line parsing 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 Persistence_representations +} // namespace Gudhi + +#endif // COMMON_PERSISTENCE_REPRESENTATIONS_H_ diff --git a/src/Persistence_representations/include/gudhi/read_persistence_from_file.h b/src/Persistence_representations/include/gudhi/read_persistence_from_file.h new file mode 100644 index 00000000..5c2d2038 --- /dev/null +++ b/src/Persistence_representations/include/gudhi/read_persistence_from_file.h @@ -0,0 +1,108 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Pawel Dlotko + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef READ_PERSISTENCE_FROM_FILE_H_ +#define READ_PERSISTENCE_FROM_FILE_H_ + +#include <gudhi/reader_utils.h> + +#include <iostream> +#include <fstream> +#include <sstream> +#include <vector> +#include <algorithm> +#include <string> +#include <utility> +#include <limits> // for std::numeric_limits<> + +namespace Gudhi { +namespace Persistence_representations { + +/** + * 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: the 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::string line; + std::vector<std::pair<double, double> > barcode_initial = + read_persistence_intervals_in_dimension(filename, static_cast<int>(dimension)); + std::vector<std::pair<double, double> > final_barcode; + final_barcode.reserve(barcode_initial.size()); + + if (dbg) { + std::cerr << "Here are the intervals that we read from the file : \n"; + for (size_t i = 0; i != barcode_initial.size(); ++i) { + std::cout << barcode_initial[i].first << " " << barcode_initial[i].second << std::endl; + } + getchar(); + } + + for (size_t i = 0; i != barcode_initial.size(); ++i) { + if (dbg) { + std::cout << "COnsidering interval : " << barcode_initial[i].first << " " << barcode_initial[i].second + << std::endl; + } + + if (barcode_initial[i].first > barcode_initial[i].second) { + // note that in this case barcode_initial[i].second != std::numeric_limits<double>::infinity() + if (dbg) std::cout << "Swap and enter \n"; + // swap them to make sure that birth < death + final_barcode.push_back(std::pair<double, double>(barcode_initial[i].second, barcode_initial[i].first)); + continue; + } else { + if (barcode_initial[i].second != std::numeric_limits<double>::infinity()) { + if (dbg) std::cout << "Simply enters\n"; + // in this case, due to the previous conditions we know that barcode_initial[i].first < + // barcode_initial[i].second, so we put them as they are + final_barcode.push_back(std::pair<double, double>(barcode_initial[i].first, barcode_initial[i].second)); + } + } + + if ((barcode_initial[i].second == std::numeric_limits<double>::infinity()) && + (what_to_substitute_for_infinite_bar != -1)) { + if (barcode_initial[i].first < what_to_substitute_for_infinite_bar) { + // if only birth < death. + final_barcode.push_back( + std::pair<double, double>(barcode_initial[i].first, what_to_substitute_for_infinite_bar)); + } + } else { + // if the variable what_to_substitute_for_infinite_bar is not set, then we ignore all the infinite bars. + } + } + + if (dbg) { + std::cerr << "Here are the final bars that we are sending further : \n"; + for (size_t i = 0; i != final_barcode.size(); ++i) { + std::cout << final_barcode[i].first << " " << final_barcode[i].second << std::endl; + } + std::cerr << "final_barcode.size() : " << final_barcode.size() << std::endl; + getchar(); + } + + return final_barcode; +} // read_persistence_intervals_in_one_dimension_from_file + +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // READ_PERSISTENCE_FROM_FILE_H_ |