23 #ifndef PERSISTENCE_VECTORS_H_ 24 #define PERSISTENCE_VECTORS_H_ 27 #include <gudhi/read_persistence_from_file.h> 28 #include <gudhi/common_persistence_representations.h> 41 namespace Persistence_representations {
44 struct Maximum_distance {
45 double operator()(
const std::pair<T, T>& f,
const std::pair<T, T>& s) {
46 return std::max(fabs(f.first - s.first), fabs(f.second - s.second));
84 unsigned dimension = std::numeric_limits<unsigned>::max());
90 friend std::ostream& operator<<(std::ostream& out, const Vector_distances_in_diagram<K>& d) {
91 for (
size_t i = 0; i != std::min(d.sorted_vector_of_distances.size(), d.where_to_cut); ++i) {
92 out << d.sorted_vector_of_distances[i] <<
" ";
101 if (position >= this->sorted_vector_of_distances.size())
102 throw(
"Wrong position in accessing Vector_distances_in_diagram::sorted_vector_of_distances\n");
103 return this->sorted_vector_of_distances[position];
109 inline size_t size()
const {
return this->sorted_vector_of_distances.size(); }
114 void write_to_file(
const char* filename)
const;
119 void print_to_file(
const char* filename)
const { this->write_to_file(filename); }
124 void load_from_file(
const char* filename);
130 if (this->sorted_vector_of_distances.size() != second.sorted_vector_of_distances.size())
return false;
131 for (
size_t i = 0; i != this->sorted_vector_of_distances.size(); ++i) {
132 if (!almost_equal(this->sorted_vector_of_distances[i], second.sorted_vector_of_distances[i]))
return false;
147 double project_to_R(
int number_of_function)
const;
157 std::vector<double> vectorize(
int number_of_function)
const;
167 void compute_average(
const std::vector<Vector_distances_in_diagram*>& to_average);
190 void plot(
const char* filename)
const {
191 std::stringstream gnuplot_script;
192 gnuplot_script << filename <<
"_GnuplotScript";
194 out.open(gnuplot_script.str().c_str());
195 out <<
"set style data histogram" << std::endl;
196 out <<
"set style histogram cluster gap 1" << std::endl;
197 out <<
"set style fill solid border -1" << std::endl;
198 out <<
"plot '-' notitle" << std::endl;
199 for (
size_t i = 0; i != this->sorted_vector_of_distances.size(); ++i) {
200 out << this->sorted_vector_of_distances[i] << std::endl;
204 std::cout <<
"To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" 205 << gnuplot_script.str().c_str() <<
"\'\"" << std::endl;
211 std::pair<double, double>
get_x_range()
const {
return std::make_pair(0, this->sorted_vector_of_distances.size()); }
217 if (this->sorted_vector_of_distances.size() == 0)
return std::make_pair(0, 0);
218 return std::make_pair(this->sorted_vector_of_distances[0], 0);
222 template <
typename Operation_type>
225 Operation_type opertion) {
228 result.sorted_vector_of_distances.reserve(
229 std::max(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size()));
230 for (
size_t i = 0; i != std::min(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size());
232 result.sorted_vector_of_distances.push_back(
233 opertion(first.sorted_vector_of_distances[i], second.sorted_vector_of_distances[i]));
235 if (first.sorted_vector_of_distances.size() ==
236 std::min(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size())) {
237 for (
size_t i = std::min(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size());
238 i != std::max(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size()); ++i) {
239 result.sorted_vector_of_distances.push_back(opertion(0, second.sorted_vector_of_distances[i]));
242 for (
size_t i = std::min(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size());
243 i != std::max(first.sorted_vector_of_distances.size(), second.sorted_vector_of_distances.size()); ++i) {
244 result.sorted_vector_of_distances.push_back(opertion(first.sorted_vector_of_distances[i], 0));
255 result.sorted_vector_of_distances.reserve(this->sorted_vector_of_distances.size());
256 for (
size_t i = 0; i != this->sorted_vector_of_distances.size(); ++i) {
257 result.sorted_vector_of_distances.push_back(scalar * this->sorted_vector_of_distances[i]);
267 return operation_on_pair_of_vectors(first, second, std::plus<double>());
274 return operation_on_pair_of_vectors(first, second, std::minus<double>());
317 if (x == 0)
throw(
"In operator /=, division by 0. Program terminated.");
318 *
this = *
this * (1 / x);
323 std::vector<std::pair<double, double> > intervals;
324 std::vector<double> sorted_vector_of_distances;
325 size_t number_of_functions_for_vectorization;
326 size_t number_of_functions_for_projections_to_reals;
329 void compute_sorted_vector_of_distances_via_heap(
size_t where_to_cut);
330 void compute_sorted_vector_of_distances_via_vector_sorting(
size_t where_to_cut);
333 : sorted_vector_of_distances(sorted_vector_of_distances_) {
334 this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
337 void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals() {
339 this->number_of_functions_for_vectorization = this->sorted_vector_of_distances.size();
340 this->number_of_functions_for_projections_to_reals = this->sorted_vector_of_distances.size();
344 template <
typename F>
346 size_t where_to_cut_)
347 : where_to_cut(where_to_cut_) {
348 std::vector<std::pair<double, double> > i(intervals_);
351 this->compute_sorted_vector_of_distances_via_vector_sorting(where_to_cut);
352 this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
355 template <
typename F>
358 : where_to_cut(where_to_cut) {
359 std::vector<std::pair<double, double> > intervals;
360 if (dimension == std::numeric_limits<unsigned>::max()) {
361 intervals = read_persistence_intervals_in_one_dimension_from_file(filename);
363 intervals = read_persistence_intervals_in_one_dimension_from_file(filename, dimension);
365 this->intervals = intervals;
366 this->compute_sorted_vector_of_distances_via_heap(where_to_cut);
368 set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
371 template <
typename F>
375 std::cerr <<
"Here are the intervals : \n";
376 for (
size_t i = 0; i != this->intervals.size(); ++i) {
377 std::cerr << this->intervals[i].first <<
" , " << this->intervals[i].second << std::endl;
380 where_to_cut = std::min(
381 where_to_cut, (
size_t)(0.5 * this->intervals.size() * (this->intervals.size() - 1) + this->intervals.size()));
383 std::vector<double> heap(where_to_cut, std::numeric_limits<int>::max());
384 std::make_heap(heap.begin(), heap.end());
389 for (
size_t i = 0; i < this->intervals.size(); ++i) {
390 for (
size_t j = i + 1; j < this->intervals.size(); ++j) {
391 double value = std::min(
392 f(this->intervals[i], this->intervals[j]),
394 f(this->intervals[i], std::make_pair(0.5 * (this->intervals[i].first + this->intervals[i].second),
395 0.5 * (this->intervals[i].first + this->intervals[i].second))),
396 f(this->intervals[j], std::make_pair(0.5 * (this->intervals[j].first + this->intervals[j].second),
397 0.5 * (this->intervals[j].first + this->intervals[j].second)))));
400 std::cerr <<
"Value : " << value << std::endl;
401 std::cerr <<
"heap.front() : " << heap.front() << std::endl;
405 if (-value < heap.front()) {
407 std::cerr <<
"Replacing : " << heap.front() <<
" with : " << -value << std::endl;
411 std::pop_heap(heap.begin(), heap.end());
415 heap[where_to_cut - 1] = -value;
416 std::push_heap(heap.begin(), heap.end());
422 for (
size_t i = 0; i < this->intervals.size(); ++i) {
423 double value = f(this->intervals[i], std::make_pair(0.5 * (this->intervals[i].first + this->intervals[i].second),
424 0.5 * (this->intervals[i].first + this->intervals[i].second)));
425 if (-value < heap.front()) {
427 std::pop_heap(heap.begin(), heap.end());
431 heap[where_to_cut - 1] = -value;
432 std::push_heap(heap.begin(), heap.end());
436 std::sort_heap(heap.begin(), heap.end());
437 for (
size_t i = 0; i != heap.size(); ++i) {
438 if (heap[i] == std::numeric_limits<int>::max()) {
446 std::cerr <<
"This is the heap after all the operations :\n";
447 for (
size_t i = 0; i != heap.size(); ++i) {
448 std::cout << heap[i] <<
" ";
450 std::cout << std::endl;
453 this->sorted_vector_of_distances = heap;
456 template <
typename F>
458 std::vector<double> distances;
459 distances.reserve((
size_t)(0.5 * this->intervals.size() * (this->intervals.size() - 1) + this->intervals.size()));
464 for (
size_t i = 0; i < this->intervals.size(); ++i) {
467 f(this->intervals[i], std::make_pair(0.5 * (this->intervals[i].first + this->intervals[i].second),
468 0.5 * (this->intervals[i].first + this->intervals[i].second))));
469 for (
size_t j = i + 1; j < this->intervals.size(); ++j) {
470 double value = std::min(
471 f(this->intervals[i], this->intervals[j]),
473 f(this->intervals[i], std::make_pair(0.5 * (this->intervals[i].first + this->intervals[i].second),
474 0.5 * (this->intervals[i].first + this->intervals[i].second))),
475 f(this->intervals[j], std::make_pair(0.5 * (this->intervals[j].first + this->intervals[j].second),
476 0.5 * (this->intervals[j].first + this->intervals[j].second)))));
477 distances.push_back(value);
480 std::sort(distances.begin(), distances.end(), std::greater<double>());
481 if (distances.size() > where_to_cut) distances.resize(where_to_cut);
483 this->sorted_vector_of_distances = distances;
487 template <
typename F>
489 if ((
size_t)number_of_function > this->number_of_functions_for_projections_to_reals)
490 throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::project_to_R";
491 if (number_of_function < 0)
492 throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::project_to_R";
495 for (
size_t i = 0; i != (size_t)number_of_function; ++i) {
496 result += sorted_vector_of_distances[i];
501 template <
typename F>
503 if (to_average.size() == 0) {
508 size_t maximal_length_of_vector = 0;
509 for (
size_t i = 0; i != to_average.size(); ++i) {
510 if (to_average[i]->sorted_vector_of_distances.size() > maximal_length_of_vector) {
511 maximal_length_of_vector = to_average[i]->sorted_vector_of_distances.size();
515 std::vector<double> av(maximal_length_of_vector, 0);
516 for (
size_t i = 0; i != to_average.size(); ++i) {
517 for (
size_t j = 0; j != to_average[i]->sorted_vector_of_distances.size(); ++j) {
518 av[j] += to_average[i]->sorted_vector_of_distances[j];
522 for (
size_t i = 0; i != maximal_length_of_vector; ++i) {
523 av[i] /=
static_cast<double>(to_average.size());
525 this->sorted_vector_of_distances = av;
526 this->where_to_cut = av.size();
529 template <
typename F>
534 std::cerr <<
"Entering double Vector_distances_in_diagram<F>::distance( const Abs_Topological_data_with_distances* " 535 "second , double power ) procedure \n";
536 std::cerr <<
"Power : " << power << std::endl;
537 std::cerr <<
"This : " << *
this << std::endl;
538 std::cerr <<
"second : " << second_ << std::endl;
542 for (
size_t i = 0; i != std::min(this->sorted_vector_of_distances.size(), second_.sorted_vector_of_distances.size());
546 std::cerr <<
"|" << this->sorted_vector_of_distances[i] <<
" - " << second_.sorted_vector_of_distances[i]
547 <<
" | : " << fabs(this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i])
550 result += fabs(this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i]);
552 if (power < std::numeric_limits<double>::max()) {
553 result += std::pow(fabs(this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i]), power);
556 if (result < fabs(this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i]))
557 result = fabs(this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i]);
560 std::cerr <<
"| " << this->sorted_vector_of_distances[i] <<
" - " << second_.sorted_vector_of_distances[i]
561 <<
" : " << fabs(this->sorted_vector_of_distances[i] - second_.sorted_vector_of_distances[i])
566 if (this->sorted_vector_of_distances.size() != second_.sorted_vector_of_distances.size()) {
567 if (this->sorted_vector_of_distances.size() > second_.sorted_vector_of_distances.size()) {
568 for (
size_t i = second_.sorted_vector_of_distances.size(); i != this->sorted_vector_of_distances.size(); ++i) {
569 result += fabs(this->sorted_vector_of_distances[i]);
573 for (
size_t i = this->sorted_vector_of_distances.size(); i != second_.sorted_vector_of_distances.size(); ++i) {
574 result += fabs(second_.sorted_vector_of_distances[i]);
580 result = std::pow(result, (1.0 / power));
585 template <
typename F>
587 if ((
size_t)number_of_function > this->number_of_functions_for_vectorization)
588 throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::vectorize";
589 if (number_of_function < 0)
throw "Wrong index of a function in a method Vector_distances_in_diagram<F>::vectorize";
591 std::vector<double> result(std::min((
size_t)number_of_function, this->sorted_vector_of_distances.size()));
592 for (
size_t i = 0; i != std::min((
size_t)number_of_function, this->sorted_vector_of_distances.size()); ++i) {
593 result[i] = this->sorted_vector_of_distances[i];
598 template <
typename F>
603 for (
size_t i = 0; i != this->sorted_vector_of_distances.size(); ++i) {
604 out << this->sorted_vector_of_distances[i] <<
" ";
610 template <
typename F>
616 std::cerr <<
"The file : " << filename <<
" do not exist. The program will now terminate \n";
617 throw "The persistence landscape file do not exist. The program will now terminate \n";
621 while (in >> number) {
622 this->sorted_vector_of_distances.push_back(number);
627 template <
typename F>
631 i != std::min(this->sorted_vector_of_distances.size(), second_vector.sorted_vector_of_distances.size()); ++i) {
632 result += this->sorted_vector_of_distances[i] * second_vector.sorted_vector_of_distances[i];
640 #endif // PERSISTENCE_VECTORS_H_ double vector_in_position(size_t position) const
Definition: Persistence_vectors.h:100
Vector_distances_in_diagram operator/=(double x)
Definition: Persistence_vectors.h:316
void load_from_file(const char *filename)
Definition: Persistence_vectors.h:611
void print_to_file(const char *filename) const
Definition: Persistence_vectors.h:119
std::vector< double > output_for_visualization() const
Definition: Persistence_vectors.h:185
std::pair< double, double > get_x_range() const
Definition: Persistence_vectors.h:211
A class implementing persistence vectors.
Definition: Persistence_vectors.h:66
Definition: SimplicialComplexForAlpha.h:26
Vector_distances_in_diagram operator*(double scalar)
Definition: Persistence_vectors.h:291
std::vector< double > vectorize(int number_of_function) const
Definition: Persistence_vectors.h:586
size_t number_of_vectorize_functions() const
Definition: Persistence_vectors.h:162
Vector_distances_in_diagram operator*=(double x)
Definition: Persistence_vectors.h:309
void write_to_file(const char *filename) const
Definition: Persistence_vectors.h:599
size_t size() const
Definition: Persistence_vectors.h:109
void plot(const char *filename) const
Definition: Persistence_vectors.h:190
size_t number_of_projections_to_R() const
Definition: Persistence_vectors.h:152
friend Vector_distances_in_diagram operator*(double scalar, const Vector_distances_in_diagram &A)
Definition: Persistence_vectors.h:279
void compute_average(const std::vector< Vector_distances_in_diagram *> &to_average)
Definition: Persistence_vectors.h:502
Global distance functions.
bool operator==(const Vector_distances_in_diagram &second) const
Definition: Persistence_vectors.h:129
Vector_distances_in_diagram operator+=(const Vector_distances_in_diagram &rhs)
Definition: Persistence_vectors.h:295
std::pair< double, double > get_y_range() const
Definition: Persistence_vectors.h:216
friend Vector_distances_in_diagram operator*(const Vector_distances_in_diagram &A, double scalar)
Definition: Persistence_vectors.h:285
Vector_distances_in_diagram()
Definition: Persistence_vectors.h:71
double compute_scalar_product(const Vector_distances_in_diagram &second) const
Definition: Persistence_vectors.h:628
double project_to_R(int number_of_function) const
Definition: Persistence_vectors.h:488
double distance(const Vector_distances_in_diagram &second, double power=1) const
Definition: Persistence_vectors.h:530
friend Vector_distances_in_diagram operator-(const Vector_distances_in_diagram &first, const Vector_distances_in_diagram &second)
Definition: Persistence_vectors.h:272
Vector_distances_in_diagram multiply_by_scalar(double scalar) const
Definition: Persistence_vectors.h:253
friend Vector_distances_in_diagram operator+(const Vector_distances_in_diagram &first, const Vector_distances_in_diagram &second)
Definition: Persistence_vectors.h:265
Vector_distances_in_diagram operator-=(const Vector_distances_in_diagram &rhs)
Definition: Persistence_vectors.h:302