Persistence_landscape_on_grid.h

/*    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):
 *      - 2025/06 Hannah Schreiber: Various small bug fixes (missing `inline`s, `DEBUG_TRACES`s etc.)
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef PERSISTENCE_LANDSCAPE_ON_GRID_H_
#define PERSISTENCE_LANDSCAPE_ON_GRID_H_
 
// standard include
#ifdef DEBUG_TRACES
#include <iostream>   // std::cerr, std::clog
#endif
#include <cstddef>    // std::size_t
#include <cstdint>    // std::uint16_t
#include <ostream>    // std::ostream
#include <fstream>    // std::ofstream, std::ifstream
#include <sstream>    // std::istringstream, std::ostringstream
#include <stdexcept>  // std::invalid_argument
#include <limits>     // std::numeric_limits
#include <algorithm>  // std::make_heap, std::pop_heap, std::push_heap
#include <cmath>      // std::pow
#include <utility>    // std::pair
#include <vector>
#include <string>
 
// gudhi include
#include <gudhi/read_persistence_from_file.h>
#include <gudhi/common_persistence_representations.h>
#include <gudhi/Debug_utils.h>
 
namespace Gudhi {
namespace Persistence_representations {
 
// pre declaration needed before C++20 for friends with templates defined inside a class
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
{
 public:
  Persistence_landscape_on_grid()
      : grid_min_(0),
        grid_max_(0),
        number_of_functions_for_vectorization_(0),
        number_of_functions_for_projections_to_reals_(0)
  {}

  Persistence_landscape_on_grid(const std::vector<std::pair<double, double> >& p,
                                double grid_min,
                                double grid_max,
                                std::size_t number_of_points)
  {
    this->_set_up_values_of_landscapes(p, grid_min, grid_max, number_of_points);
  }

  Persistence_landscape_on_grid(const std::vector<std::pair<double, double> >& p,
                                double grid_min,
                                double grid_max,
                                std::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(const char* filename,
                                double grid_min,
                                double grid_max,
                                std::size_t number_of_points,
                                unsigned number_of_levels_of_landscape,
                                std::uint16_t dimension = std::numeric_limits<std::uint16_t>::max())
  {
    std::vector<std::pair<double, double> > p;
    if (dimension == std::numeric_limits<std::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(const char* filename,
                                double grid_min,
                                double grid_max,
                                std::size_t number_of_points,
                                std::uint16_t dimension = std::numeric_limits<std::uint16_t>::max())
  {
    std::vector<std::pair<double, double> > p;
    if (dimension == std::numeric_limits<std::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(const char* filename,
                                std::size_t number_of_points,
                                unsigned number_of_levels_of_landscape,
                                std::uint16_t dimension = std::numeric_limits<std::uint16_t>::max())
  {
    std::vector<std::pair<double, double> > p;
    if (dimension == std::numeric_limits<std::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 (std::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);
  }

  Persistence_landscape_on_grid(const char* filename,
                                std::size_t number_of_points,
                                std::uint16_t dimension = std::numeric_limits<std::uint16_t>::max())
  {
    std::vector<std::pair<double, double> > p;
    if (dimension == std::numeric_limits<std::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 (std::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);
  }

  void load_landscape_from_file(const char* filename);

  void print_to_file(const char* filename) const;

  double compute_integral_of_landscape() const
  {
    std::size_t maximal_level = this->number_of_nonzero_levels();
    double result = 0;
    for (std::size_t i = 0; i != maximal_level; ++i) {
      result += this->compute_integral_of_landscape(i);
    }
    return result;
  }

  double compute_integral_of_landscape(std::size_t level) const
  {
    double result = 0;
    double dx = (this->grid_max_ - this->grid_min_) / static_cast<double>(this->values_of_landscapes_.size() - 1);
 
#ifdef DEBUG_TRACES
    std::clog << "this->grid_max : " << this->grid_max_ << std::endl;
    std::clog << "this->grid_min : " << this->grid_min_ << std::endl;
    std::clog << "this->values_of_landscapes.size() : " << this->values_of_landscapes_.size() << std::endl;
#endif
 
    double previous_x = this->grid_min_ - dx;
    double previous_y = 0;
    for (std::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];
 
#ifdef DEBUG_TRACES
      std::clog << "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::clog << "this->values_of_landscapes[i][level] : " << this->values_of_landscapes_[i][level] << std::endl;
      std::clog << "previous_y : " << previous_y << std::endl;
      std::clog << "current_y : " << current_y << std::endl;
      std::clog << "dx : " << dx << std::endl;
      std::clog << "0.5*dx*( previous_y + current_y ); " << 0.5 * dx * (previous_y + current_y) << std::endl;
#endif
 
      result += 0.5 * dx * (previous_y + current_y);
      previous_x = current_x;
      previous_y = current_y;
    }
    return result;
  }

  double compute_integral_of_landscape(double p) const
  {
    std::size_t maximal_level = this->number_of_nonzero_levels();
    double result = 0;
    for (std::size_t i = 0; i != maximal_level; ++i) {
      result += this->compute_integral_of_landscape(p, i);
    }
    return result;
  }

  double compute_integral_of_landscape(double p, std::size_t level) const
  {
    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];
 
#ifdef DEBUG_TRACES
    std::clog << "dx : " << dx << std::endl;
    std::clog << "previous_x : " << previous_x << std::endl;
    std::clog << "previous_y : " << previous_y << std::endl;
    std::clog << "power : " << p << std::endl;
#endif
 
    for (std::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];
 
#ifdef DEBUG_TRACES
      std::clog << "current_y : " << current_y << std::endl;
#endif
 
      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;
 
#ifdef DEBUG_TRACES
      std::clog << "A line passing through points : (" << previous_x << "," << previous_y << ") and (" << current_x
                << "," << current_y << ") is : " << a << "x+" << b << std::endl;
#endif
 
      // 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)) * (std::pow((a * current_x + b), p + 1) - std::pow((a * previous_x + b), p + 1));
      } else {
        value_to_add = (current_x - previous_x) * (std::pow(b, p));
      }
      result += value_to_add;
#ifdef DEBUG_TRACES
      std::clog << "Increasing result by : " << value_to_add << std::endl;
      std::clog << "result : " << result << std::endl;
#endif
      previous_x = current_x;
      previous_y = current_y;
    }
#ifdef DEBUG_TRACES
    std::clog << "The total result is : " << result << std::endl;
#endif
    return result;
  }

  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 (std::size_t i = 0; i != land.values_of_landscapes_.size(); ++i) {
      out << x << " : ";
      for (std::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;
  }

  double compute_value_at_a_given_point(unsigned level, double x) const
  {
    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);
    std::size_t position = std::size_t((x - this->grid_min_) / dx);
 
#ifdef DEBUG_TRACES
    std::clog << "This is a procedure compute_value_at_a_given_point \n";
    std::clog << "level : " << level << std::endl;
    std::clog << "x : " << x << std::endl;
    std::clog << "position : " << position << std::endl;
#endif
 
    // 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;
  }
 
  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;
  }

  template <typename Operation>
  friend Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid(
      const Persistence_landscape_on_grid& land1,
      const Persistence_landscape_on_grid& land2)
  {
    // TODO: replace this with a GUDHI_CHECK?
    // first we need to check if the domains are the same:
    if (!check_if_defined_on_the_same_domain(land1, land2)) throw std::invalid_argument("Two grids are not compatible");
 
    Operation 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 (std::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 (std::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;
  }
 
  // TODO: `add_two_landscapes` and `subtract_two_landscapes` do not seem to be used outside of resp.
  // `operator+` and `operator-` and their doc is marked private. Are they really necessary?
  // Can `operation_on_pair_of_landscapes_on_grid` not be directly called in `operator+` and `operator-`?

  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);
  }

  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);
  }

  friend Persistence_landscape_on_grid operator+(const Persistence_landscape_on_grid& first,
                                                 const Persistence_landscape_on_grid& second)
  {
    return add_two_landscapes(first, second);
  }

  friend Persistence_landscape_on_grid operator-(const Persistence_landscape_on_grid& first,
                                                 const Persistence_landscape_on_grid& second)
  {
    return subtract_two_landscapes(first, second);
  }

  friend Persistence_landscape_on_grid operator*(const Persistence_landscape_on_grid& first, double con)
  {
    return first._multiply_landscape_by_real_number_not_overwrite(con);
  }

  friend Persistence_landscape_on_grid operator*(double con, const Persistence_landscape_on_grid& first)
  {
    return first._multiply_landscape_by_real_number_not_overwrite(con);
  }

  Persistence_landscape_on_grid operator+=(const Persistence_landscape_on_grid& rhs)
  {
    *this = *this + rhs;
    return *this;
  }

  Persistence_landscape_on_grid operator-=(const Persistence_landscape_on_grid& rhs)
  {
    *this = *this - rhs;
    return *this;
  }

  Persistence_landscape_on_grid operator*=(double x)
  {
    *this = *this * x;
    return *this;
  }

  Persistence_landscape_on_grid operator/=(double x)
  {
    if (x == 0) throw("In operator /=, division by 0. Program terminated.");
    *this = *this * (1 / x);
    return *this;
  }

  bool operator==(const Persistence_landscape_on_grid& rhs) const
  {
    if (this->values_of_landscapes_.size() != rhs.values_of_landscapes_.size()) {
#ifdef DEBUG_TRACES
      std::clog << "values_of_landscapes of incompatible sizes\n";
#endif
      return false;
    }
    if (!almost_equal(this->grid_min_, rhs.grid_min_)) {
#ifdef DEBUG_TRACES
      std::clog << "grid_min not equal\n";
#endif
      return false;
    }
    if (!almost_equal(this->grid_max_, rhs.grid_max_)) {
#ifdef DEBUG_TRACES
      std::clog << "grid_max not equal\n";
#endif
      return false;
    }
    for (std::size_t i = 0; i != this->values_of_landscapes_.size(); ++i) {
      for (std::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])) {
#ifdef DEBUG_TRACES
          std::clog << "Problem in the position : " << i << " of values_of_landscapes. \n";
          std::clog << this->values_of_landscapes_[i][aa] << " " << rhs.values_of_landscapes_[i][aa] << std::endl;
#endif
          return false;
        }
      }
    }
    return true;
  }

  bool operator!=(const Persistence_landscape_on_grid& rhs) const { return !((*this) == rhs); }

  double compute_maximum() const
  {
    // since the function can only be entirely positive or negative, the maximal value will be an extrema value in the
    // arrays:
    double max_value = -std::numeric_limits<double>::max();
    for (std::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;
  }

  std::pair<double, double> compute_minimum_maximum() const
  {
    // since the function can only be entirely positive or negative, the maximal value will be an extrema value in the
    // arrays:
    double max_value = -std::numeric_limits<double>::max();
    double min_value = 0;
    for (std::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);
  }

  std::pair<double, double> get_x_range(std::size_t level = 0) const
  {
    return std::make_pair(this->grid_min_, this->grid_max_);
  }

  std::pair<double, double> get_y_range(std::size_t level = 0) const { return this->compute_minimum_maximum(); }

  std::size_t number_of_nonzero_levels() const
  {
    std::size_t result = 0;
    for (std::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;
  }

  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);
    }
  }

  double operator()(unsigned int level, double x) const { return this->compute_value_at_a_given_point(level, x); }

  friend 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 (std::size_t i = 0; i != first.values_of_landscapes_.size(); ++i) {
      for (std::size_t j = 0;
           j != std::min(first.values_of_landscapes_[i].size(), second.values_of_landscapes_[i].size());
           ++j) {
        if (result < std::abs(first.values_of_landscapes_[i][j] - second.values_of_landscapes_[i][j])) {
          result = std::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 (std::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 (std::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;
  }

  void abs()
  {
    for (std::size_t i = 0; i != this->values_of_landscapes_.size(); ++i) {
      for (std::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]);
      }
    }
  }

  std::size_t size() const { return this->number_of_nonzero_levels(); }

  double find_max(unsigned lambda) const
  {
    double max_value = -std::numeric_limits<double>::max();
    for (std::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;
  }

  friend double compute_inner_product(const Persistence_landscape_on_grid& l1, const Persistence_landscape_on_grid& l2)
  {
    // TODO: replace with a GUDHI_CHECK?
    if (!check_if_defined_on_the_same_domain(l1, l2))
      throw std::invalid_argument("Landscapes are not defined on the same grid, the program will now terminate");
    std::size_t maximal_level = l1.number_of_nonzero_levels();
    double result = 0;
    for (std::size_t i = 0; i != maximal_level; ++i) {
      result += compute_inner_product(l1, l2, i);
    }
    return result;
  }

  friend double compute_inner_product(const Persistence_landscape_on_grid& l1,
                                      const Persistence_landscape_on_grid& l2,
                                      std::size_t level)
  {
    // TODO: replace with a GUDHI_CHECK?
    if (!check_if_defined_on_the_same_domain(l1, l2))
      throw std::invalid_argument("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 (std::size_t i = 0; i != l1.values_of_landscapes_.size(); ++i) {
#ifdef DEBUG_TRACES
      std::clog << "i : " << i << std::endl;
#endif
 
      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];
 
#ifdef DEBUG_TRACES
      std::clog << "previous_x  : " << previous_x << std::endl;
      std::clog << "previous_y_l1 : " << previous_y_l1 << std::endl;
      std::clog << "current_y_l1 : " << current_y_l1 << std::endl;
      std::clog << "previous_y_l2 : " << previous_y_l2 << std::endl;
      std::clog << "current_y_l2 : " << current_y_l2 << std::endl;
#endif
 
      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;
 
#ifdef DEBUG_TRACES
      std::clog << "Here are the formulas for a line: \n";
      std::clog << "a : " << a << std::endl;
      std::clog << "b : " << b << std::endl;
      std::clog << "c : " << c << std::endl;
      std::clog << "d : " << d << std::endl;
#endif
 
      // 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);
 
#ifdef DEBUG_TRACES
      std::clog << "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::clog << "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;
#endif
 
      result += added_value;
 
#ifdef DEBUG_TRACES
      std::clog << "added_value : " << added_value << std::endl;
      std::clog << "result : " << result << std::endl;
#endif
 
      previous_x = current_x;
      previous_y_l1 = current_y_l1;
      previous_y_l2 = current_y_l2;
    }
    return result;
  }

  friend double compute_distance_of_landscapes_on_grid(const Persistence_landscape_on_grid& first,
                                                       const Persistence_landscape_on_grid& second,
                                                       double p)
  {
    // This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p).
    // We will do it one step at a time:
 
#ifdef DEBUG_TRACES
    std::clog << "first : " << first << std::endl;
    std::clog << "second : " << second << std::endl;
#endif
 
    // first-second :
    Persistence_landscape_on_grid lan = first - second;
 
#ifdef DEBUG_TRACES
    std::clog << "Difference : " << lan << std::endl;
#endif
 
    //| first-second |:
    lan.abs();
 
#ifdef DEBUG_TRACES
    std::clog << "Abs : " << lan << std::endl;
#endif
 
    if (p < std::numeric_limits<double>::max()) {
      // \int_{- \infty}^{+\infty}| first-second |^p
      double result;
      if (p != 1) {
#ifdef DEBUG_TRACES
        std::clog << "p : " << p << std::endl;
#endif
        result = lan.compute_integral_of_landscape(p);
#ifdef DEBUG_TRACES
        std::clog << "integral : " << result << std::endl;
#endif
      } else {
        result = lan.compute_integral_of_landscape();
#ifdef DEBUG_TRACES
        std::clog << "integral, without power : " << result << std::endl;
#endif
      }
      // (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p)
      return std::pow(result, 1.0 / p);
    } else {
      // p == infty
      return lan.compute_maximum();
    }
  }
 
  // Functions that are needed for that class to implement the concept.

  double project_to_R(int number_of_function) const
  {
    return this->compute_integral_of_landscape(static_cast<std::size_t>(number_of_function));
  }

  std::size_t number_of_projections_to_R() const { return number_of_functions_for_projections_to_reals_; }

  std::vector<double> vectorize(int number_of_function) const
  {
    // TODO: GUDHI_CHECK instead?
    // TODO(PD) think of something smarter over here
    if (number_of_function < 0 || static_cast<std::size_t>(number_of_function) >= this->values_of_landscapes_.size()) {
      throw std::invalid_argument("Wrong number of function.");
    }
    std::vector<double> v(this->values_of_landscapes_.size());
    for (std::size_t i = 0; i != this->values_of_landscapes_.size(); ++i) {
      v[i] = 0;
      if (this->values_of_landscapes_[i].size() > static_cast<std::size_t>(number_of_function)) {
        v[i] = this->values_of_landscapes_[i][number_of_function];
      }
    }
    return v;
  }

  std::size_t number_of_vectorize_functions() const { return number_of_functions_for_vectorization_; }

  void compute_average(const std::vector<Persistence_landscape_on_grid*>& to_average)
  {
    // 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 (std::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_;
 
#ifdef DEBUG_TRACES
    std::clog << "Computations of average. The data from the current landscape have been cleared. We are ready to do "
                 "the computations. \n";
#endif
 
    // for every point in the grid:
    for (std::size_t grid_point = 0; grid_point != (to_average[0])->values_of_landscapes_.size(); ++grid_point) {
      // set up a vector of the correct size:
      std::size_t maximal_size_of_vector = 0;
      for (std::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);
 
#ifdef DEBUG_TRACES
      std::clog << "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";
#endif
 
      // and compute an arithmetic average:
      for (std::size_t land_no = 0; land_no != to_average.size(); ++land_no) {
        // summing:
        for (std::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 (std::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

  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);
    }
  }

  double compute_scalar_product(const Persistence_landscape_on_grid& second)
  {
    return compute_inner_product((*this), second);
  }
 
  // end of implementation of functions needed for concepts.

  std::vector<std::vector<double> > output_for_visualization() const { return this->values_of_landscapes_; }

  void plot(const char* filename, std::size_t from, std::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);
  }

  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(),
            std::size_t from = std::numeric_limits<std::size_t>::max(),
            std::size_t to = std::numeric_limits<std::size_t>::max()) const;
 
 private:
  double grid_min_;
  double grid_max_;
  std::vector<std::vector<double> > values_of_landscapes_;
  std::size_t number_of_functions_for_vectorization_;
  std::size_t number_of_functions_for_projections_to_reals_;
 
  void _set_up_values_of_landscapes(const std::vector<std::pair<double, double> >& p,
                                    double grid_min,
                                    double grid_max,
                                    std::size_t number_of_points,
                                    unsigned number_of_levels = std::numeric_limits<unsigned>::max());
  Persistence_landscape_on_grid _multiply_landscape_by_real_number_not_overwrite(double x) const;
};
 
inline void Persistence_landscape_on_grid::_set_up_values_of_landscapes(
    const std::vector<std::pair<double, double> >& p,
    double grid_min,
    double grid_max,
    std::size_t number_of_points,
    unsigned number_of_levels)
{
#ifdef DEBUG_TRACES
  std::clog << "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::clog << "Here are the intervals at our disposal : \n";
  for (std::size_t i = 0; i != p.size(); ++i) {
    std::clog << p[i].first << " , " << p[i].second << std::endl;
  }
#endif
 
  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 (std::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<std::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<std::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 std::invalid_argument(
        "Wrong parameters of grid_min and grid_max given to the procedure. The program will now terminate.");
  }
 
  double dx = (grid_max - grid_min) / static_cast<double>(number_of_points);
  // for every interval in the diagram:
  for (std::size_t int_no = 0; int_no != p.size(); ++int_no) {
    std::size_t grid_interval_begin = (p[int_no].first - grid_min) / dx;
    std::size_t grid_interval_end = (p[int_no].second - grid_min) / dx;
    std::size_t grid_interval_midpoint = (std::size_t)(0.5 * (grid_interval_begin + grid_interval_end));
 
#ifdef DEBUG_TRACES
    std::clog << "Considering an interval : " << p[int_no].first << "," << p[int_no].second << std::endl;
 
    std::clog << "grid_interval_begin : " << grid_interval_begin << std::endl;
    std::clog << "grid_interval_end : " << grid_interval_end << std::endl;
    std::clog << "grid_interval_midpoint : " << grid_interval_midpoint << std::endl;
#endif
 
    double landscape_value = dx;
    for (std::size_t i = grid_interval_begin + 1; i < grid_interval_midpoint; ++i) {
#ifdef DEBUG_TRACES
      std::clog << "Adding landscape value (going up) for a point : " << i << " equal : " << landscape_value
                << std::endl;
#endif
      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 (std::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);
        }
 
#ifdef DEBUG_TRACES
        std::clog << "Adding landscape value (going down) for a point : " << i << " equal : " << landscape_value
                  << std::endl;
#endif
      }
      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 (std::size_t pt = 0; pt != this->values_of_landscapes_.size(); ++pt) {
      for (std::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 (std::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
 
inline 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()) {
#ifdef DEBUG_TRACES
    std::cerr << "The file : " << filename << " do not exist. The program will now terminate \n";
#endif
    throw std::invalid_argument("The persistence landscape file do not exist.");
  }
 
  std::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 (std::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::getline(in, line);
    std::istringstream stream(line);
    while (stream >> number) {
      v[i].push_back(number);
    }
  }
  this->values_of_landscapes_ = v;
  in.close();
}
 
inline 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 (std::size_t i = 0; i != this->values_of_landscapes_.size(); ++i) {
    for (std::size_t j = 0; j != this->values_of_landscapes_[i].size(); ++j) {
      out << this->values_of_landscapes_[i][j] << " ";
    }
    out << std::endl;
  }
 
  out.close();
}
 
inline void Persistence_landscape_on_grid::plot(const char* filename,
                                                double min_x,
                                                double max_x,
                                                double min_y,
                                                double max_y,
                                                std::size_t from_,
                                                std::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;
  }
 
  std::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);
 
  std::size_t from = 0;
  if (from_ != std::numeric_limits<std::size_t>::max()) {
    if (from_ < number_of_nonzero_levels) {
      from = from_;
    } else {
      return;
    }
  }
  std::size_t to = number_of_nonzero_levels;
  if (to_ != std::numeric_limits<std::size_t>::max()) {
    if (to_ < number_of_nonzero_levels) {
      to = to_;
    }
  }
 
  out << "plot ";
  for (std::size_t lambda = from; lambda != to; ++lambda) {
    out << "     '-' using 1:2 notitle with lp";
    if (lambda + 1 != to) {
      out << ", \\";
    }
    out << std::endl;
  }
 
  for (std::size_t lambda = from; lambda != to; ++lambda) {
    double point = this->grid_min_;
    for (std::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;
  }
#ifdef DEBUG_TRACES
  std::clog << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'"
            << gnuplot_script.str().c_str() << "\'\"" << std::endl;
#endif
}
 
inline Persistence_landscape_on_grid Persistence_landscape_on_grid::_multiply_landscape_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 (std::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 (std::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;
}
 
}  // namespace Persistence_representations
}  // namespace Gudhi
 
#endif  // PERSISTENCE_LANDSCAPE_ON_GRID_H_