Loading...
Searching...
No Matches
Persistence_landscape.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Pawel Dlotko
4 *
5 * Copyright (C) 2016 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
11#ifndef PERSISTENCE_LANDSCAPE_H_
12#define PERSISTENCE_LANDSCAPE_H_
13
14// gudhi include
15#include <gudhi/read_persistence_from_file.h>
16#include <gudhi/common_persistence_representations.h>
17
18// standard include
19#include <cmath>
20#include <iostream>
21#include <vector>
22#include <limits>
23#include <fstream>
24#include <sstream>
25#include <algorithm>
26#include <string>
27#include <utility>
28#include <functional>
29
30namespace Gudhi {
31namespace Persistence_representations {
32
33// pre declaration
35template <typename operation>
36Persistence_landscape operation_on_pair_of_landscapes(const Persistence_landscape& land1,
37 const Persistence_landscape& land2);
38
61 public:
65 Persistence_landscape() { this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals(); }
66
70 Persistence_landscape(const std::vector<std::pair<double, double> >& p,
71 size_t number_of_levels = std::numeric_limits<size_t>::max());
72
78 Persistence_landscape(const char* filename, size_t dimension = std::numeric_limits<unsigned>::max(),
79 size_t number_of_levels = std::numeric_limits<size_t>::max());
80
84 void load_landscape_from_file(const char* filename);
85
89 void print_to_file(const char* filename) const;
90
95 double compute_integral_of_landscape() const;
96
100 double compute_integral_of_a_level_of_a_landscape(size_t level) const;
101
106 double compute_integral_of_landscape(double p) const; // this function compute integral of p-th power of landscape.
107
113 double compute_value_at_a_given_point(unsigned level, double x) const;
114
120 friend std::ostream& operator<<(std::ostream& out, Persistence_landscape& land);
121
122 template <typename operation>
123 friend Persistence_landscape operation_on_pair_of_landscapes(const Persistence_landscape& land1,
124 const Persistence_landscape& land2);
125
130 const Persistence_landscape& land2) {
131 return operation_on_pair_of_landscapes<std::plus<double> >(land1, land2);
132 }
133
138 const Persistence_landscape& land2) {
139 return operation_on_pair_of_landscapes<std::minus<double> >(land1, land2);
140 }
141
146 return add_two_landscapes(first, second);
147 }
148
153 return subtract_two_landscapes(first, second);
154 }
155
159 friend Persistence_landscape operator*(const Persistence_landscape& first, double con) {
160 return first.multiply_lanscape_by_real_number_not_overwrite(con);
161 }
162
166 friend Persistence_landscape operator*(double con, const Persistence_landscape& first) {
167 return first.multiply_lanscape_by_real_number_not_overwrite(con);
168 }
169
174 *this = *this + rhs;
175 return *this;
176 }
177
182 *this = *this - rhs;
183 return *this;
184 }
185
191 *this = *this * x;
192 return *this;
193 }
194
199 if (x == 0) throw("In operator /=, division by 0. Program terminated.");
200 *this = *this * (1 / x);
201 return *this;
202 }
203
207 bool operator==(const Persistence_landscape& rhs) const;
208
212 bool operator!=(const Persistence_landscape& rhs) const { return !((*this) == rhs); }
213
217 double compute_maximum() const {
218 double maxValue = 0;
219 if (this->land.size()) {
220 maxValue = -std::numeric_limits<int>::max();
221 for (size_t i = 0; i != this->land[0].size(); ++i) {
222 if (this->land[0][i].second > maxValue) maxValue = this->land[0][i].second;
223 }
224 }
225 return maxValue;
226 }
227
231 double compute_minimum() const {
232 double minValue = 0;
233 if (this->land.size()) {
234 minValue = std::numeric_limits<int>::max();
235 for (size_t i = 0; i != this->land[0].size(); ++i) {
236 if (this->land[0][i].second < minValue) minValue = this->land[0][i].second;
237 }
238 }
239 return minValue;
240 }
241
245 double compute_norm_of_landscape(double i) {
247 if (i < std::numeric_limits<double>::max()) {
248 return compute_distance_of_landscapes(*this, l, i);
249 } else {
251 }
252 }
253
257 double operator()(unsigned level, double x) const { return this->compute_value_at_a_given_point(level, x); }
258
263 const Persistence_landscape& second);
264
268 friend double compute_distance_of_landscapes(const Persistence_landscape& first, const Persistence_landscape& second,
269 double p);
270
278
279 Persistence_landscape* new_abs();
280
284 size_t size() const { return this->land.size(); }
285
289 double find_max(unsigned lambda) const;
290
294 friend double compute_inner_product(const Persistence_landscape& l1, const Persistence_landscape& l2);
295
296 // Implementations of functions for various concepts.
297
306 double project_to_R(int number_of_function) const {
307 return this->compute_integral_of_a_level_of_a_landscape((size_t)number_of_function);
308 }
309
314 size_t number_of_projections_to_R() const { return this->number_of_functions_for_projections_to_reals; }
315
320 std::vector<double> vectorize(int number_of_function) const {
321 // TODO(PD) think of something smarter over here
322 std::vector<double> v;
323 if ((size_t)number_of_function > this->land.size()) {
324 return v;
325 }
326 v.reserve(this->land[number_of_function].size());
327 for (size_t i = 0; i != this->land[number_of_function].size(); ++i) {
328 v.push_back(this->land[number_of_function][i].second);
329 }
330 return v;
331 }
336 size_t number_of_vectorize_functions() const { return this->number_of_functions_for_vectorization; }
337
342 void compute_average(const std::vector<Persistence_landscape*>& to_average) {
343 bool dbg = false;
344
345 if (dbg) {
346 std::clog << "to_average.size() : " << to_average.size() << std::endl;
347 }
348
349 std::vector<Persistence_landscape*> nextLevelMerge(to_average.size());
350 for (size_t i = 0; i != to_average.size(); ++i) {
351 nextLevelMerge[i] = to_average[i];
352 }
353 bool is_this_first_level = true; // in the loop, we will create dynamically a number of intermediate complexes. We
354 // have to clean that up, but we cannot erase the initial landscapes we have
355 // to average. In this case, we simply check if the nextLevelMerge are the input landscapes or the ones created in
356 // that loop by using this extra variable.
357
358 while (nextLevelMerge.size() != 1) {
359 if (dbg) {
360 std::clog << "nextLevelMerge.size() : " << nextLevelMerge.size() << std::endl;
361 }
362 std::vector<Persistence_landscape*> nextNextLevelMerge;
363 nextNextLevelMerge.reserve(to_average.size());
364 for (size_t i = 0; i < nextLevelMerge.size(); i = i + 2) {
365 if (dbg) {
366 std::clog << "i : " << i << std::endl;
367 }
369 if (i + 1 != nextLevelMerge.size()) {
370 (*l) = (*nextLevelMerge[i]) + (*nextLevelMerge[i + 1]);
371 } else {
372 (*l) = *nextLevelMerge[i];
373 }
374 nextNextLevelMerge.push_back(l);
375 }
376 if (dbg) {
377 std::clog << "After this iteration \n";
378 getchar();
379 }
380
381 if (!is_this_first_level) {
382 // deallocate the memory if the vector nextLevelMerge do not consist of the initial landscapes
383 for (size_t i = 0; i != nextLevelMerge.size(); ++i) {
384 delete nextLevelMerge[i];
385 }
386 }
387 is_this_first_level = false;
388 nextLevelMerge.swap(nextNextLevelMerge);
389 }
390 (*this) = (*nextLevelMerge[0]);
391 (*this) *= 1 / static_cast<double>(to_average.size());
392 }
393
400 double distance(const Persistence_landscape& second, double power = 1) const {
401 if (power < std::numeric_limits<double>::max()) {
402 return compute_distance_of_landscapes(*this, second, power);
403 } else {
404 return compute_max_norm_distance_of_landscapes(*this, second);
405 }
406 }
407
413 double compute_scalar_product(const Persistence_landscape& second) const {
414 return compute_inner_product((*this), second);
415 }
416 // end of implementation of functions needed for concepts.
417
422 std::pair<double, double> get_y_range(size_t level = 0) const {
423 std::pair<double, double> result;
424 if (level < this->land.size()) {
425 double maxx = this->compute_maximum();
426 double minn = this->compute_minimum();
427 result = std::make_pair(minn, maxx);
428 } else {
429 result = std::make_pair(0, 0);
430 }
431 return result;
432 }
433
434 // a function used to create a gnuplot script for visualization of landscapes
435 void plot(const char* filename, double xRangeBegin = std::numeric_limits<double>::max(),
436 double xRangeEnd = std::numeric_limits<double>::max(),
437 double yRangeBegin = std::numeric_limits<double>::max(),
438 double yRangeEnd = std::numeric_limits<double>::max(), int from = std::numeric_limits<int>::max(),
439 int to = std::numeric_limits<int>::max());
440
441 protected:
442 std::vector<std::vector<std::pair<double, double> > > land;
443 size_t number_of_functions_for_vectorization;
444 size_t number_of_functions_for_projections_to_reals;
445
446 void construct_persistence_landscape_from_barcode(const std::vector<std::pair<double, double> >& p,
447 size_t number_of_levels = std::numeric_limits<size_t>::max());
448 Persistence_landscape multiply_lanscape_by_real_number_not_overwrite(double x) const;
449 void multiply_lanscape_by_real_number_overwrite(double x);
450 friend double compute_maximal_distance_non_symmetric(const Persistence_landscape& pl1,
451 const Persistence_landscape& pl2);
452
453 void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals() {
454 // warning, this function can be only called after filling in the intervals vector.
455 this->number_of_functions_for_vectorization = this->land.size();
456 this->number_of_functions_for_projections_to_reals = this->land.size();
457 }
458};
459
460Persistence_landscape::Persistence_landscape(const char* filename, size_t dimension, size_t number_of_levels) {
461 std::vector<std::pair<double, double> > barcode;
462 if (dimension < std::numeric_limits<double>::max()) {
463 barcode = read_persistence_intervals_in_one_dimension_from_file(filename, dimension);
464 } else {
465 barcode = read_persistence_intervals_in_one_dimension_from_file(filename);
466 }
467 this->construct_persistence_landscape_from_barcode(barcode, number_of_levels);
468 this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
469}
470
471bool operatorEqualDbg = false;
473 if (this->land.size() != rhs.land.size()) {
474 if (operatorEqualDbg) std::clog << "1\n";
475 return false;
476 }
477 for (size_t level = 0; level != this->land.size(); ++level) {
478 if (this->land[level].size() != rhs.land[level].size()) {
479 if (operatorEqualDbg) std::clog << "this->land[level].size() : " << this->land[level].size() << "\n";
480 if (operatorEqualDbg) std::clog << "rhs.land[level].size() : " << rhs.land[level].size() << "\n";
481 if (operatorEqualDbg) std::clog << "2\n";
482 return false;
483 }
484 for (size_t i = 0; i != this->land[level].size(); ++i) {
485 if (!(almost_equal(this->land[level][i].first, rhs.land[level][i].first) &&
486 almost_equal(this->land[level][i].second, rhs.land[level][i].second))) {
487 if (operatorEqualDbg)
488 std::clog << "this->land[level][i] : " << this->land[level][i].first << " " << this->land[level][i].second
489 << "\n";
490 if (operatorEqualDbg)
491 std::clog << "rhs.land[level][i] : " << rhs.land[level][i].first << " " << rhs.land[level][i].second << "\n";
492 if (operatorEqualDbg) std::clog << "3\n";
493 return false;
494 }
495 }
496 }
497 return true;
498}
499
500Persistence_landscape::Persistence_landscape(const std::vector<std::pair<double, double> >& p,
501 size_t number_of_levels) {
502 this->construct_persistence_landscape_from_barcode(p, number_of_levels);
503 this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
504}
505
506void Persistence_landscape::construct_persistence_landscape_from_barcode(
507 const std::vector<std::pair<double, double> >& p, size_t number_of_levels) {
508 bool dbg = false;
509 if (dbg) {
510 std::clog << "Persistence_landscape::Persistence_landscape( const std::vector< std::pair< double , double > >& p )"
511 << std::endl;
512 }
513
514 // this is a general algorithm to construct persistence landscapes.
515 std::vector<std::pair<double, double> > bars;
516 bars.insert(bars.begin(), p.begin(), p.end());
517 std::sort(bars.begin(), bars.end(), compare_points_sorting);
518
519 if (dbg) {
520 std::clog << "Bars : \n";
521 for (size_t i = 0; i != bars.size(); ++i) {
522 std::clog << bars[i].first << " " << bars[i].second << "\n";
523 }
524 getchar();
525 }
526
527 std::vector<std::pair<double, double> > characteristicPoints(p.size());
528 for (size_t i = 0; i != bars.size(); ++i) {
529 characteristicPoints[i] =
530 std::make_pair((bars[i].first + bars[i].second) / 2.0, (bars[i].second - bars[i].first) / 2.0);
531 }
532 std::vector<std::vector<std::pair<double, double> > > Persistence_landscape;
533 size_t number_of_levels_in_the_landscape = 0;
534 while (!characteristicPoints.empty()) {
535 if (dbg) {
536 for (size_t i = 0; i != characteristicPoints.size(); ++i) {
537 std::clog << "(" << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n";
538 }
539 std::cin.ignore();
540 }
541
542 std::vector<std::pair<double, double> > lambda_n;
543 lambda_n.push_back(std::make_pair(-std::numeric_limits<int>::max(), 0));
544 lambda_n.push_back(std::make_pair(minus_length(characteristicPoints[0]), 0));
545 lambda_n.push_back(characteristicPoints[0]);
546
547 if (dbg) {
548 std::clog << "1 Adding to lambda_n : (" << -std::numeric_limits<int>::max() << " " << 0 << ") , ("
549 << minus_length(characteristicPoints[0]) << " " << 0 << ") , (" << characteristicPoints[0].first << " "
550 << characteristicPoints[0].second << ") \n";
551 }
552
553 size_t i = 1;
554 std::vector<std::pair<double, double> > newCharacteristicPoints;
555 while (i < characteristicPoints.size()) {
556 size_t p = 1;
557 if ((minus_length(characteristicPoints[i]) >= minus_length(lambda_n[lambda_n.size() - 1])) &&
558 (birth_plus_deaths(characteristicPoints[i]) > birth_plus_deaths(lambda_n[lambda_n.size() - 1]))) {
559 if (minus_length(characteristicPoints[i]) < birth_plus_deaths(lambda_n[lambda_n.size() - 1])) {
560 std::pair<double, double> point = std::make_pair(
561 (minus_length(characteristicPoints[i]) + birth_plus_deaths(lambda_n[lambda_n.size() - 1])) / 2,
562 (birth_plus_deaths(lambda_n[lambda_n.size() - 1]) - minus_length(characteristicPoints[i])) / 2);
563 lambda_n.push_back(point);
564 if (dbg) {
565 std::clog << "2 Adding to lambda_n : (" << point.first << " " << point.second << ")\n";
566 }
567
568 if (dbg) {
569 std::clog << "characteristicPoints[i+p] : " << characteristicPoints[i + p].first << " "
570 << characteristicPoints[i + p].second << "\n";
571 std::clog << "point : " << point.first << " " << point.second << "\n";
572 getchar();
573 }
574
575 while ((i + p < characteristicPoints.size()) &&
576 (almost_equal(minus_length(point), minus_length(characteristicPoints[i + p]))) &&
577 (birth_plus_deaths(point) <= birth_plus_deaths(characteristicPoints[i + p]))) {
578 newCharacteristicPoints.push_back(characteristicPoints[i + p]);
579 if (dbg) {
580 std::clog << "3.5 Adding to newCharacteristicPoints : (" << characteristicPoints[i + p].first << " "
581 << characteristicPoints[i + p].second << ")\n";
582 getchar();
583 }
584 ++p;
585 }
586
587 newCharacteristicPoints.push_back(point);
588 if (dbg) {
589 std::clog << "4 Adding to newCharacteristicPoints : (" << point.first << " " << point.second << ")\n";
590 }
591
592 while ((i + p < characteristicPoints.size()) &&
593 (minus_length(point) <= minus_length(characteristicPoints[i + p])) &&
594 (birth_plus_deaths(point) >= birth_plus_deaths(characteristicPoints[i + p]))) {
595 newCharacteristicPoints.push_back(characteristicPoints[i + p]);
596 if (dbg) {
597 std::clog << "characteristicPoints[i+p] : " << characteristicPoints[i + p].first << " "
598 << characteristicPoints[i + p].second << "\n";
599 std::clog << "point : " << point.first << " " << point.second << "\n";
600 std::clog << "characteristicPoints[i+p] birth and death : " << minus_length(characteristicPoints[i + p])
601 << " , " << birth_plus_deaths(characteristicPoints[i + p]) << "\n";
602 std::clog << "point birth and death : " << minus_length(point) << " , " << birth_plus_deaths(point)
603 << "\n";
604
605 std::clog << "3 Adding to newCharacteristicPoints : (" << characteristicPoints[i + p].first << " "
606 << characteristicPoints[i + p].second << ")\n";
607 getchar();
608 }
609 ++p;
610 }
611
612 } else {
613 lambda_n.push_back(std::make_pair(birth_plus_deaths(lambda_n[lambda_n.size() - 1]), 0));
614 lambda_n.push_back(std::make_pair(minus_length(characteristicPoints[i]), 0));
615 if (dbg) {
616 std::clog << "5 Adding to lambda_n : (" << birth_plus_deaths(lambda_n[lambda_n.size() - 1]) << " " << 0
617 << ")\n";
618 std::clog << "5 Adding to lambda_n : (" << minus_length(characteristicPoints[i]) << " " << 0 << ")\n";
619 }
620 }
621 lambda_n.push_back(characteristicPoints[i]);
622 if (dbg) {
623 std::clog << "6 Adding to lambda_n : (" << characteristicPoints[i].first << " "
624 << characteristicPoints[i].second << ")\n";
625 }
626 } else {
627 newCharacteristicPoints.push_back(characteristicPoints[i]);
628 if (dbg) {
629 std::clog << "7 Adding to newCharacteristicPoints : (" << characteristicPoints[i].first << " "
630 << characteristicPoints[i].second << ")\n";
631 }
632 }
633 i = i + p;
634 }
635 lambda_n.push_back(std::make_pair(birth_plus_deaths(lambda_n[lambda_n.size() - 1]), 0));
636 lambda_n.push_back(std::make_pair(std::numeric_limits<int>::max(), 0));
637
638 characteristicPoints = newCharacteristicPoints;
639
640 lambda_n.erase(std::unique(lambda_n.begin(), lambda_n.end()), lambda_n.end());
641 this->land.push_back(lambda_n);
642
643 ++number_of_levels_in_the_landscape;
644 if (number_of_levels == number_of_levels_in_the_landscape) {
645 break;
646 }
647 }
648}
649
650// this function find maximum of lambda_n
651double Persistence_landscape::find_max(unsigned lambda) const {
652 if (this->land.size() < lambda) return 0;
653 double maximum = -std::numeric_limits<int>::max();
654 for (size_t i = 0; i != this->land[lambda].size(); ++i) {
655 if (this->land[lambda][i].second > maximum) maximum = this->land[lambda][i].second;
656 }
657 return maximum;
658}
659
661 double result = 0;
662 for (size_t i = 0; i != this->land.size(); ++i) {
663 for (size_t nr = 2; nr != this->land[i].size() - 1; ++nr) {
664 // it suffices to compute every planar integral and then sum them up for each lambda_n
665 result += 0.5 * (this->land[i][nr].first - this->land[i][nr - 1].first) *
666 (this->land[i][nr].second + this->land[i][nr - 1].second);
667 }
668 }
669 return result;
670}
671
673 double result = 0;
674 if (level >= this->land.size()) {
675 // this landscape function is constantly equal 0, so is the integral.
676 return result;
677 }
678 // also negative landscapes are assumed to be zero.
679 if (level < 0) return 0;
680
681 for (size_t nr = 2; nr != this->land[level].size() - 1; ++nr) {
682 // it suffices to compute every planar integral and then sum them up for each lambda_n
683 result += 0.5 * (this->land[level][nr].first - this->land[level][nr - 1].first) *
684 (this->land[level][nr].second + this->land[level][nr - 1].second);
685 }
686
687 return result;
688}
689
691 bool dbg = false;
692 double result = 0;
693 for (size_t i = 0; i != this->land.size(); ++i) {
694 for (size_t nr = 2; nr != this->land[i].size() - 1; ++nr) {
695 if (dbg) std::clog << "nr : " << nr << "\n";
696 // In this interval, the landscape has a form f(x) = ax+b. We want to compute integral of (ax+b)^p = 1/a *
697 // (ax+b)^{p+1}/(p+1)
698 std::pair<double, double> coef = compute_parameters_of_a_line(this->land[i][nr], this->land[i][nr - 1]);
699 double a = coef.first;
700 double b = coef.second;
701
702 if (dbg)
703 std::clog << "(" << this->land[i][nr].first << "," << this->land[i][nr].second << ") , "
704 << this->land[i][nr - 1].first << "," << this->land[i][nr].second << ")" << std::endl;
705 if (this->land[i][nr].first == this->land[i][nr - 1].first) continue;
706 if (a != 0) {
707 result += 1 / (a * (p + 1)) *
708 (pow((a * this->land[i][nr].first + b), p + 1) - pow((a * this->land[i][nr - 1].first + b), p + 1));
709 } else {
710 result += (this->land[i][nr].first - this->land[i][nr - 1].first) * (pow(this->land[i][nr].second, p));
711 }
712 if (dbg) {
713 std::clog << "a : " << a << " , b : " << b << std::endl;
714 std::clog << "result : " << result << std::endl;
715 }
716 }
717 }
718 return result;
719}
720
721// this is O(log(n)) algorithm, where n is number of points in this->land.
722double Persistence_landscape::compute_value_at_a_given_point(unsigned level, double x) const {
723 bool compute_value_at_a_given_pointDbg = false;
724 // in such a case lambda_level = 0.
725 if (level >= this->land.size()) return 0;
726
727 // we know that the points in this->land[level] are ordered according to x coordinate. Therefore, we can find the
728 // point by using bisection:
729 unsigned coordBegin = 1;
730 unsigned coordEnd = this->land[level].size() - 2;
731
732 if (compute_value_at_a_given_pointDbg) {
733 std::clog << "Here \n";
734 std::clog << "x : " << x << "\n";
735 std::clog << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n";
736 std::clog << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n";
737 }
738
739 // in this case x is outside the support of the landscape, therefore the value of the landscape is 0.
740 if (x <= this->land[level][coordBegin].first) return 0;
741 if (x >= this->land[level][coordEnd].first) return 0;
742
743 if (compute_value_at_a_given_pointDbg) std::clog << "Entering to the while loop \n";
744
745 while (coordBegin + 1 != coordEnd) {
746 if (compute_value_at_a_given_pointDbg) {
747 std::clog << "coordBegin : " << coordBegin << "\n";
748 std::clog << "coordEnd : " << coordEnd << "\n";
749 std::clog << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n";
750 std::clog << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n";
751 }
752
753 unsigned newCord = (unsigned)floor((coordEnd + coordBegin) / 2.0);
754
755 if (compute_value_at_a_given_pointDbg) {
756 std::clog << "newCord : " << newCord << "\n";
757 std::clog << "this->land[level][newCord].first : " << this->land[level][newCord].first << "\n";
758 std::cin.ignore();
759 }
760
761 if (this->land[level][newCord].first <= x) {
762 coordBegin = newCord;
763 if (this->land[level][newCord].first == x) return this->land[level][newCord].second;
764 } else {
765 coordEnd = newCord;
766 }
767 }
768
769 if (compute_value_at_a_given_pointDbg) {
770 std::clog << "x : " << x << " is between : " << this->land[level][coordBegin].first << " a "
771 << this->land[level][coordEnd].first << "\n";
772 std::clog << "the y coords are : " << this->land[level][coordBegin].second << " a "
773 << this->land[level][coordEnd].second << "\n";
774 std::clog << "coordBegin : " << coordBegin << "\n";
775 std::clog << "coordEnd : " << coordEnd << "\n";
776 std::cin.ignore();
777 }
778 return function_value(this->land[level][coordBegin], this->land[level][coordEnd], x);
779}
780
781std::ostream& operator<<(std::ostream& out, Persistence_landscape& land) {
782 for (size_t level = 0; level != land.land.size(); ++level) {
783 out << "Lambda_" << level << ":" << std::endl;
784 for (size_t i = 0; i != land.land[level].size(); ++i) {
785 if (land.land[level][i].first == -std::numeric_limits<int>::max()) {
786 out << "-inf";
787 } else {
788 if (land.land[level][i].first == std::numeric_limits<int>::max()) {
789 out << "+inf";
790 } else {
791 out << land.land[level][i].first;
792 }
793 }
794 out << " , " << land.land[level][i].second << std::endl;
795 }
796 }
797 return out;
798}
799
800void Persistence_landscape::multiply_lanscape_by_real_number_overwrite(double x) {
801 for (size_t dim = 0; dim != this->land.size(); ++dim) {
802 for (size_t i = 0; i != this->land[dim].size(); ++i) {
803 this->land[dim][i].second *= x;
804 }
805 }
806}
807
808bool AbsDbg = false;
811 for (size_t level = 0; level != this->land.size(); ++level) {
812 if (AbsDbg) {
813 std::clog << "level: " << level << std::endl;
814 }
815 std::vector<std::pair<double, double> > lambda_n;
816 lambda_n.push_back(std::make_pair(-std::numeric_limits<int>::max(), 0));
817 for (size_t i = 1; i != this->land[level].size(); ++i) {
818 if (AbsDbg) {
819 std::clog << "this->land[" << level << "][" << i << "] : " << this->land[level][i].first << " "
820 << this->land[level][i].second << std::endl;
821 }
822 // if a line segment between this->land[level][i-1] and this->land[level][i] crosses the x-axis, then we have to
823 // add one landscape point t o result
824 if ((this->land[level][i - 1].second) * (this->land[level][i].second) < 0) {
825 double zero =
826 find_zero_of_a_line_segment_between_those_two_points(this->land[level][i - 1], this->land[level][i]);
827
828 lambda_n.push_back(std::make_pair(zero, 0));
829 lambda_n.push_back(std::make_pair(this->land[level][i].first, fabs(this->land[level][i].second)));
830 if (AbsDbg) {
831 std::clog << "Adding pair : (" << zero << ",0)" << std::endl;
832 std::clog << "In the same step adding pair : (" << this->land[level][i].first << ","
833 << fabs(this->land[level][i].second) << ") " << std::endl;
834 std::cin.ignore();
835 }
836 } else {
837 lambda_n.push_back(std::make_pair(this->land[level][i].first, fabs(this->land[level][i].second)));
838 if (AbsDbg) {
839 std::clog << "Adding pair : (" << this->land[level][i].first << "," << fabs(this->land[level][i].second)
840 << ") " << std::endl;
841 std::cin.ignore();
842 }
843 }
844 }
845 result.land.push_back(lambda_n);
846 }
847 return result;
848}
849
850Persistence_landscape* Persistence_landscape::new_abs() {
851 Persistence_landscape* result = new Persistence_landscape(*this);
852 for (size_t level = 0; level != this->land.size(); ++level) {
853 if (AbsDbg) {
854 std::clog << "level: " << level << std::endl;
855 }
856 std::vector<std::pair<double, double> > lambda_n;
857 lambda_n.push_back(std::make_pair(-std::numeric_limits<int>::max(), 0));
858 for (size_t i = 1; i != this->land[level].size(); ++i) {
859 if (AbsDbg) {
860 std::clog << "this->land[" << level << "][" << i << "] : " << this->land[level][i].first << " "
861 << this->land[level][i].second << std::endl;
862 }
863 // if a line segment between this->land[level][i-1] and this->land[level][i] crosses the x-axis, then we have to
864 // add one landscape point t o result
865 if ((this->land[level][i - 1].second) * (this->land[level][i].second) < 0) {
866 double zero =
867 find_zero_of_a_line_segment_between_those_two_points(this->land[level][i - 1], this->land[level][i]);
868
869 lambda_n.push_back(std::make_pair(zero, 0));
870 lambda_n.push_back(std::make_pair(this->land[level][i].first, fabs(this->land[level][i].second)));
871 if (AbsDbg) {
872 std::clog << "Adding pair : (" << zero << ",0)" << std::endl;
873 std::clog << "In the same step adding pair : (" << this->land[level][i].first << ","
874 << fabs(this->land[level][i].second) << ") " << std::endl;
875 std::cin.ignore();
876 }
877 } else {
878 lambda_n.push_back(std::make_pair(this->land[level][i].first, fabs(this->land[level][i].second)));
879 if (AbsDbg) {
880 std::clog << "Adding pair : (" << this->land[level][i].first << "," << fabs(this->land[level][i].second)
881 << ") " << std::endl;
882 std::cin.ignore();
883 }
884 }
885 }
886 result->land.push_back(lambda_n);
887 }
888 return result;
889}
890
891Persistence_landscape Persistence_landscape::multiply_lanscape_by_real_number_not_overwrite(double x) const {
892 std::vector<std::vector<std::pair<double, double> > > result(this->land.size());
893 for (size_t dim = 0; dim != this->land.size(); ++dim) {
894 std::vector<std::pair<double, double> > lambda_dim(this->land[dim].size());
895 for (size_t i = 0; i != this->land[dim].size(); ++i) {
896 lambda_dim[i] = std::make_pair(this->land[dim][i].first, x * this->land[dim][i].second);
897 }
898 result[dim] = lambda_dim;
899 }
901 // CHANGE
902 // res.land = result;
903 res.land.swap(result);
904 return res;
905} // multiply_lanscape_by_real_number_overwrite
906
907void Persistence_landscape::print_to_file(const char* filename) const {
908 std::ofstream write;
909 write.open(filename);
910 for (size_t dim = 0; dim != this->land.size(); ++dim) {
911 write << "#lambda_" << dim << std::endl;
912 for (size_t i = 1; i != this->land[dim].size() - 1; ++i) {
913 write << this->land[dim][i].first << " " << this->land[dim][i].second << std::endl;
914 }
915 }
916 write.close();
917}
918
920 bool dbg = false;
921 // removing the current content of the persistence landscape.
922 this->land.clear();
923
924 // this constructor reads persistence landscape form a file. This file have to be created by this software before head
925 std::ifstream in;
926 in.open(filename);
927 if (!in.good()) {
928 std::cerr << "The file : " << filename << " do not exist. The program will now terminate \n";
929 throw "The persistence landscape file do not exist. The program will now terminate \n";
930 }
931
932 std::string line;
933 std::vector<std::pair<double, double> > landscapeAtThisLevel;
934
935 bool isThisAFirsLine = true;
936 while (in.good()) {
937 getline(in, line);
938 if (!(line.length() == 0 || line[0] == '#')) {
939 std::stringstream lineSS;
940 lineSS << line;
941 double beginn, endd;
942 lineSS >> beginn;
943 lineSS >> endd;
944 landscapeAtThisLevel.push_back(std::make_pair(beginn, endd));
945 if (dbg) {
946 std::clog << "Reading a point : " << beginn << " , " << endd << std::endl;
947 }
948 } else {
949 if (dbg) {
950 std::clog << "IGNORE LINE\n";
951 getchar();
952 }
953 if (!isThisAFirsLine) {
954 landscapeAtThisLevel.push_back(std::make_pair(std::numeric_limits<int>::max(), 0));
955 this->land.push_back(landscapeAtThisLevel);
956 std::vector<std::pair<double, double> > newLevelOdLandscape;
957 landscapeAtThisLevel.swap(newLevelOdLandscape);
958 }
959 landscapeAtThisLevel.push_back(std::make_pair(-std::numeric_limits<int>::max(), 0));
960 isThisAFirsLine = false;
961 }
962 }
963 if (landscapeAtThisLevel.size() > 1) {
964 // seems that the last line of the file is not finished with the newline sign. We need to put what we have in
965 // landscapeAtThisLevel to the constructed landscape.
966 landscapeAtThisLevel.push_back(std::make_pair(std::numeric_limits<int>::max(), 0));
967 this->land.push_back(landscapeAtThisLevel);
968 }
969
970 in.close();
971}
972
973template <typename T>
974Persistence_landscape operation_on_pair_of_landscapes(const Persistence_landscape& land1,
975 const Persistence_landscape& land2) {
976 bool operation_on_pair_of_landscapesDBG = false;
977 if (operation_on_pair_of_landscapesDBG) {
978 std::clog << "operation_on_pair_of_landscapes\n";
979 std::cin.ignore();
980 }
982 std::vector<std::vector<std::pair<double, double> > > land(std::max(land1.land.size(), land2.land.size()));
983 result.land = land;
984 T oper;
985
986 if (operation_on_pair_of_landscapesDBG) {
987 for (size_t i = 0; i != std::min(land1.land.size(), land2.land.size()); ++i) {
988 std::clog << "land1.land[" << i << "].size() : " << land1.land[i].size() << std::endl;
989 std::clog << "land2.land[" << i << "].size() : " << land2.land[i].size() << std::endl;
990 }
991 getchar();
992 }
993
994 for (size_t i = 0; i != std::min(land1.land.size(), land2.land.size()); ++i) {
995 std::vector<std::pair<double, double> > lambda_n;
996 size_t p = 0;
997 size_t q = 0;
998 while ((p + 1 < land1.land[i].size()) && (q + 1 < land2.land[i].size())) {
999 if (operation_on_pair_of_landscapesDBG) {
1000 std::clog << "p : " << p << "\n";
1001 std::clog << "q : " << q << "\n";
1002 std::clog << "land1.land.size() : " << land1.land.size() << std::endl;
1003 std::clog << "land2.land.size() : " << land2.land.size() << std::endl;
1004 std::clog << "land1.land[" << i << "].size() : " << land1.land[i].size() << std::endl;
1005 std::clog << "land2.land[" << i << "].size() : " << land2.land[i].size() << std::endl;
1006 std::clog << "land1.land[i][p].first : " << land1.land[i][p].first << "\n";
1007 std::clog << "land2.land[i][q].first : " << land2.land[i][q].first << "\n";
1008 }
1009
1010 if (land1.land[i][p].first < land2.land[i][q].first) {
1011 if (operation_on_pair_of_landscapesDBG) {
1012 std::clog << "first \n";
1013 std::clog << " function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) : "
1014 << function_value(land2.land[i][q - 1], land2.land[i][q], land1.land[i][p].first) << "\n";
1015 }
1016 lambda_n.push_back(
1017 std::make_pair(land1.land[i][p].first,
1018 oper(static_cast<double>(land1.land[i][p].second),
1019 function_value(land2.land[i][q - 1], land2.land[i][q], land1.land[i][p].first))));
1020 ++p;
1021 continue;
1022 }
1023 if (land1.land[i][p].first > land2.land[i][q].first) {
1024 if (operation_on_pair_of_landscapesDBG) {
1025 std::clog << "Second \n";
1026 std::clog << "function_value(" << land1.land[i][p - 1].first << " " << land1.land[i][p - 1].second << " ,"
1027 << land1.land[i][p].first << " " << land1.land[i][p].second << ", " << land2.land[i][q].first
1028 << " ) : " << function_value(land1.land[i][p - 1], land1.land[i][p - 1], land2.land[i][q].first)
1029 << "\n";
1030 std::clog << "oper( " << function_value(land1.land[i][p], land1.land[i][p - 1], land2.land[i][q].first) << ","
1031 << land2.land[i][q].second << " : "
1032 << oper(land2.land[i][q].second,
1033 function_value(land1.land[i][p], land1.land[i][p - 1], land2.land[i][q].first))
1034 << "\n";
1035 }
1036 lambda_n.push_back(std::make_pair(
1037 land2.land[i][q].first, oper(function_value(land1.land[i][p], land1.land[i][p - 1], land2.land[i][q].first),
1038 land2.land[i][q].second)));
1039 ++q;
1040 continue;
1041 }
1042 if (land1.land[i][p].first == land2.land[i][q].first) {
1043 if (operation_on_pair_of_landscapesDBG) std::clog << "Third \n";
1044 lambda_n.push_back(
1045 std::make_pair(land2.land[i][q].first, oper(land1.land[i][p].second, land2.land[i][q].second)));
1046 ++p;
1047 ++q;
1048 }
1049 if (operation_on_pair_of_landscapesDBG) {
1050 std::clog << "Next iteration \n";
1051 }
1052 }
1053 while ((p + 1 < land1.land[i].size()) && (q + 1 >= land2.land[i].size())) {
1054 if (operation_on_pair_of_landscapesDBG) {
1055 std::clog << "New point : " << land1.land[i][p].first
1056 << " oper(land1.land[i][p].second,0) : " << oper(land1.land[i][p].second, 0) << std::endl;
1057 }
1058 lambda_n.push_back(std::make_pair(land1.land[i][p].first, oper(land1.land[i][p].second, 0)));
1059 ++p;
1060 }
1061 while ((p + 1 >= land1.land[i].size()) && (q + 1 < land2.land[i].size())) {
1062 if (operation_on_pair_of_landscapesDBG) {
1063 std::clog << "New point : " << land2.land[i][q].first
1064 << " oper(0,land2.land[i][q].second) : " << oper(0, land2.land[i][q].second) << std::endl;
1065 }
1066 lambda_n.push_back(std::make_pair(land2.land[i][q].first, oper(0, land2.land[i][q].second)));
1067 ++q;
1068 }
1069 lambda_n.push_back(std::make_pair(std::numeric_limits<int>::max(), 0));
1070 // CHANGE
1071 // result.land[i] = lambda_n;
1072 result.land[i].swap(lambda_n);
1073 }
1074 if (land1.land.size() > std::min(land1.land.size(), land2.land.size())) {
1075 if (operation_on_pair_of_landscapesDBG) {
1076 std::clog << "land1.land.size() > std::min( land1.land.size() , land2.land.size() )" << std::endl;
1077 }
1078 for (size_t i = std::min(land1.land.size(), land2.land.size()); i != std::max(land1.land.size(), land2.land.size());
1079 ++i) {
1080 std::vector<std::pair<double, double> > lambda_n(land1.land[i]);
1081 for (size_t nr = 0; nr != land1.land[i].size(); ++nr) {
1082 lambda_n[nr] = std::make_pair(land1.land[i][nr].first, oper(land1.land[i][nr].second, 0));
1083 }
1084 // CHANGE
1085 // result.land[i] = lambda_n;
1086 result.land[i].swap(lambda_n);
1087 }
1088 }
1089 if (land2.land.size() > std::min(land1.land.size(), land2.land.size())) {
1090 if (operation_on_pair_of_landscapesDBG) {
1091 std::clog << "( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) ) " << std::endl;
1092 }
1093 for (size_t i = std::min(land1.land.size(), land2.land.size()); i != std::max(land1.land.size(), land2.land.size());
1094 ++i) {
1095 std::vector<std::pair<double, double> > lambda_n(land2.land[i]);
1096 for (size_t nr = 0; nr != land2.land[i].size(); ++nr) {
1097 lambda_n[nr] = std::make_pair(land2.land[i][nr].first, oper(0, land2.land[i][nr].second));
1098 }
1099 // CHANGE
1100 // result.land[i] = lambda_n;
1101 result.land[i].swap(lambda_n);
1102 }
1103 }
1104 if (operation_on_pair_of_landscapesDBG) {
1105 std::clog << "operation_on_pair_of_landscapes END\n";
1106 std::cin.ignore();
1107 }
1108 return result;
1109} // operation_on_pair_of_landscapes
1110
1111double compute_maximal_distance_non_symmetric(const Persistence_landscape& pl1, const Persistence_landscape& pl2) {
1112 bool dbg = false;
1113 if (dbg) std::clog << " compute_maximal_distance_non_symmetric \n";
1114 // this distance is not symmetric. It compute ONLY distance between inflection points of pl1 and pl2.
1115 double maxDist = 0;
1116 size_t minimalNumberOfLevels = std::min(pl1.land.size(), pl2.land.size());
1117 for (size_t level = 0; level != minimalNumberOfLevels; ++level) {
1118 if (dbg) {
1119 std::clog << "Level : " << level << std::endl;
1120 std::clog << "PL1 : \n";
1121 for (size_t i = 0; i != pl1.land[level].size(); ++i) {
1122 std::clog << "(" << pl1.land[level][i].first << "," << pl1.land[level][i].second << ") \n";
1123 }
1124 std::clog << "PL2 : \n";
1125 for (size_t i = 0; i != pl2.land[level].size(); ++i) {
1126 std::clog << "(" << pl2.land[level][i].first << "," << pl2.land[level][i].second << ") \n";
1127 }
1128 std::cin.ignore();
1129 }
1130
1131 int p2Count = 0;
1132 // In this case, I consider points at the infinity
1133 for (size_t i = 1; i != pl1.land[level].size() - 1; ++i) {
1134 while (true) {
1135 if ((pl1.land[level][i].first >= pl2.land[level][p2Count].first) &&
1136 (pl1.land[level][i].first <= pl2.land[level][p2Count + 1].first))
1137 break;
1138 p2Count++;
1139 }
1140 double val =
1141 fabs(function_value(pl2.land[level][p2Count], pl2.land[level][p2Count + 1], pl1.land[level][i].first) -
1142 pl1.land[level][i].second);
1143 if (maxDist <= val) maxDist = val;
1144
1145 if (dbg) {
1146 std::clog << pl1.land[level][i].first << "in [" << pl2.land[level][p2Count].first << ","
1147 << pl2.land[level][p2Count + 1].first << "] \n";
1148 std::clog << "pl1[level][i].second : " << pl1.land[level][i].second << std::endl;
1149 std::clog << "function_value( pl2[level][p2Count] , pl2[level][p2Count+1] , pl1[level][i].first ) : "
1150 << function_value(pl2.land[level][p2Count], pl2.land[level][p2Count + 1], pl1.land[level][i].first)
1151 << std::endl;
1152 std::clog << "val : " << val << std::endl;
1153 std::cin.ignore();
1154 }
1155 }
1156 }
1157
1158 if (dbg) std::clog << "minimalNumberOfLevels : " << minimalNumberOfLevels << std::endl;
1159
1160 if (minimalNumberOfLevels < pl1.land.size()) {
1161 for (size_t level = minimalNumberOfLevels; level != pl1.land.size(); ++level) {
1162 for (size_t i = 0; i != pl1.land[level].size(); ++i) {
1163 if (dbg) std::clog << "pl1[level][i].second : " << pl1.land[level][i].second << std::endl;
1164 if (maxDist < pl1.land[level][i].second) maxDist = pl1.land[level][i].second;
1165 }
1166 }
1167 }
1168 return maxDist;
1169}
1170
1171double compute_distance_of_landscapes(const Persistence_landscape& first, const Persistence_landscape& second,
1172 double p) {
1173 bool dbg = false;
1174 // This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a
1175 // time:
1176
1177 // first-second :
1178 Persistence_landscape lan = first - second;
1179
1180 //| first-second |:
1181 lan = lan.abs();
1182
1183 if (dbg) {
1184 std::clog << "Abs of difference ; " << lan << std::endl;
1185 getchar();
1186 }
1187
1188 if (p < std::numeric_limits<double>::max()) {
1189 // \int_{- \infty}^{+\infty}| first-second |^p
1190 double result;
1191 if (p != 1) {
1192 if (dbg) std::clog << "Power != 1, compute integral to the power p\n";
1193 result = lan.compute_integral_of_landscape(p);
1194 } else {
1195 if (dbg) std::clog << "Power = 1, compute integral \n";
1196 result = lan.compute_integral_of_landscape();
1197 }
1198 // (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p)
1199 return pow(result, 1.0 / p);
1200 } else {
1201 // p == infty
1202 if (dbg) std::clog << "Power = infty, compute maximum \n";
1203 return lan.compute_maximum();
1204 }
1205}
1206
1207double compute_max_norm_distance_of_landscapes(const Persistence_landscape& first,
1208 const Persistence_landscape& second) {
1209 return std::max(compute_maximal_distance_non_symmetric(first, second),
1210 compute_maximal_distance_non_symmetric(second, first));
1211}
1212
1213bool comparePairsForMerging(std::pair<double, unsigned> first, std::pair<double, unsigned> second) {
1214 return (first.first < second.first);
1215}
1216
1217double compute_inner_product(const Persistence_landscape& l1, const Persistence_landscape& l2) {
1218 bool dbg = false;
1219 double result = 0;
1220
1221 for (size_t level = 0; level != std::min(l1.size(), l2.size()); ++level) {
1222 if (dbg) {
1223 std::clog << "Computing inner product for a level : " << level << std::endl;
1224 getchar();
1225 }
1226 auto&& l1_land_level = l1.land[level];
1227 auto&& l2_land_level = l2.land[level];
1228
1229 if (l1_land_level.size() * l2_land_level.size() == 0) continue;
1230
1231 // endpoints of the interval on which we will compute the inner product of two locally linear functions:
1232 double x1 = -std::numeric_limits<int>::max();
1233 double x2;
1234 if (l1_land_level[1].first < l2_land_level[1].first) {
1235 x2 = l1_land_level[1].first;
1236 } else {
1237 x2 = l2_land_level[1].first;
1238 }
1239
1240 // iterators for the landscapes l1 and l2
1241 size_t l1It = 0;
1242 size_t l2It = 0;
1243
1244 while ((l1It < l1_land_level.size() - 1) && (l2It < l2_land_level.size() - 1)) {
1245 // compute the value of a inner product on a interval [x1,x2]
1246
1247 double a, b, c, d;
1248
1249 if (l1_land_level[l1It + 1].first != l1_land_level[l1It].first) {
1250 a = (l1_land_level[l1It + 1].second - l1_land_level[l1It].second) /
1251 (l1_land_level[l1It + 1].first - l1_land_level[l1It].first);
1252 } else {
1253 a = 0;
1254 }
1255 b = l1_land_level[l1It].second - a * l1_land_level[l1It].first;
1256 if (l2_land_level[l2It + 1].first != l2_land_level[l2It].first) {
1257 c = (l2_land_level[l2It + 1].second - l2_land_level[l2It].second) /
1258 (l2_land_level[l2It + 1].first - l2_land_level[l2It].first);
1259 } else {
1260 c = 0;
1261 }
1262 d = l2_land_level[l2It].second - c * l2_land_level[l2It].first;
1263
1264 double contributionFromThisPart = (a * c * x2 * x2 * x2 / 3 + (a * d + b * c) * x2 * x2 / 2 + b * d * x2) -
1265 (a * c * x1 * x1 * x1 / 3 + (a * d + b * c) * x1 * x1 / 2 + b * d * x1);
1266
1267 result += contributionFromThisPart;
1268
1269 if (dbg) {
1270 std::clog << "[l1_land_level[l1It].first,l1_land_level[l1It+1].first] : " << l1_land_level[l1It].first
1271 << " , " << l1_land_level[l1It + 1].first << std::endl;
1272 std::clog << "[l2_land_level[l2It].first,l2_land_level[l2It+1].first] : " << l2_land_level[l2It].first
1273 << " , " << l2_land_level[l2It + 1].first << std::endl;
1274 std::clog << "a : " << a << ", b : " << b << " , c: " << c << ", d : " << d << std::endl;
1275 std::clog << "x1 : " << x1 << " , x2 : " << x2 << std::endl;
1276 std::clog << "contributionFromThisPart : " << contributionFromThisPart << std::endl;
1277 std::clog << "result : " << result << std::endl;
1278 getchar();
1279 }
1280
1281 // we have two intervals in which functions are constant:
1282 // [l1_land_level[l1It].first , l1_land_level[l1It+1].first]
1283 // and
1284 // [l2_land_level[l2It].first , l2_land_level[l2It+1].first]
1285 // We also have an interval [x1,x2]. Since the intervals in the landscapes cover the whole R, then it is clear
1286 // that x2
1287 // is either l1_land_level[l1It+1].first of l2_land_level[l2It+1].first or both. Lets test it.
1288 if (x2 == l1_land_level[l1It + 1].first) {
1289 if (x2 == l2_land_level[l2It + 1].first) {
1290 // in this case, we increment both:
1291 ++l2It;
1292 if (dbg) {
1293 std::clog << "Incrementing both \n";
1294 }
1295 } else {
1296 if (dbg) {
1297 std::clog << "Incrementing first \n";
1298 }
1299 }
1300 ++l1It;
1301 } else {
1302 // in this case we increment l2It
1303 ++l2It;
1304 if (dbg) {
1305 std::clog << "Incrementing second \n";
1306 }
1307 }
1308
1309 if ( l1It + 1 >= l1_land_level.size() )break;
1310 if ( l2It + 1 >= l2_land_level.size() )break;
1311
1312 // Now, we shift x1 and x2:
1313 x1 = x2;
1314 if (l1_land_level[l1It + 1].first < l2_land_level[l2It + 1].first) {
1315 x2 = l1_land_level[l1It + 1].first;
1316 } else {
1317 x2 = l2_land_level[l2It + 1].first;
1318 }
1319 }
1320 }
1321 return result;
1322}
1323
1324void Persistence_landscape::plot(const char* filename, double xRangeBegin, double xRangeEnd, double yRangeBegin,
1325 double yRangeEnd, int from, int to) {
1326 // this program create a gnuplot script file that allows to plot persistence diagram.
1327 std::ofstream out;
1328
1329 std::ostringstream gnuplot_script;
1330 gnuplot_script << filename << "_GnuplotScript";
1331 out.open(gnuplot_script.str().c_str());
1332
1333 if ((xRangeBegin != std::numeric_limits<double>::max()) || (xRangeEnd != std::numeric_limits<double>::max()) ||
1334 (yRangeBegin != std::numeric_limits<double>::max()) || (yRangeEnd != std::numeric_limits<double>::max())) {
1335 out << "set xrange [" << xRangeBegin << " : " << xRangeEnd << "]" << std::endl;
1336 out << "set yrange [" << yRangeBegin << " : " << yRangeEnd << "]" << std::endl;
1337 }
1338
1339 if (from == std::numeric_limits<int>::max()) {
1340 from = 0;
1341 }
1342 if (to == std::numeric_limits<int>::max()) {
1343 to = this->land.size();
1344 }
1345
1346 out << "plot ";
1347 for (size_t lambda = std::min((size_t)from, this->land.size()); lambda != std::min((size_t)to, this->land.size());
1348 ++lambda) {
1349 // out << " '-' using 1:2 title 'l" << lambda << "' with lp";
1350 out << " '-' using 1:2 notitle with lp";
1351 if (lambda + 1 != std::min((size_t)to, this->land.size())) {
1352 out << ", \\";
1353 }
1354 out << std::endl;
1355 }
1356
1357 for (size_t lambda = std::min((size_t)from, this->land.size()); lambda != std::min((size_t)to, this->land.size());
1358 ++lambda) {
1359 for (size_t i = 1; i != this->land[lambda].size() - 1; ++i) {
1360 out << this->land[lambda][i].first << " " << this->land[lambda][i].second << std::endl;
1361 }
1362 out << "EOF" << std::endl;
1363 }
1364 std::clog << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'"
1365 << gnuplot_script.str().c_str() << "\'\"" << std::endl;
1366}
1367
1368} // namespace Persistence_representations
1369} // namespace Gudhi
1370
1371#endif // PERSISTENCE_LANDSCAPE_H_
A class implementing persistence landscapes data structures.
Definition: Persistence_landscape.h:60
friend double compute_inner_product(const Persistence_landscape &l1, const Persistence_landscape &l2)
Definition: Persistence_landscape.h:1217
bool operator==(const Persistence_landscape &rhs) const
Definition: Persistence_landscape.h:472
friend double compute_distance_of_landscapes(const Persistence_landscape &first, const Persistence_landscape &second, double p)
Definition: Persistence_landscape.h:1171
friend Persistence_landscape operator+(const Persistence_landscape &first, const Persistence_landscape &second)
Definition: Persistence_landscape.h:145
size_t number_of_projections_to_R() const
Definition: Persistence_landscape.h:314
double compute_scalar_product(const Persistence_landscape &second) const
Definition: Persistence_landscape.h:413
friend Persistence_landscape operator-(const Persistence_landscape &first, const Persistence_landscape &second)
Definition: Persistence_landscape.h:152
friend double compute_max_norm_distance_of_landscapes(const Persistence_landscape &first, const Persistence_landscape &second)
Definition: Persistence_landscape.h:1207
size_t number_of_vectorize_functions() const
Definition: Persistence_landscape.h:336
Persistence_landscape abs()
Definition: Persistence_landscape.h:809
friend Persistence_landscape add_two_landscapes(const Persistence_landscape &land1, const Persistence_landscape &land2)
Definition: Persistence_landscape.h:129
double find_max(unsigned lambda) const
Definition: Persistence_landscape.h:651
double compute_integral_of_a_level_of_a_landscape(size_t level) const
Definition: Persistence_landscape.h:672
friend Persistence_landscape subtract_two_landscapes(const Persistence_landscape &land1, const Persistence_landscape &land2)
Definition: Persistence_landscape.h:137
double compute_integral_of_landscape() const
Definition: Persistence_landscape.h:660
Persistence_landscape()
Definition: Persistence_landscape.h:65
Persistence_landscape operator-=(const Persistence_landscape &rhs)
Definition: Persistence_landscape.h:181
double distance(const Persistence_landscape &second, double power=1) const
Definition: Persistence_landscape.h:400
Persistence_landscape operator+=(const Persistence_landscape &rhs)
Definition: Persistence_landscape.h:173
double compute_value_at_a_given_point(unsigned level, double x) const
Definition: Persistence_landscape.h:722
double project_to_R(int number_of_function) const
Definition: Persistence_landscape.h:306
Persistence_landscape operator/=(double x)
Definition: Persistence_landscape.h:198
friend Persistence_landscape operator*(const Persistence_landscape &first, double con)
Definition: Persistence_landscape.h:159
std::vector< double > vectorize(int number_of_function) const
Definition: Persistence_landscape.h:320
double compute_maximum() const
Definition: Persistence_landscape.h:217
std::pair< double, double > get_y_range(size_t level=0) const
Definition: Persistence_landscape.h:422
void compute_average(const std::vector< Persistence_landscape * > &to_average)
Definition: Persistence_landscape.h:342
Persistence_landscape operator*=(double x)
Definition: Persistence_landscape.h:190
void print_to_file(const char *filename) const
Definition: Persistence_landscape.h:907
double operator()(unsigned level, double x) const
Definition: Persistence_landscape.h:257
friend Persistence_landscape operator*(double con, const Persistence_landscape &first)
Definition: Persistence_landscape.h:166
friend std::ostream & operator<<(std::ostream &out, Persistence_landscape &land)
Definition: Persistence_landscape.h:781
size_t size() const
Definition: Persistence_landscape.h:284
bool operator!=(const Persistence_landscape &rhs) const
Definition: Persistence_landscape.h:212
void load_landscape_from_file(const char *filename)
Definition: Persistence_landscape.h:919