a few more missing files

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/gudhi_stat@1617 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 4adab1d225550cfb24f65fba48d26af24860c83e

3 files changed, 2899 insertions, 0 deletions

diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h
new file mode 100644
index 00000000..82ee5878
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape.h
@@ -0,0 +1,1508 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+#ifndef Persistence_landscapes_H
+#define Persistence_landscapes_H
+
+//standard include
+#include <cmath>
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <fstream>
+#include <sstream>
+#include <algorithm>
+#include <unistd.h>
+
+
+//gudhi include
+#include <gudhi/abstract_classes/Abs_Vectorized_topological_data.h>
+#include <gudhi/abstract_classes/Abs_Topological_data_with_averages.h>
+#include <gudhi/abstract_classes/Abs_Topological_data_with_distances.h>
+#include <gudhi/abstract_classes/Abs_Real_valued_topological_data.h>
+#include <gudhi/abstract_classes/Abs_Topological_data_with_scalar_product.h>
+#include <gudhi/concretizations/read_persitence_from_file.h>
+#include <gudhi/common.h>
+using namespace std;
+
+
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+
+
+
+
+
+/**
+ * Extra functions needed in construction of barcodes.
+**/
+double minus_length( std::pair<double,double> a )
+{
+    return a.first-a.second;
+}
+double birth_plus_deaths( std::pair<double,double> a )
+{
+    return a.first+a.second;
+}
+
+
+
+/**
+ * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two points.
+**/
+std::pair<double,double> compute_parameters_of_a_line( std::pair<double,double> p1 , std::pair<double,double> p2 )
+{
+    double a = (p2.second-p1.second)/( p2.first - p1.first );
+    double b = p1.second - a*p1.first;
+    return std::make_pair(a,b);
+}
+
+/**
+ * This procedure given two points which lies on the opposide sides of x axis, compute x for which the line connecting those two points crosses x axis.
+**/
+double find_zero_of_a_line_segment_between_those_two_points ( std::pair<double,double> p1, std::pair<double,double> p2 )
+{
+    if ( p1.first == p2.first )return p1.first;
+    if ( p1.second*p2.second > 0 )
+    {
+        std::ostringstream errMessage;
+        errMessage <<"In function find_zero_of_a_line_segment_between_those_two_points the agguments are: (" << p1.first << "," << p1.second << ") and (" << p2.first << "," << p2.second << "). There is no zero in line between those two points. Program terminated.";
+        std::string errMessageStr = errMessage.str();
+        const char* err = errMessageStr.c_str();
+        throw(err);
+    }
+    //we assume here, that x \in [ p1.first, p2.first ] and p1 and p2 are points between which we will put the line segment
+    double a = (p2.second - p1.second)/(p2.first - p1.first);
+    double b = p1.second - a*p1.first;
+    //cerr << "Line crossing points : (" << p1.first << "," << p1.second << ") oraz (" << p2.first << "," << p2.second << ") : \n";
+    //cerr << "a : " << a << " , b : " << b << " , x : " << x << endl;
+    return -b/a;
+}
+
+
+
+
+/**
+ * Lexicographical ordering of points	.
+**/
+bool compare_points_sorting( std::pair<double,double> f, std::pair<double,double> s )
+{
+    if ( f.first < s.first )
+    {
+        return true;
+    }
+    else
+    {//f.first >= s.first
+        if ( f.first > s.first )
+        {
+            return false;
+        }
+        else
+        {//f.first == s.first
+            if ( f.second > s.second )
+            {
+                return true;
+            }
+            else
+            {
+                return false;
+            }
+        }
+    }
+}
+
+
+
+/**
+ * This procedure takes two points in R^2 and a double value x. It conputes the line pasing through those two points and return the value of that linear function at x.
+**/
+double function_value ( std::pair<double,double> p1, std::pair<double,double> p2 , double x )
+{
+    //we assume here, that x \in [ p1.first, p2.first ] and p1 and p2 are points between which we will put the line segment
+    double a = (p2.second - p1.second)/(p2.first - p1.first);
+    double b = p1.second - a*p1.first;
+    return (a*x+b);
+}
+
+
+
+
+
+/**
+ * A clas implementing persistence landascpes data structures. For theroretical desciritpion, please consult a paper ''Statistical topological data analysis using persistence landscapes'' by Peter Bubenik.
+ * For details of algorithms, please consult ''A persistence landscapes toolbox for topological statistics'' by Peter Bubenik and Pawel Dlotko.
+ * Persistence landscapes allow vertorization, computations of distances, computations of projections to Real, computations of averages and scalar products. Therefore they implement suitable interfaces.
+**/
+class Persistence_landscape :
+								public Abs_Vectorized_topological_data ,
+								public Abs_Topological_data_with_distances,
+								public Abs_Real_valued_topological_data,
+								public Abs_Topological_data_with_averages,
+								public Abs_Topological_data_with_scalar_product
+{
+public:
+	/**
+	 * Default constructor.
+	**/
+    Persistence_landscape()
+    {
+		this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+	}
+
+    /**
+	 * Constructor that takes as an input a vector of birth-death pairs.
+	**/
+    Persistence_landscape( const std::vector< std::pair< double , double > >& p );
+
+    /**
+	 * Assignement operator.
+	**/
+    Persistence_landscape& operator=( const Persistence_landscape& org );
+
+    /**
+	 * Copy constructor.
+	**/
+    Persistence_landscape(const Persistence_landscape&);
+
+    /**
+	 * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed
+	 * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read.
+	**/
+    Persistence_landscape(const char* filename , size_t dimension = 0);
+
+
+
+    /**
+     * This procedure loads a landscape from file. It erase all the data that was previously stored in this landscape.
+    **/
+    void load_landscape_from_file( const char* filename );
+
+
+    /**
+     * The procedure stores a landscape to a file. The file can be later used by a procedure load_landscape_from_file.
+    **/
+    void print_to_file( const char* filename )const;
+
+
+
+	/**
+	 * This function compute integral of the landscape (defined formally as sum of integrals on R of all landscape functions)
+	**/
+    double compute_integral_of_landscape()const;
+
+
+    /**
+	 * This function compute integral of the 'level'-level of a landscape.
+	**/
+    double compute_integral_of_landscape( size_t level )const;
+
+
+    /**
+	 * This function compute integral of the landscape p-th power of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions)
+	**/
+	double compute_integral_of_landscape( double p )const;//this function compute integral of p-th power of landscape.
+
+
+    /**
+     * A function that computes the value of a landscape at a given point. The parameters of the function are: unsigned level and double x.
+     * The procedure will compute the value of the level-landscape at the point x.
+    **/
+    double compute_value_at_a_given_point( unsigned level , double x )const;
+
+    /**
+     * Writing landscape into a stream. A i-th level landscape starts with a string "lambda_i". Then the discontinuity points of the landscapes follows.
+     * Shall those points be joined with lines, we will obtain the i-th landscape function.
+    **/
+    friend std::ostream& operator<<(std::ostream& out, Persistence_landscape& land );
+
+
+
+
+
+	/**
+	 * A function that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape add_two_landscapes ( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 )
+    {
+        return operation_on_pair_of_landscapes< std::plus<double> >(land1,land2);
+    }
+
+    /**
+	 * A function that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape subtract_two_landscapes ( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 )
+    {
+        return operation_on_pair_of_landscapes< std::minus<double> >(land1,land2);
+    }
+
+	/**
+	 * An operator +, that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape operator+( const Persistence_landscape& first , const Persistence_landscape& second )
+    {
+        return add_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator -, that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape operator-( const Persistence_landscape& first , const Persistence_landscape& second )
+    {
+        return subtract_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator * that allows multipilication of a landscape by a real number.
+	**/
+    friend Persistence_landscape operator*( const Persistence_landscape& first , double con )
+    {
+        return first.multiply_lanscape_by_real_number_not_overwrite(con);
+    }
+
+	/**
+	 * An operator * that allows multipilication of a landscape by a real number (order of parameters swapped).
+	**/
+    friend Persistence_landscape operator*( double con , const Persistence_landscape& first  )
+    {
+        return first.multiply_lanscape_by_real_number_not_overwrite(con);
+    }
+
+	/**
+	 * Operator +=. The second parameter is persistnece landwscape.
+	**/
+    Persistence_landscape operator += ( const Persistence_landscape& rhs )
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+	/**
+	 * Operator -=. The second parameter is persistnece landwscape.
+	**/
+    Persistence_landscape operator -= ( const Persistence_landscape& rhs )
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+
+	/**
+	 * Operator *=. The second parameter is a real number by which the y values of all landscape functions are multiplied. The x-values remain unchanged. 
+	**/
+    Persistence_landscape operator *= ( double x )
+    {
+        *this = *this*x;
+        return *this;
+    }
+
+	/**
+	 * Operator /=. The second parameter is a real number.
+	**/
+    Persistence_landscape operator /= ( double x )
+    {
+        if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." );
+        *this = *this * (1/x);
+        return *this;
+    }
+
+	/**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator == ( const Persistence_landscape& rhs  )const;
+
+
+    /**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator != ( const Persistence_landscape& rhs  )const
+    {
+		return !((*this) == rhs);
+	}
+
+
+	/**
+	 * Computations of maximum (y) value of landscape.
+	**/
+    double compute_maximum()const
+    {
+        double maxValue = 0;
+        if ( this->land.size() )
+        {
+            maxValue = -std::numeric_limits<int>::max();
+            for ( size_t i = 0 ; i != this->land[0].size() ; ++i )
+            {
+                if ( this->land[0][i].second > maxValue )maxValue = this->land[0][i].second;
+            }
+        }
+        return maxValue;
+    }
+
+	/**
+	 * Computations of a L^i norm of landscape, where i is the input parameter.
+	**/
+    double compute_norm_of_landscape( double i )
+    {
+        Persistence_landscape l;
+        if ( i != -1 )
+        {
+            return compute_discance_of_landscapes(*this,l,i);
+        }
+        else
+        {
+            return compute_max_norm_discance_of_landscapes(*this,l);
+        }
+    }
+
+	/**
+ 	 * An operator to compute the value of a landscape in the level 'level' at the argument 'x'.
+	**/
+    double operator()(unsigned level,double x)const{return this->compute_value_at_a_given_point(level,x);}
+
+	/**
+	 * Computations of L^{\infty} distance between two landscapes.
+	**/
+    friend double compute_max_norm_discance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second );
+    //friend double compute_max_norm_discance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , unsigned& nrOfLand , double&x , double& y1, double& y2 );
+
+
+	/**
+	 * Computations of L^{p} distance between two landscapes. p is the parameter of the procedure.
+	**/
+    friend double compute_discance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , int p );
+
+
+
+	/**
+	 * Function to compute absolute value of a PL function. The representation of persistence landscapes allow to store general PL-function. When computing distance betwen two landscapes, we compute difference between
+	 * them. In this case, a general PL-function with negative value can appear as a result. Then in order to compute distance, we need to take its absolute value. This is the purpose of this procedure.
+	**/
+    Persistence_landscape abs();
+
+	/**
+	 * Computes the number of landscape functions.
+	**/
+    size_t size()const{return this->land.size(); }
+
+	/**
+	 *  Computate maximal value of lambda-level landscape.
+	**/
+    double find_max( unsigned lambda )const;
+
+	/**
+	 * Function to compute inner (scalar) product of two landscapes.
+	**/
+    friend double compute_inner_product( const Persistence_landscape& l1 , const Persistence_landscape& l2 );
+
+
+
+
+	//concretization of abstract functions:
+
+	/**
+	 * The number of projections to R is defined to the number of nonzero landscape functions. I-th projection is an integral of i-th landscape function over whole R.
+	**/
+    double project_to_R( int number_of_function )
+    {
+		return this->compute_integral_of_landscape( (size_t)number_of_function );
+	}
+
+    std::vector<double> vectorize( int number_of_function )
+    {
+		//TODO, think of something smarter over here
+		std::vector<double> v;
+		if ( (size_t)number_of_function > this->land.size() )
+		{
+			return v;
+		}
+		v.reserve( this->land[number_of_function].size() );
+		for ( size_t i = 0 ; i != this->land[number_of_function].size() ; ++i )
+		{
+			v.push_back( this->land[number_of_function][i].second );
+		}
+		return v;
+	}
+    void compute_average( std::vector< Abs_Topological_data_with_averages* > to_average )
+    {
+		bool dbg = false;
+
+		std::vector< Persistence_landscape* > nextLevelMerge( to_average.size() );
+        for ( size_t i = 0 ; i != to_average.size() ; ++i )
+        {
+            nextLevelMerge[i] = (Persistence_landscape*)to_average[i];
+        }
+        bool is_this_first_level = true;//in the loop, we will create dynamically a unmber of intermediate complexes. We have to clean that up, but we cannot erase the initial andscapes we have
+										//to average. In this case, we simply check if the nextLevelMerge are the input landscapes or the ones created in that loop by usig this extra variable.
+
+        while ( nextLevelMerge.size() != 1 )
+        {
+            if ( dbg ){cerr << "nextLevelMerge.size() : " << nextLevelMerge.size() << endl;}
+            std::vector< Persistence_landscape* > nextNextLevelMerge;
+            nextNextLevelMerge.reserve( to_average.size() );
+            for ( size_t i = 0 ; i < nextLevelMerge.size() ; i=i+2 )
+            {
+                if ( dbg ){cerr << "i : " << i << endl;}
+                Persistence_landscape* l = new Persistence_landscape;
+                if ( i+1 != nextLevelMerge.size() )
+                {
+                    (*l) = (*nextLevelMerge[i])+(*nextLevelMerge[i+1]);
+                }
+                else
+                {
+                    (*l) = *nextLevelMerge[i];
+                }
+                nextNextLevelMerge.push_back( l );
+            }
+            if ( dbg ){cerr << "After this iteration \n";}
+
+            if ( !is_this_first_level )
+            {
+				//deallocate the memory if the vector nextLevelMerge do not consist of the initial landscapes
+				for ( size_t i = 0 ; i != nextLevelMerge.size() ; ++i )
+				{
+					delete nextLevelMerge[i];
+				}
+			}
+			is_this_first_level = false;
+            nextLevelMerge.swap(nextNextLevelMerge);
+        }
+        (*this) = (*nextLevelMerge[0]);
+        (*this) *= 1/( (double)to_average.size() );
+	}
+
+
+    double distance( const Abs_Topological_data_with_distances* second , double power = 1 )
+    {
+		if ( power != -1 )
+		{
+			return compute_discance_of_landscapes( *this , *((Persistence_landscape*)second) , power );
+		}
+		else
+		{
+			return compute_max_norm_discance_of_landscapes( *this , *((Persistence_landscape*)second) );
+		}
+	}
+
+
+	double compute_scalar_product( const Abs_Topological_data_with_scalar_product* second )
+	{
+		return compute_inner_product( (*this) , *((Persistence_landscape*)second) );
+	}
+
+
+	std::vector< std::vector< std::pair<double,double> > > output_for_visualization()
+	{
+		return this->land;
+	}
+	
+	
+	//a function used to create a gnuplot script for visualization of landscapes
+	void plot( const char* filename ,int from = -1, int to = -1 ,  double xRangeBegin = -1 , double xRangeEnd = -1 , double yRangeBegin = -1 , double yRangeEnd = -1 );
+
+
+protected:
+    std::vector< std::vector< std::pair<double,double> > > land;
+
+    void construct_persistence_landscape_from_barcode( const std::vector< std::pair< double , double > > & p );
+    Persistence_landscape multiply_lanscape_by_real_number_not_overwrite( double x )const;
+    void multiply_lanscape_by_real_number_overwrite( double x );
+    template < typename oper > friend Persistence_landscape operation_on_pair_of_landscapes ( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 );
+    friend double compute_maximal_distance_non_symmetric( const Persistence_landscape& pl1, const Persistence_landscape& pl2 );
+
+    void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals()
+	{
+		//warning, this function can be only called after filling in the intervals vector.
+		this->number_of_functions_for_vectorization = this->land.size();
+		this->number_of_functions_for_projections_to_reals = this->land.size();
+	}
+};
+
+
+Persistence_landscape::Persistence_landscape(const Persistence_landscape& oryginal)
+{
+    //std::cerr << "Running copy constructor \n";
+    std::vector< std::vector< std::pair<double,double> > > land( oryginal.land.size() );
+    for ( size_t i = 0 ; i != oryginal.land.size() ; ++i )
+    {
+        land[i].insert( land[i].end() , oryginal.land[i].begin() , oryginal.land[i].end() );
+    }
+    this->land = land;
+    this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+}
+
+
+
+
+
+Persistence_landscape::Persistence_landscape(const char* filename , size_t dimension)
+{
+    bool dbg = false;
+
+    if ( dbg )
+    {
+        std::cerr << "Using constructor : Persistence_landscape(char* filename)" << std::endl;
+    }   
+    //standard file with barcode
+    std::vector< std::pair< double , double > > barcode = read_standard_file( filename );    
+    //gudhi file with barcode
+    //std::vector< std::pair< double , double > > barcode = read_gudhi_file( filename , dimension );        
+	this->construct_persistence_landscape_from_barcode( barcode );
+	this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+}
+
+
+bool operatorEqualDbg = false;
+bool Persistence_landscape::operator == ( const Persistence_landscape& rhs  )const
+{
+    if ( this->land.size() != rhs.land.size() )
+    {
+        if (operatorEqualDbg)std::cerr << "1\n";
+        return false;
+    }
+    for ( size_t level = 0 ; level != this->land.size() ; ++level )
+    {
+        if ( this->land[level].size() != rhs.land[level].size() )
+        {
+            if (operatorEqualDbg)std::cerr << "this->land[level].size() : " << this->land[level].size() <<  "\n";
+            if (operatorEqualDbg)std::cerr << "rhs.land[level].size() : " << rhs.land[level].size() <<  "\n";
+            if (operatorEqualDbg)std::cerr << "2\n";
+            return false;
+        }
+        for ( size_t i = 0 ; i != this->land[level].size() ; ++i )
+        {
+            if ( !( almost_equal(this->land[level][i].first , rhs.land[level][i].first) && almost_equal(this->land[level][i].second , rhs.land[level][i].second) ) )
+            {
+				//cerr<< this->land[level][i].first << " , " <<  rhs.land[level][i].first << " and " << this->land[level][i].second << " , " << rhs.land[level][i].second << endl;
+                if (operatorEqualDbg)std::cerr << "this->land[level][i] : " << this->land[level][i].first << " " << this->land[level][i].second << "\n";
+                if (operatorEqualDbg)std::cerr << "rhs.land[level][i] : " << rhs.land[level][i].first << " " << rhs.land[level][i].second	 << "\n";
+                if (operatorEqualDbg)std::cerr << "3\n";
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+
+Persistence_landscape& Persistence_landscape::operator=( const Persistence_landscape& oryginal )
+{
+    std::vector< std::vector< std::pair<double,double> > > land( oryginal.land.size() );
+    for ( size_t i = 0 ; i != oryginal.land.size() ; ++i )
+    {
+        land[i].insert( land[i].end() , oryginal.land[i].begin() , oryginal.land[i].end() );
+    }
+    this->land = land;
+    return *this;
+}
+
+
+
+Persistence_landscape::Persistence_landscape( const std::vector< std::pair< double , double > > & p )
+{
+    this->construct_persistence_landscape_from_barcode( p );
+    this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+}
+
+
+void Persistence_landscape::construct_persistence_landscape_from_barcode( const std::vector< std::pair< double , double > > & p )
+{
+	bool dbg = false;
+    if ( dbg ){cerr << "Persistence_landscape::Persistence_landscape( const std::vector< std::pair< double , double > >& p )" << endl;}
+
+	//this is a general algorithm to construct persistence landscapes.
+	std::vector< std::pair<double,double> > bars;
+	bars.insert( bars.begin() , p.begin() , p.end() );
+	std::sort( bars.begin() , bars.end() , compare_points_sorting );
+
+	if (dbg)
+	{
+		std::cerr << "Bars : \n";
+		for ( size_t i = 0 ; i != bars.size() ; ++i )
+		{
+			std::cerr << bars[i].first << " " << bars[i].second << "\n";
+		}
+		getchar();
+	}
+
+	std::vector< std::pair<double,double> > characteristicPoints(p.size());
+	for ( size_t i = 0 ; i != bars.size() ; ++i )
+	{
+		characteristicPoints[i] = std::make_pair((bars[i].first+bars[i].second)/2.0 , (bars[i].second - bars[i].first)/2.0);
+	}
+	std::vector< std::vector< std::pair<double,double> > > Persistence_landscape;
+	while ( !characteristicPoints.empty() )
+	{
+		if(dbg)
+		{
+			for ( size_t i = 0 ; i != characteristicPoints.size() ; ++i )
+			{
+				std::cout << "("  << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n";
+			}
+			std::cin.ignore();
+		}
+
+		std::vector< std::pair<double,double> > lambda_n;
+		lambda_n.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) );
+		lambda_n.push_back( std::make_pair(minus_length(characteristicPoints[0]),0) );
+		lambda_n.push_back( characteristicPoints[0] );
+
+		if (dbg)
+		{
+			std::cerr << "1 Adding to lambda_n : (" << -std::numeric_limits<int>::max() << " " << 0  << ") , (" << minus_length(characteristicPoints[0]) << " " << 0 << ") , (" << characteristicPoints[0].first << " " << characteristicPoints[0].second << ") \n";
+		}
+
+		size_t i = 1;
+		std::vector< std::pair<double,double> >  newCharacteristicPoints;
+		while ( i < characteristicPoints.size() )
+		{
+			size_t p = 1;
+			if ( (minus_length(characteristicPoints[i]) >= minus_length(lambda_n[lambda_n.size()-1])) && (birth_plus_deaths(characteristicPoints[i]) > birth_plus_deaths(lambda_n[lambda_n.size()-1])) )
+			{
+				if ( minus_length(characteristicPoints[i]) < birth_plus_deaths(lambda_n[lambda_n.size()-1]) )
+				{
+					std::pair<double,double> point = std::make_pair( (minus_length(characteristicPoints[i])+birth_plus_deaths(lambda_n[lambda_n.size()-1]))/2 , (birth_plus_deaths(lambda_n[lambda_n.size()-1])-minus_length(characteristicPoints[i]))/2 );
+					lambda_n.push_back( point );
+					if (dbg)
+					{
+						std::cerr << "2 Adding to lambda_n : (" << point.first << " " << point.second << ")\n";
+					}
+
+
+					if ( dbg )
+					{
+						std::cerr << "characteristicPoints[i+p] : " << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << "\n";
+						std::cerr << "point : " << point.first << " " << point.second << "\n";
+						getchar();
+					}
+
+					while ( (i+p < characteristicPoints.size() ) && ( almost_equal(minus_length(point),minus_length(characteristicPoints[i+p])) ) && ( birth_plus_deaths(point) <= birth_plus_deaths(characteristicPoints[i+p]) ) )
+					{
+						newCharacteristicPoints.push_back( characteristicPoints[i+p] );
+						if (dbg)
+						{
+							std::cerr << "3.5 Adding to newCharacteristicPoints : (" << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << ")\n";
+							getchar();
+						}
+						++p;
+					}
+
+
+					newCharacteristicPoints.push_back( point );
+					if (dbg)
+					{
+						std::cerr << "4 Adding to newCharacteristicPoints : (" << point.first << " " << point.second << ")\n";
+					}
+
+
+					while ( (i+p < characteristicPoints.size() ) && ( minus_length(point) <= minus_length(characteristicPoints[i+p]) ) && (birth_plus_deaths(point)>=birth_plus_deaths(characteristicPoints[i+p])) )
+					{
+						newCharacteristicPoints.push_back( characteristicPoints[i+p] );
+						if (dbg)
+						{
+							std::cerr << "characteristicPoints[i+p] : " << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << "\n";
+							std::cerr << "point : " << point.first << " " << point.second << "\n";
+							std::cerr << "characteristicPoints[i+p] birth and death : " << minus_length(characteristicPoints[i+p]) << " , " << birth_plus_deaths(characteristicPoints[i+p]) << "\n";
+							std::cerr << "point birth and death : " << minus_length(point) << " , " << birth_plus_deaths(point) << "\n";
+
+							std::cerr << "3 Adding to newCharacteristicPoints : (" << characteristicPoints[i+p].first << " " << characteristicPoints[i+p].second << ")\n";
+							getchar();
+						}
+						++p;
+					}
+
+				}
+				else
+				{
+					lambda_n.push_back( std::make_pair( birth_plus_deaths(lambda_n[lambda_n.size()-1]) , 0 ) );
+					lambda_n.push_back( std::make_pair( minus_length(characteristicPoints[i]) , 0 ) );
+					if (dbg)
+					{
+						std::cerr << "5 Adding to lambda_n : (" << birth_plus_deaths(lambda_n[lambda_n.size()-1]) << " " <<  0 << ")\n";
+						std::cerr << "5 Adding to lambda_n : (" << minus_length(characteristicPoints[i]) << " " << 0  << ")\n";
+					}
+				}
+				lambda_n.push_back( characteristicPoints[i] );
+				if (dbg)
+				{
+					std::cerr << "6 Adding to lambda_n : (" << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n";
+				}
+			}
+			else
+			{
+				newCharacteristicPoints.push_back( characteristicPoints[i] );
+				if (dbg)
+				{
+					std::cerr << "7 Adding to newCharacteristicPoints : (" << characteristicPoints[i].first << " " << characteristicPoints[i].second << ")\n";
+				}
+			}
+			i = i+p;
+		}
+		lambda_n.push_back( std::make_pair(birth_plus_deaths(lambda_n[lambda_n.size()-1]),0) );
+		lambda_n.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) );
+
+		characteristicPoints = newCharacteristicPoints;
+
+		lambda_n.erase(std::unique(lambda_n.begin(), lambda_n.end()), lambda_n.end());
+		this->land.push_back( lambda_n );
+	}
+}
+
+
+
+//this function find maximum of lambda_n
+double Persistence_landscape::find_max( unsigned lambda )const
+{
+    if ( this->land.size() < lambda )return 0;
+    double maximum = -std::numeric_limits<int>::max();
+    for ( size_t i = 0 ; i != this->land[lambda].size() ; ++i )
+    {
+        if ( this->land[lambda][i].second > maximum )maximum = this->land[lambda][i].second;
+    }
+    return maximum;
+}
+
+
+double Persistence_landscape::compute_integral_of_landscape()const
+{
+    double result = 0;
+    for ( size_t i = 0 ; i != this->land.size() ; ++i )
+    {
+        for ( size_t nr = 2 ; nr != this->land[i].size()-1 ; ++nr )
+        {
+            //it suffices to compute every planar integral and then sum them ap for each lambda_n
+            result += 0.5*( this->land[i][nr].first - this->land[i][nr-1].first )*(this->land[i][nr].second + this->land[i][nr-1].second);
+        }
+    }
+    return result;
+}
+
+double Persistence_landscape::compute_integral_of_landscape( size_t  level )const
+{
+    double result = 0;
+    if ( level >= this->land.size() )
+    {
+		//this landscape function is constantly equal 0, so is the intergral.
+		return result;
+	}
+	//also negative landscapes are assumed to be zero.
+	if ( level < 0 )return 0;
+
+	for ( size_t nr = 2 ; nr != this->land[ level ].size()-1 ; ++nr )
+	{
+		//it suffices to compute every planar integral and then sum them ap for each lambda_n
+		result += 0.5*( this->land[ level ][nr].first - this->land[ level ][nr-1].first )*(this->land[ level ][nr].second + this->land[ level ][nr-1].second);
+	}
+
+    return result;
+}
+
+
+double Persistence_landscape::compute_integral_of_landscape( double p )const
+{
+	bool dbg = false;
+    double result = 0;
+    for ( size_t i = 0 ; i != this->land.size() ; ++i )
+    {
+        for ( size_t nr = 2 ; nr != this->land[i].size()-1 ; ++nr )
+        {
+            if (dbg)std::cout << "nr : " << nr << "\n";
+            //In this interval, the landscape has a form f(x) = ax+b. We want to compute integral of (ax+b)^p = 1/a * (ax+b)^{p+1}/(p+1)
+            std::pair<double,double> coef = compute_parameters_of_a_line( this->land[i][nr] , this->land[i][nr-1] );
+            double a = coef.first;
+            double b = coef.second;
+
+            if (dbg)std::cout << "(" << this->land[i][nr].first << "," << this->land[i][nr].second << ") , " << this->land[i][nr-1].first << "," << this->land[i][nr].second << ")" << std::endl;
+            if ( this->land[i][nr].first == this->land[i][nr-1].first )continue;
+            if ( a != 0 )
+            {
+                result += 1/(a*(p+1)) * ( pow((a*this->land[i][nr].first+b),p+1) - pow((a*this->land[i][nr-1].first+b),p+1));
+            }
+            else
+            {
+                result += ( this->land[i][nr].first - this->land[i][nr-1].first )*( pow(this->land[i][nr].second,p) );
+            }
+            if ( dbg )
+            {
+                std::cout << "a : " <<a << " , b : " << b << std::endl;
+                std::cout << "result : " << result << std::endl;
+            }
+        }
+        //if (compute_integral_of_landscapeDbg) std::cin.ignore();
+    }
+    return result;
+}
+
+
+//this is O(log(n)) algorithm, where n is number of points in this->land.
+double Persistence_landscape::compute_value_at_a_given_point( unsigned level , double x )const
+{
+	bool compute_value_at_a_given_pointDbg = false;
+    //in such a case lambda_level = 0.
+    if ( level > this->land.size() ) return 0;
+
+    //we know that the points in this->land[level] are ordered according to x coordinate. Therefore, we can find the point by using bisection:
+    unsigned coordBegin = 1;
+    unsigned coordEnd = this->land[level].size()-2;
+
+    if ( compute_value_at_a_given_pointDbg )
+    {
+        std::cerr << "Tutaj \n";
+        std::cerr << "x : " << x << "\n";
+        std::cerr << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n";
+        std::cerr << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n";
+    }
+
+    //in this case x is outside the support of the landscape, therefore the value of the landscape is 0.
+    if ( x <= this->land[level][coordBegin].first )return 0;
+    if ( x >= this->land[level][coordEnd].first )return 0;
+
+    if (compute_value_at_a_given_pointDbg)std::cerr << "Entering to the while loop \n";
+
+    while ( coordBegin+1 != coordEnd )
+    {
+        if (compute_value_at_a_given_pointDbg)
+        {
+            std::cerr << "coordBegin : " << coordBegin << "\n";
+            std::cerr << "coordEnd : " << coordEnd << "\n";
+            std::cerr << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n";
+            std::cerr << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n";
+        }
+
+
+        unsigned newCord = (unsigned)floor((coordEnd+coordBegin)/2.0);
+
+        if (compute_value_at_a_given_pointDbg)
+        {
+            std::cerr << "newCord : " << newCord << "\n";
+            std::cerr << "this->land[level][newCord].first : " << this->land[level][newCord].first << "\n";
+            std::cin.ignore();
+        }
+
+        if ( this->land[level][newCord].first <= x )
+        {
+            coordBegin = newCord;
+            if ( this->land[level][newCord].first == x )return this->land[level][newCord].second;
+        }
+        else
+        {
+            coordEnd = newCord;
+        }
+    }
+
+    if (compute_value_at_a_given_pointDbg)
+    {
+        std::cout << "x : " << x << " is between : " << this->land[level][coordBegin].first << " a  " << this->land[level][coordEnd].first << "\n";
+        std::cout << "the y coords are : " << this->land[level][coordBegin].second << " a  " << this->land[level][coordEnd].second << "\n";
+        std::cerr << "coordBegin : " << coordBegin << "\n";
+        std::cerr << "coordEnd : " << coordEnd << "\n";
+        std::cin.ignore();
+    }
+    return function_value( this->land[level][coordBegin] , this->land[level][coordEnd] , x );
+}
+
+std::ostream& operator<<(std::ostream& out, Persistence_landscape& land )
+{
+    for ( size_t level = 0 ; level != land.land.size() ; ++level )
+    {
+        out << "Lambda_" << level << ":" << std::endl;
+        for ( size_t i = 0 ; i != land.land[level].size() ; ++i )
+        {
+            if ( land.land[level][i].first == -std::numeric_limits<int>::max() )
+            {
+                out << "-inf";
+            }
+            else
+            {
+                if ( land.land[level][i].first == std::numeric_limits<int>::max() )
+                {
+                    out << "+inf";
+                }
+                else
+                {
+                    out << land.land[level][i].first;
+                }
+            }
+            out << " , " << land.land[level][i].second << std::endl;
+        }
+    }
+    return out;
+}
+
+
+
+
+void Persistence_landscape::multiply_lanscape_by_real_number_overwrite( double x )
+{
+    for ( size_t dim = 0 ; dim != this->land.size() ; ++dim )
+    {
+        for ( size_t i = 0 ; i != this->land[dim].size() ; ++i )
+        {
+             this->land[dim][i].second *= x;
+        }
+    }
+}
+
+bool AbsDbg = false;
+Persistence_landscape Persistence_landscape::abs()
+{
+    Persistence_landscape result;
+    for ( size_t level = 0 ; level != this->land.size() ; ++level )
+    {
+        if ( AbsDbg ){ std::cout << "level: " << level << std::endl; }
+        std::vector< std::pair<double,double> > lambda_n;
+        lambda_n.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) );
+        for ( size_t i = 1 ; i != this->land[level].size() ; ++i )
+        {
+            if ( AbsDbg ){std::cout << "this->land[" << level << "][" << i << "] : " << this->land[level][i].first << " " << this->land[level][i].second << std::endl;}
+            //if a line segment between this->land[level][i-1] and this->land[level][i] crosses the x-axis, then we have to add one landscape point t oresult
+            if ( (this->land[level][i-1].second)*(this->land[level][i].second)  < 0 )
+            {
+                double zero = find_zero_of_a_line_segment_between_those_two_points( this->land[level][i-1] , this->land[level][i] );
+
+                lambda_n.push_back( std::make_pair(zero , 0) );
+                lambda_n.push_back( std::make_pair(this->land[level][i].first , fabs(this->land[level][i].second)) );
+                if ( AbsDbg )
+                {
+                    std::cout << "Adding pair : (" << zero << ",0)" << std::endl;
+                    std::cout << "In the same step adding pair : (" << this->land[level][i].first << "," << fabs(this->land[level][i].second) << ") " << std::endl;
+                    std::cin.ignore();
+                }
+            }
+            else
+            {
+                lambda_n.push_back( std::make_pair(this->land[level][i].first , fabs(this->land[level][i].second)) );
+                if ( AbsDbg )
+                {
+                    std::cout << "Adding pair : (" << this->land[level][i].first << "," << fabs(this->land[level][i].second) << ") " << std::endl;
+                    std::cin.ignore();
+                }
+            }
+        }
+        result.land.push_back( lambda_n );
+    }
+    return result;
+}
+
+
+Persistence_landscape Persistence_landscape::multiply_lanscape_by_real_number_not_overwrite( double x )const
+{
+    std::vector< std::vector< std::pair<double,double> > > result(this->land.size());
+    for ( size_t dim = 0 ; dim != this->land.size() ; ++dim )
+    {
+        std::vector< std::pair<double,double> > lambda_dim( this->land[dim].size() );
+        for ( size_t i = 0 ; i != this->land[dim].size() ; ++i )
+        {
+            lambda_dim[i] = std::make_pair( this->land[dim][i].first , x*this->land[dim][i].second );
+        }
+        result[dim] = lambda_dim;
+    }
+    Persistence_landscape res;
+    //CHANGE
+    //res.land = result;
+    res.land.swap(result);
+    return res;
+}//multiply_lanscape_by_real_number_overwrite
+
+
+void Persistence_landscape::print_to_file( const char* filename )const
+{
+    std::ofstream write;
+    write.open(filename);
+    for ( size_t dim = 0 ; dim != this->land.size() ; ++dim )
+    {
+        write << "#lambda_" << dim << std::endl;
+        for ( size_t i = 1 ; i != this->land[dim].size()-1 ; ++i )
+        {
+            write << this->land[dim][i].first << "  " << this->land[dim][i].second << std::endl;
+        }
+    }
+    write.close();
+}
+
+void Persistence_landscape::load_landscape_from_file( const char* filename )
+{
+	bool dbg = false;
+	//removing the current content of the persistence landscape.
+	this->land.clear();
+
+
+	//this constructor reads persistence landscape form a file. This file have to be created by this software beforehead
+    std::ifstream in;
+    in.open( filename );
+	if ( !( access( filename, F_OK ) != -1 ) )
+	{
+		cerr << "The file : " << filename << " do not exist. The program will now terminate \n";
+		throw "The file from which you are trying to read the persistence landscape do not exist. The program will now terminate \n";
+	}
+
+    std::string line;
+    std::vector< std::pair<double,double> > landscapeAtThisLevel;
+
+    bool isThisAFirsLine = true;
+    while ( !in.eof() )
+    {
+        getline(in,line);
+        if ( !(line.length() == 0 || line[0] == '#') )
+        {
+            std::stringstream lineSS;
+            lineSS << line;
+            double beginn, endd;
+            lineSS >> beginn;
+            lineSS >> endd;
+            landscapeAtThisLevel.push_back( std::make_pair( beginn , endd ) );
+            if (dbg){std::cerr << "Reading a pont : " << beginn << " , " << endd << std::endl;}
+        }
+        else
+        {
+            if (dbg)
+            {
+                std::cout << "IGNORE LINE\n";
+                getchar();
+            }
+            if ( !isThisAFirsLine )
+            {
+                landscapeAtThisLevel.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) );
+                this->land.push_back(landscapeAtThisLevel);
+                std::vector< std::pair<double,double> > newLevelOdLandscape;
+                landscapeAtThisLevel.swap(newLevelOdLandscape);
+            }
+            landscapeAtThisLevel.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) );
+            isThisAFirsLine = false;
+        }
+	}
+	if ( landscapeAtThisLevel.size() > 1 )
+    {
+        //seems that the last line of the file is not finished with the newline sign. We need to put what we have in landscapeAtThisLevel to the constructed landscape.
+        landscapeAtThisLevel.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) );
+        this->land.push_back(landscapeAtThisLevel);
+    }
+
+    in.close();
+}
+
+
+template < typename T >
+Persistence_landscape operation_on_pair_of_landscapes ( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 )
+{
+	bool operation_on_pair_of_landscapesDBG = false;
+    if ( operation_on_pair_of_landscapesDBG ){std::cout << "operation_on_pair_of_landscapes\n";std::cin.ignore();}
+    Persistence_landscape result;
+    std::vector< std::vector< std::pair<double,double> > > land( std::max( land1.land.size() , land2.land.size() ) );
+    result.land = land;
+    T oper;
+
+    for ( size_t i = 0 ; i != std::min( land1.land.size() , land2.land.size() ) ; ++i )
+    {
+        std::vector< std::pair<double,double> > lambda_n;
+        size_t p = 0;
+        size_t q = 0;
+        while ( (p+1 < land1.land[i].size()) && (q+1 < land2.land[i].size()) )
+        {
+            if ( operation_on_pair_of_landscapesDBG )
+            {
+                std::cerr << "p : " << p << "\n";
+                std::cerr << "q : " << q << "\n";
+                std::cout << "land1.land[i][p].first : " << land1.land[i][p].first << "\n";
+                std::cout << "land2.land[i][q].first : " << land2.land[i][q].first << "\n";
+            }
+
+            if ( land1.land[i][p].first < land2.land[i][q].first )
+            {
+                if ( operation_on_pair_of_landscapesDBG )
+                {
+                    std::cout << "first \n";
+                    std::cout << " function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) : "<<  function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) << "\n";
+                    //std::cout << "oper( " << land1.land[i][p].second <<"," << function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) << " : " << oper( land1.land[i][p].second , function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) ) << "\n";
+                }
+                lambda_n.push_back( 
+									std::make_pair( 
+									              land1.land[i][p].first , 
+									              oper( (double)land1.land[i][p].second , function_value(land2.land[i][q-1],land2.land[i][q],land1.land[i][p].first) ) 
+												  ) 
+								  );
+                ++p;
+                continue;
+            }
+            if ( land1.land[i][p].first > land2.land[i][q].first )
+            {
+                if ( operation_on_pair_of_landscapesDBG )
+                {
+                    std::cout << "Second \n";
+                    std::cout << "function_value("<< land1.land[i][p-1].first << " " << land1.land[i][p-1].second <<" ,"<< land1.land[i][p].first << " " << land1.land[i][p].second <<", " << land2.land[i][q].first<<" ) : " << function_value( land1.land[i][p-1] , land1.land[i][p-1] ,land2.land[i][q].first ) << "\n";
+                    std::cout << "oper( " << function_value( land1.land[i][p] , land1.land[i][p-1] ,land2.land[i][q].first ) <<"," << land2.land[i][q].second <<" : " << oper( land2.land[i][q].second , function_value( land1.land[i][p] , land1.land[i][p-1] ,land2.land[i][q].first ) ) << "\n";
+                }
+                lambda_n.push_back( std::make_pair( land2.land[i][q].first , oper( function_value( land1.land[i][p] , land1.land[i][p-1] ,land2.land[i][q].first ) , land2.land[i][q].second )  )  );
+                ++q;
+                continue;
+            }
+            if ( land1.land[i][p].first == land2.land[i][q].first )
+            {
+                if (operation_on_pair_of_landscapesDBG)std::cout << "Third \n";
+                lambda_n.push_back( std::make_pair( land2.land[i][q].first , oper( land1.land[i][p].second , land2.land[i][q].second ) ) );
+                ++p;++q;
+            }
+            if (operation_on_pair_of_landscapesDBG){std::cout << "Next iteration \n";getchar();}
+        }
+        while ( (p+1 < land1.land[i].size())&&(q+1 >= land2.land[i].size()) )
+        {
+            if (operation_on_pair_of_landscapesDBG)
+            {
+                std::cout << "New point : " << land1.land[i][p].first << "  oper(land1.land[i][p].second,0) : " <<  oper(land1.land[i][p].second,0) << std::endl;
+            }
+            lambda_n.push_back( std::make_pair(land1.land[i][p].first , oper(land1.land[i][p].second,0) ) );
+            ++p;
+        }
+        while ( (p+1 >= land1.land[i].size())&&(q+1 < land2.land[i].size()) )
+        {
+            if (operation_on_pair_of_landscapesDBG)
+            {
+                std::cout << "New point : " << land2.land[i][q].first << " oper(0,land2.land[i][q].second) : " <<  oper(0,land2.land[i][q].second) << std::endl;
+            }
+            lambda_n.push_back( std::make_pair(land2.land[i][q].first , oper(0,land2.land[i][q].second) ) );
+            ++q;
+        }
+        lambda_n.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) );
+        //CHANGE
+        //result.land[i] = lambda_n;
+        result.land[i].swap(lambda_n);
+    }
+    if ( land1.land.size() > std::min( land1.land.size() , land2.land.size() ) )
+    {
+        if (operation_on_pair_of_landscapesDBG){std::cout << "land1.land.size() > std::min( land1.land.size() , land2.land.size() )" << std::endl;}
+        for ( size_t i = std::min( land1.land.size() , land2.land.size() ) ; i != std::max( land1.land.size() , land2.land.size() ) ; ++i )
+        {
+            std::vector< std::pair<double,double> > lambda_n( land1.land[i] );
+            for ( size_t nr = 0 ; nr != land1.land[i].size() ; ++nr )
+            {
+                lambda_n[nr] = std::make_pair( land1.land[i][nr].first , oper( land1.land[i][nr].second , 0 ) );
+            }
+            //CHANGE
+            //result.land[i] = lambda_n;
+            result.land[i].swap(lambda_n);
+        }
+    }
+    if ( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) )
+    {
+        if (operation_on_pair_of_landscapesDBG){std::cout << "( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) ) " << std::endl;}
+        for ( size_t i = std::min( land1.land.size() , land2.land.size() ) ; i != std::max( land1.land.size() , land2.land.size() ) ; ++i )
+        {
+            std::vector< std::pair<double,double> > lambda_n( land2.land[i] );
+            for ( size_t nr = 0 ; nr != land2.land[i].size() ; ++nr )
+            {
+                lambda_n[nr] = std::make_pair( land2.land[i][nr].first , oper( 0 , land2.land[i][nr].second ) );
+            }
+            //CHANGE
+            //result.land[i] = lambda_n;
+            result.land[i].swap(lambda_n);
+        }
+    }
+    if ( operation_on_pair_of_landscapesDBG ){std::cout << "operation_on_pair_of_landscapes\n";std::cin.ignore();}
+    return result;
+}//operation_on_pair_of_landscapes
+
+
+
+double compute_maximal_distance_non_symmetric( const Persistence_landscape& pl1, const Persistence_landscape& pl2 )
+{
+    bool dbg = false;
+    if (dbg)std::cerr << " compute_maximal_distance_non_symmetric \n";
+    //this distance is not symmetric. It compute ONLY distance between inflection points of pl1 and pl2.
+    double maxDist = 0;
+    size_t minimalNumberOfLevels = std::min( pl1.land.size() , pl2.land.size() );
+    for ( size_t level = 0 ; level != minimalNumberOfLevels ; ++ level )
+    {
+        if (dbg)
+        {
+            std::cerr << "Level : " << level << std::endl;
+            std::cerr << "PL1 : \n";
+            for ( size_t i = 0 ; i  != pl1.land[level].size() ; ++i )
+            {
+                std::cerr << "(" <<pl1.land[level][i].first << "," << pl1.land[level][i].second << ") \n";
+            }
+            std::cerr << "PL2 : \n";
+            for ( size_t i = 0 ; i  != pl2.land[level].size() ; ++i )
+            {
+                std::cerr << "(" <<pl2.land[level][i].first << "," << pl2.land[level][i].second << ") \n";
+            }
+            std::cin.ignore();
+        }
+
+        int p2Count = 0;
+        for ( size_t i = 1 ; i != pl1.land[level].size()-1 ; ++i ) //w tym przypadku nie rozwarzam punktow w nieskocznosci
+        {
+            while ( true )
+            {
+                if (  (pl1.land[level][i].first>=pl2.land[level][p2Count].first) && (pl1.land[level][i].first<=pl2.land[level][p2Count+1].first)  )break;
+                p2Count++;
+            }
+            double val = fabs( function_value( pl2.land[level][p2Count] , pl2.land[level][p2Count+1] , pl1.land[level][i].first ) - pl1.land[level][i].second);
+            if ( maxDist <= val )maxDist = val;
+
+            if (dbg)
+            {
+                std::cerr << pl1.land[level][i].first <<"in [" << pl2.land[level][p2Count].first << "," <<  pl2.land[level][p2Count+1].first <<"] \n";
+                std::cerr << "pl1[level][i].second : " << pl1.land[level][i].second << std::endl;
+                std::cerr << "function_value( pl2[level][p2Count] , pl2[level][p2Count+1] , pl1[level][i].first ) : " << function_value( pl2.land[level][p2Count] , pl2.land[level][p2Count+1] , pl1.land[level][i].first ) << std::endl;
+                std::cerr << "val : "  << val << std::endl;
+                std::cin.ignore();
+            }
+        }
+    }
+
+    if (dbg)std::cerr << "minimalNumberOfLevels : " << minimalNumberOfLevels << std::endl;
+
+    if ( minimalNumberOfLevels < pl1.land.size() )
+    {
+        for ( size_t level = minimalNumberOfLevels ; level != pl1.land.size() ; ++ level )
+        {
+            for ( size_t i = 0 ; i != pl1.land[level].size() ; ++i )
+            {
+                if (dbg)std::cerr << "pl1[level][i].second  : " << pl1.land[level][i].second << std::endl;
+                if ( maxDist < pl1.land[level][i].second )maxDist = pl1.land[level][i].second;
+            }
+        }
+    }
+    return maxDist;
+}
+
+
+
+
+double compute_discance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , int p )
+{
+	bool dbg = false;
+    //This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a time:
+
+    //first-second :
+    Persistence_landscape lan = first-second;
+
+    //| first-second |:
+    lan = lan.abs();
+    
+    if ( dbg ){cerr << "Abs of difference ; " << lan << endl;getchar();}
+    
+	if ( p != -1 )
+	{
+		//\int_{- \infty}^{+\infty}| first-second |^p
+		double result;
+		if ( p != 1 )
+		{
+			if ( dbg )std::cerr << "Power != 1, compute integral to the power p\n";
+			result = lan.compute_integral_of_landscape( (double)p );
+		}
+		else
+		{
+			if ( dbg )std::cerr << "Power = 1, compute integral \n";
+			result = lan.compute_integral_of_landscape();
+		}
+		//(\int_{- \infty}^{+\infty}| first-second |^p)^(1/p)
+		return pow( result , 1/(double)p );
+	}
+	else
+	{
+		//p == -1
+		if ( dbg )std::cerr << "Power = -1, compute maximum \n";
+		return lan.compute_maximum();
+	}
+}
+
+double compute_max_norm_discance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second )
+{
+    return std::max( compute_maximal_distance_non_symmetric(first,second) , compute_maximal_distance_non_symmetric(second,first) );
+}
+
+
+bool comparePairsForMerging( std::pair< double , unsigned > first , std::pair< double , unsigned > second )
+{
+    return (first.first < second.first);
+}
+
+
+
+
+double compute_inner_product( const Persistence_landscape& l1 , const Persistence_landscape& l2 )
+{
+    bool dbg = false;
+    double result = 0;
+
+    for ( size_t level = 0 ; level != std::min( l1.size() , l2.size() ) ; ++level )
+    {
+        if ( dbg ){cerr << "Computing inner product for a level : " << level << endl;getchar();}
+        if ( l1.land[level].size() * l2.land[level].size() == 0 )continue;
+
+        //endpoints of the interval on which we will compute the inner product of two locally linear functions:
+        double x1 = -std::numeric_limits<int>::max();
+        double x2;
+        if ( l1.land[level][1].first < l2.land[level][1].first )
+        {
+            x2 = l1.land[level][1].first;
+        }
+        else
+        {
+            x2 = l2.land[level][1].first;
+        }
+
+        //iterators for the landscapes l1 and l2
+        size_t l1It = 0;
+        size_t l2It = 0;
+
+        while ( (l1It < l1.land[level].size()-1) && (l2It < l2.land[level].size()-1) )
+        {
+            //compute the value of a inner product on a interval [x1,x2]
+
+            double a,b,c,d;
+
+			if ( l1.land[level][l1It+1].first != l1.land[level][l1It].first )
+			{
+				a = (l1.land[level][l1It+1].second - l1.land[level][l1It].second)/(l1.land[level][l1It+1].first - l1.land[level][l1It].first);
+			}
+			else
+			{
+				a = 0;
+			}
+            b = l1.land[level][l1It].second - a*l1.land[level][l1It].first;
+            if ( l2.land[level][l2It+1].first != l2.land[level][l2It].first )
+            {
+				c = (l2.land[level][l2It+1].second - l2.land[level][l2It].second)/(l2.land[level][l2It+1].first - l2.land[level][l2It].first);
+			}
+			else
+			{
+				c = 0;
+			}
+            d = l2.land[level][l2It].second - c*l2.land[level][l2It].first;
+
+            double contributionFromThisPart
+            =
+            (a*c*x2*x2*x2/3 + (a*d+b*c)*x2*x2/2 + b*d*x2) - (a*c*x1*x1*x1/3 + (a*d+b*c)*x1*x1/2 + b*d*x1);
+
+            result += contributionFromThisPart;
+
+            if ( dbg )
+            {
+                cerr << "[l1.land[level][l1It].first,l1.land[level][l1It+1].first] : " << l1.land[level][l1It].first << " , " << l1.land[level][l1It+1].first << endl;
+                cerr << "[l2.land[level][l2It].first,l2.land[level][l2It+1].first] : " << l2.land[level][l2It].first << " , " << l2.land[level][l2It+1].first << endl;
+                cerr << "a : " << a << ", b : " << b << " , c: " << c << ", d : " << d << endl;
+                cerr << "x1 : " << x1 << " , x2 : " << x2 << endl;
+                cerr << "contributionFromThisPart : " << contributionFromThisPart << endl;
+                cerr << "result : " << result << endl;
+                getchar();
+            }
+
+            //we have two intervals in which functions are constant:
+            //[l1.land[level][l1It].first , l1.land[level][l1It+1].first]
+            //and
+            //[l2.land[level][l2It].first , l2.land[level][l2It+1].first]
+            //We also have an interval [x1,x2]. Since the intervals in the landscapes cover the whole R, then it is clear that x2
+            //is either l1.land[level][l1It+1].first of l2.land[level][l2It+1].first or both. Lets test it.
+            if ( x2 == l1.land[level][l1It+1].first )
+            {
+                if ( x2 == l2.land[level][l2It+1].first )
+                {
+                    //in this case, we increment both:
+                    ++l2It;
+                    if ( dbg ){cerr << "Incrementing both \n";}
+                }
+                else
+                {
+                    if ( dbg ){cerr << "Incrementing first \n";}
+                }
+                ++l1It;
+            }
+            else
+            {
+                //in this case we increment l2It
+                ++l2It;
+                if ( dbg ){cerr << "Incrementing second \n";}
+            }
+            //Now, we shift x1 and x2:
+            x1 = x2;
+            if ( l1.land[level][l1It+1].first < l2.land[level][l2It+1].first )
+            {
+                x2 = l1.land[level][l1It+1].first;
+            }
+            else
+            {
+                x2 = l2.land[level][l2It+1].first;
+            }
+
+        }
+
+    }
+    return result;
+}
+
+
+void Persistence_landscape::plot( const char* filename , int from, int to , double xRangeBegin , double xRangeEnd , double yRangeBegin , double yRangeEnd )
+{
+    //this program create a gnuplot script file that allows to plot persistence diagram.
+    ofstream out;
+
+    std::ostringstream nameSS;
+    nameSS << filename << "_GnuplotScript";
+    std::string nameStr = nameSS.str();
+    out.open( (char*)nameStr.c_str() );
+
+    if ( (xRangeBegin != -1) || (xRangeEnd != -1) || (yRangeBegin != -1) || (yRangeEnd != -1)  )
+    {
+        out << "set xrange [" << xRangeBegin << " : " << xRangeEnd << "]" << endl;
+        out << "set yrange [" << yRangeBegin << " : " << yRangeEnd << "]" << endl;
+    }
+
+    if ( from == -1 ){from = 0;}
+    if ( to == -1 ){to = this->land.size();}
+
+    out << "plot ";
+    for ( size_t lambda= std::min((size_t)from,this->land.size()) ; lambda != std::min((size_t)to,this->land.size()) ; ++lambda )
+    {
+        //out << "     '-' using 1:2 title 'l" << lambda << "' with lp";
+        out << "     '-' using 1:2 notitle with lp";
+        if ( lambda+1 != std::min((size_t)to,this->land.size()) )
+        {
+            out << ", \\";
+        }
+        out << endl;
+    }
+
+    for ( size_t lambda= std::min((size_t)from,this->land.size()) ; lambda != std::min((size_t)to,this->land.size()) ; ++lambda )
+    {
+        for ( size_t i = 1 ; i != this->land[lambda].size()-1 ; ++i )
+        {
+            out << this->land[lambda][i].first << " " << this->land[lambda][i].second << endl;
+        }
+        out << "EOF" << endl;
+    }
+    cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << endl;
+}
+
+
+
+
+}//namespace gudhi stat
+}//namespace gudhi
+
+
+#endif
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h
new file mode 100644
index 00000000..7afd515d
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_landscape_on_grid.h
@@ -0,0 +1,1255 @@
+#ifndef Persistence_landscape_on_grid_H_
+#define Persistence_landscape_on_grid_H_
+
+ 
+//standard include
+#include <iostream>
+#include <vector>
+#include <limits>
+#include <fstream>
+#include <sstream>
+#include <algorithm>
+#include <unistd.h>
+#include <cmath>
+
+
+//gudhi include
+#include <gudhi/abstract_classes/Abs_Vectorized_topological_data.h>
+#include <gudhi/abstract_classes/Abs_Topological_data_with_averages.h>
+#include <gudhi/abstract_classes/Abs_Topological_data_with_distances.h>
+#include <gudhi/abstract_classes/Abs_Real_valued_topological_data.h>
+#include <gudhi/abstract_classes/Abs_Topological_data_with_scalar_product.h>
+#include <gudhi/concretizations/read_persitence_from_file.h>
+
+
+
+using namespace std;
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+/**
+ * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two points.
+**/
+std::pair<double,double> compute_parameters_of_a_line( std::pair<double,double> p1 , std::pair<double,double> p2 )
+{
+    double a = (p2.second-p1.second)/( p2.first - p1.first );
+    double b = p1.second - a*p1.first;
+    return std::make_pair(a,b);
+}
+
+struct greater
+{
+    template<class T>
+    bool operator()(T const &a, T const &b) const { return a > b; }
+};
+
+//double epsi = std::numeric_limits<double>::epsilon();
+double epsi = 0.000005;
+
+
+/**
+ *  A procedure used to compare doubles. Typically gien two doubles A and B, comparing A == B is not good idea. In this case, we use the procedure almostEqual with the epsi defined at
+ *  the top of the file. Setting up the epsi give the user a tolerance on what should be consider equal.
+**/
+inline bool almost_equal( double a , double b )
+{
+    if ( fabs(a-b) < epsi )
+        return true;
+    return false;
+}
+
+
+
+
+class Persistence_landscape_on_grid  :
+								public Abs_Vectorized_topological_data ,
+								public Abs_Topological_data_with_distances,
+								public Abs_Real_valued_topological_data,
+								public Abs_Topological_data_with_averages,
+								public Abs_Topological_data_with_scalar_product
+{
+public:
+	/**
+	 * Default constructor.
+	**/
+    Persistence_landscape_on_grid()
+    {
+		this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+		this->grid_min = this->grid_max = 0;
+	}
+
+    /**
+	 * Constructor that takes as an input a vector of birth-death pairs.
+	**/
+    Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ );
+
+    /**
+	 * Assignement operator.
+	**/
+    Persistence_landscape_on_grid& operator=( const Persistence_landscape_on_grid& org );
+
+    /**
+	 * Copy constructor.
+	**/
+    Persistence_landscape_on_grid(const Persistence_landscape_on_grid&);
+
+    /**
+	 * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed
+	 * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. The additional parameters of this procedure are: ranges of grid, resoltion of a grid
+	 * and the dimension of intervals that are need to be read from a file (in case of Gudhi format files). 
+	**/
+    Persistence_landscape_on_grid(const char* filename , double grid_min_, double grid_max_ , size_t number_of_points_ , size_t dimension_ = 0 );
+        
+    /**
+	 * Constructor that reads persistence intervals from file and creates persistence landscape. The format of the input file is the following: in each line we put birth-death pair. Last line is assumed
+	 * to be empty. Even if the points within a line are not ordered, they will be ordered while the input is read. The additional parameter is the resoution of a grid. The remaning parameters are 
+	 * calculated based on data. 
+	**/
+    Persistence_landscape_on_grid(const char* filename , size_t number_of_points );
+
+
+    /**
+     * This procedure loads a landscape from file. It erase all the data that was previously stored in this landscape.
+    **/
+    void load_landscape_from_file( const char* filename );
+
+
+    /**
+     * The procedure stores a landscape to a file. The file can be later used by a procedure load_landscape_from_file.
+    **/
+    void print_to_file( const char* filename )const;
+
+
+
+	/**
+	 * This function compute integral of the landscape (defined formally as sum of integrals on R of all landscape functions)
+	**/
+    double compute_integral_of_landscape()const
+    {
+		size_t maximal_level = this->number_of_nonzero_levels();
+		double result = 0;
+		for ( size_t i = 0 ; i != maximal_level ; ++i )
+		{
+			result += this->compute_integral_of_landscape(i);
+		}
+		return result;
+	}
+
+
+    /**
+	 * This function compute integral of the 'level'-level of a landscape.
+	**/
+    double compute_integral_of_landscape( size_t level )const
+    {
+		bool dbg = false;
+		double result = 0;		
+		double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1);
+		
+		if ( dbg )
+		{			
+			cerr << "this->grid_max : " << this->grid_max << endl;
+			cerr << "this->grid_min : " << this->grid_min << endl;
+			cerr << "this->values_of_landscapes.size() : " << this->values_of_landscapes.size() << endl;
+			getchar();
+		}
+		
+		
+		double previous_x = this->grid_min-dx;
+		double previous_y = 0;
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			double current_x = previous_x + dx;
+			double current_y = 0;
+			if ( this->values_of_landscapes[i].size() > level )current_y = this->values_of_landscapes[i][level];
+		
+			if ( dbg )
+			{
+				cerr << "this->values_of_landscapes[i].size() : " << this->values_of_landscapes[i].size() << " , level : " << level << endl;	
+				if ( this->values_of_landscapes[i].size() > level )cerr << "this->values_of_landscapes[i][level] : " << this->values_of_landscapes[i][level] << endl;
+				cerr << "previous_y : " << previous_y << endl;
+				cerr << "current_y : " << current_y << endl;
+				cerr << "dx : " << dx << endl;
+				cerr << "0.5*dx*( previous_y + current_y ); " << 0.5*dx*( previous_y + current_y ) << endl;
+			}
+			
+			result += 0.5*dx*( previous_y + current_y );			
+			previous_x = current_x;
+			previous_y = current_y;
+		}
+		return result;
+	}
+
+	/**
+	 * This function compute integral of the landscape p-th power of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions)
+	**/
+	double compute_integral_of_landscape( double p )const
+	{
+		size_t maximal_level = this->number_of_nonzero_levels();
+		double result = 0;
+		for ( size_t i = 0 ; i != maximal_level ; ++i )
+		{
+			result += this->compute_integral_of_landscape(p,i);
+		}
+		return result;
+	}
+
+    /**
+	 * This function compute integral of the landscape p-th power of a level of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions)
+	**/	
+	double compute_integral_of_landscape( double p , size_t level )const
+	{
+		bool dbg = false;
+		
+		double result = 0;				
+		double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1);		
+		double previous_x = this->grid_min;
+		double previous_y = 0;
+		if ( this->values_of_landscapes[0].size() > level )previous_y = this->values_of_landscapes[0][level];
+		
+		if ( dbg )
+		{
+			cerr << "dx : " << dx << endl;
+			cerr << "previous_x : " << previous_x << endl;
+			cerr << "previous_y : " << previous_y << endl;
+			cerr << "power : " << p << endl;
+			getchar();
+		}
+		
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			double current_x = previous_x + dx;
+			double current_y = 0;
+			if ( this->values_of_landscapes[i].size() > level )current_y = this->values_of_landscapes[i][level];
+			
+			if ( dbg )cerr << "current_y : " << current_y << endl;
+			
+			if ( current_y == previous_y )continue; 
+			
+			std::pair<double,double> coef = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y ) , std::make_pair( current_x , current_y ) );
+			double a = coef.first;
+			double b = coef.second;
+			
+			if ( dbg )
+			{
+				cerr << "A line passing through points : (" << previous_x << "," << previous_y << ") and (" << current_x << "," << current_y << ") is : " << a << "x+" << b << endl;
+			}
+			
+			//In this interval, the landscape has a form f(x) = ax+b. We want to compute integral of (ax+b)^p = 1/a * (ax+b)^{p+1}/(p+1)
+			double value_to_add = 0;
+			if ( a != 0 )
+			{			
+				value_to_add = 1/(a*(p+1)) * ( pow((a*current_x+b),p+1) - pow((a*previous_x+b),p+1));				
+			}
+			else
+			{
+				value_to_add = ( current_x - previous_x )*( pow(b,p) );
+			}				
+			result += value_to_add;
+			if ( dbg )
+			{
+				cerr << "Increasing result by : " << value_to_add << endl;
+				cerr << "restult : " << result << endl;
+				getchar();
+			}		
+			previous_x = current_x;
+			previous_y = current_y;
+		}		
+		if ( dbg )cerr << "The total result is : " << result << endl;
+		return result;
+	}
+	
+	/**
+     * Writing landscape into a stream. A i-th level landscape starts with a string "lambda_i". Then the discontinuity points of the landscapes follows.
+     * Shall those points be joined with lines, we will obtain the i-th landscape function.
+    **/
+    friend std::ostream& operator<<(std::ostream& out, const Persistence_landscape_on_grid& land )
+	{
+		double dx = (land.grid_max - land.grid_min)/(double)(land.values_of_landscapes.size()-1);
+		double x = land.grid_min;
+		for ( size_t i = 0 ; i != land.values_of_landscapes.size() ; ++i )
+		{
+			out << x << " : ";
+			for ( size_t j = 0 ; j != land.values_of_landscapes[i].size() ; ++j )
+			{
+				out << land.values_of_landscapes[i][j] << " ";
+			}
+			out << endl;
+			x += dx;
+		}
+		return out;
+	}
+
+
+    /**
+     * A function that computes the value of a landscape at a given point. The parameters of the function are: unsigned level and double x.
+     * The procedure will compute the value of the level-landscape at the point x.
+    **/
+    double compute_value_at_a_given_point( unsigned level , double x )const
+    {
+		bool dbg = false;
+		if ( (x < this->grid_min) || (x > this->grid_max) )return 0;
+		
+		//find a position of a vector closest to x:
+		double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1);
+		size_t position = size_t((x-this->grid_min)/dx);
+		
+		if ( dbg )
+		{
+			std::cerr << "This is a procedure compute_value_at_a_given_point \n";
+			std::cerr << "level : " << level << endl;
+			std::cerr << "x : " << x << endl;
+			std::cerr << "psoition : " << position << endl;
+		}
+		//check if we are not exacly in the grid point:
+		if ( almost_equal( position*dx+ this->grid_min ,  x) )
+		{
+			if ( this->values_of_landscapes[position].size() < level )
+			{
+				return this->values_of_landscapes[position][level];
+			}
+			else
+			{
+				return 0;
+			}
+		}
+		//in the other case, approximate with a line:
+		std::pair<double,double> line;
+		if ( (this->values_of_landscapes[position].size() > level) && (this->values_of_landscapes[position+1].size() > level) )
+		{
+			line = compute_parameters_of_a_line( std::make_pair( position*dx+ this->grid_min , this->values_of_landscapes[position][level] ) , std::make_pair( (position+1)*dx+ this->grid_min , this->values_of_landscapes[position+1][level] ) );
+		}
+		else
+		{
+			if ( (this->values_of_landscapes[position].size() > level) || (this->values_of_landscapes[position+1].size() > level) )
+			{
+				if ( (this->values_of_landscapes[position].size() > level) )
+				{
+						line = compute_parameters_of_a_line( std::make_pair( position*dx+ this->grid_min , this->values_of_landscapes[position][level] ) , std::make_pair( (position+1)*dx+ this->grid_min , 0 ) );
+				}
+				else
+				{
+					//(this->values_of_landscapes[position+1].size() > level)
+					line = compute_parameters_of_a_line( std::make_pair( position*dx+ this->grid_min , 0 ) , std::make_pair( (position+1)*dx+ this->grid_min , this->values_of_landscapes[position+1][level] ) );	
+				}
+			}
+			else
+			{
+				return 0;
+			}
+		}
+		//compute the value of the linear function parametrized by line on a point x:
+		return line.first*x+line.second;
+	}
+    
+
+	/**
+	 * A function that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid add_two_landscapes ( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 )
+    {
+        return operation_on_pair_of_landscapes_on_grid< std::plus<double> >(land1,land2);
+    }
+
+    /**
+	 * A function that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid subtract_two_landscapes ( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 )
+    {
+        return operation_on_pair_of_landscapes_on_grid< std::minus<double> >(land1,land2);
+    }
+
+	/**
+	 * An operator +, that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid operator+( const Persistence_landscape_on_grid& first , const Persistence_landscape_on_grid& second )
+    {
+        return add_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator -, that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid operator-( const Persistence_landscape_on_grid& first , const Persistence_landscape_on_grid& second )
+    {
+        return subtract_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator * that allows multipilication of a landscape by a real number.
+	**/
+    friend Persistence_landscape_on_grid operator*( const Persistence_landscape_on_grid& first , double con )
+    {
+        return first.multiply_lanscape_by_real_number_not_overwrite(con);
+    }
+
+	/**
+	 * An operator * that allows multipilication of a landscape by a real number (order of parameters swapped).
+	**/
+    friend Persistence_landscape_on_grid operator*( double con , const Persistence_landscape_on_grid& first  )
+    {
+        return first.multiply_lanscape_by_real_number_not_overwrite(con);
+    }
+
+	friend bool check_if_defined_on_the_same_domain( const Persistence_landscape_on_grid& land1, const Persistence_landscape_on_grid& land2 )
+	{
+		if ( land1.values_of_landscapes.size() != land2.values_of_landscapes.size() )return false;
+		if ( land1.grid_min != land2.grid_min )return false;
+		if ( land1.grid_max != land2.grid_max )return false;
+		return true;
+	}
+
+	/**
+	 * Operator +=. The second parameter is persistnece landwscape.
+	**/
+    Persistence_landscape_on_grid operator += ( const Persistence_landscape_on_grid& rhs )
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+	/**
+	 * Operator -=. The second parameter is persistnece landwscape.
+	**/
+    Persistence_landscape_on_grid operator -= ( const Persistence_landscape_on_grid& rhs )
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+
+	/**
+	 * Operator *=. The second parameter is a real number by which the y values of all landscape functions are multiplied. The x-values remain unchanged. 
+	**/
+    Persistence_landscape_on_grid operator *= ( double x )
+    {
+        *this = *this*x;
+        return *this;
+    }
+
+	/**
+	 * Operator /=. The second parameter is a real number.
+	**/
+    Persistence_landscape_on_grid operator /= ( double x )
+    {
+        if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." );
+        *this = *this * (1/x);
+        return *this;
+    }
+
+	/**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator == ( const Persistence_landscape_on_grid& rhs  )const
+    {
+		bool dbg = true;
+		if ( ! this->values_of_landscapes.size() == rhs.values_of_landscapes.size() )
+		{
+			if (dbg) cerr << "values_of_landscapes of incompatable sizes\n";
+			return false;
+		}
+		if ( !almost_equal( this->grid_min , rhs.grid_min ) )
+		{
+			if (dbg) cerr << "grid_min not equal\n";
+			return false; 
+		}
+		if ( !almost_equal(this->grid_max,rhs.grid_max ) )
+		{	
+			if (dbg) cerr << "grid_max not equal\n";
+			return false;		
+		}
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{					
+			for ( size_t aa = 0 ; aa != this->values_of_landscapes[i].size() ; ++aa )
+			{
+				if ( !almost_equal( this->values_of_landscapes[i][aa]  , rhs.values_of_landscapes[i][aa]  ) )
+				{					
+					if (dbg) 
+					{
+						cerr << "Problem in the position : " << i << " of values_of_landscapes. \n";
+						cerr << this->values_of_landscapes[i][aa]   << " " <<  rhs.values_of_landscapes[i][aa] << endl;
+					}
+					return false;
+				}
+			}				
+		}	
+		return true;
+	}
+
+
+    /**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator != ( const Persistence_landscape_on_grid& rhs  )const
+    {
+		return !((*this) == rhs);
+	}
+
+
+	/**
+	 * Computations of maximum (y) value of landscape.
+	**/
+    double compute_maximum()const
+    {
+       //since the function can only be entirely positive or negative, the maximal value will be an extremal value in the arrays:
+		double max_value = -std::numeric_limits<double>::max();
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			if ( this->values_of_landscapes[i].size() )
+			{
+				if ( this->values_of_landscapes[i][0] > max_value )max_value = this->values_of_landscapes[i][0];
+				if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] > max_value )max_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ];
+			}
+		}
+		return max_value;        
+    }
+    
+    /**
+	 * Computations of minimum and maximum value of landscape.
+	**/
+    std::pair<double,double> compute_minimum_maximum()const
+    {        
+       //since the function can only be entirely positive or negative, the maximal value will be an extremal value in the arrays:
+		double max_value = -std::numeric_limits<double>::max();
+		double min_value = 0;
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			if ( this->values_of_landscapes[i].size() )
+			{
+				if ( this->values_of_landscapes[i][0] > max_value )max_value = this->values_of_landscapes[i][0];
+				if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] > max_value )max_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ];
+				
+				if ( this->values_of_landscapes[i][0] < min_value )min_value = this->values_of_landscapes[i][0];
+				if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] < min_value )min_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ];
+			}
+		}
+		return std::make_pair(min_value , max_value);
+    }
+    
+    /**
+     * This function computes maximal lambda for which lambda-level landscape is nonzero.
+    **/
+    size_t number_of_nonzero_levels()const
+    {
+		size_t result = 0;
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			if ( this->values_of_landscapes[i].size() > result )result = this->values_of_landscapes[i].size();
+		}
+		return result;
+	}
+
+	/**
+	 * Computations of a L^i norm of landscape, where i is the input parameter.
+	**/
+    double compute_norm_of_landscape( double i )
+    {
+		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 != -1 )
+        {
+            return compute_discance_of_landscapes_on_grid(*this,l,i);
+        }
+        else
+        {
+            return compute_max_norm_discance_of_landscapes(*this,l);
+        }
+    }
+
+	/**
+ 	 * An operator to compute the value of a landscape in the level 'level' at the argument 'x'.
+	**/
+    double operator()(unsigned level,double x)const{return this->compute_value_at_a_given_point(level,x);}
+
+	/**
+	 * Computations of L^{\infty} distance between two landscapes.
+	**/
+    friend double compute_max_norm_discance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second );
+    //friend double compute_max_norm_discance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second , unsigned& nrOfLand , double&x , double& y1, double& y2 );
+
+
+
+
+	/**
+	 * Function to compute absolute value of a PL function. The representation of persistence landscapes allow to store general PL-function. When computing distance betwen two landscapes, we compute difference between
+	 * them. In this case, a general PL-function with negative value can appear as a result. Then in order to compute distance, we need to take its absolute value. This is the purpose of this procedure.
+	**/
+    void abs()
+    {
+		//Be careful here. We assume that the functions are either entirely positive or negative. They do not change signs. That is why I can implemnt abs in the way presented below:
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			for ( size_t j = 0 ; j != this->values_of_landscapes[i].size() ; ++j )
+			{
+				this->values_of_landscapes[i][j] = std::abs( this->values_of_landscapes[i][j] );
+			}
+		}
+	}
+
+	/**
+	 * Computes the number of landscape functions.
+	**/
+    size_t size()const{return this->number_of_nonzero_levels(); }
+
+	/**
+	 *  Computate maximal value of lambda-level landscape.
+	**/
+    double find_max( unsigned lambda )const
+    {		
+		double max_value = -std::numeric_limits<double>::max();
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			if ( this->values_of_landscapes[i].size() > lambda )
+			{
+				if ( this->values_of_landscapes[i][lambda] > max_value )max_value = this->values_of_landscapes[i][lambda];				
+			}
+		}
+		return max_value;
+	}
+
+	/**
+	 * Function to compute inner (scalar) product of two landscapes.
+	**/
+    friend double compute_inner_product( const Persistence_landscape_on_grid& l1 , const Persistence_landscape_on_grid& l2 )
+    {
+		if ( !check_if_defined_on_the_same_domain(l1,l2) )throw "Landscapes are not defined on the same grid, the program will now terminate";
+		size_t maximal_level = l1.number_of_nonzero_levels();
+		double result = 0;
+		for ( size_t i = 0 ; i != maximal_level ; ++i )
+		{
+			result += compute_inner_product(l1,l2,i);
+		}
+		return result;
+	}
+
+	
+	
+	/**
+	 * Function to compute inner (scalar) product of given levels of two landscapes.
+	**/
+    friend double compute_inner_product( const Persistence_landscape_on_grid& l1 , const Persistence_landscape_on_grid& l2 , size_t level )
+    {
+		bool dbg = false;
+		
+		if ( !check_if_defined_on_the_same_domain(l1,l2) )throw "Landscapes are not defined on the same grid, the program will now terminate";
+		double result = 0;
+		
+		double dx = (l1.grid_max - l1.grid_min)/(double)(l1.values_of_landscapes.size()-1);
+		
+		double previous_x = l1.grid_min-dx;
+		double previous_y_l1 = 0;
+		double previous_y_l2 = 0;
+		for ( size_t i = 0 ; i != l1.values_of_landscapes.size() ; ++i )
+		{
+			if ( dbg )std::cerr << "i : " << i << std::endl;
+			
+			double current_x = previous_x + dx;
+			double current_y_l1 = 0;
+			if ( l1.values_of_landscapes[i].size() > level )current_y_l1 = l1.values_of_landscapes[i][level];
+			
+			double current_y_l2 = 0;
+			if ( l2.values_of_landscapes[i].size() > level )current_y_l2 = l2.values_of_landscapes[i][level];
+			
+			if ( dbg )
+			{
+				std::cerr << "previous_x  : " << previous_x << std::endl;
+				std::cerr << "previous_y_l1 : " << previous_y_l1 << std::endl;
+				std::cerr << "current_y_l1 : " << current_y_l1 << std::endl;
+				std::cerr << "previous_y_l2 : " << previous_y_l2 << std::endl;
+				std::cerr << "current_y_l2 : " << current_y_l2 << std::endl;
+			}
+			
+			std::pair<double,double> l1_coords = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y_l1 ) , std::make_pair( current_x , current_y_l1 ) );			
+			std::pair<double,double> l2_coords = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y_l2 ) , std::make_pair( current_x , current_y_l2 ) );
+			
+			//let us assume that the first line is of a form y = ax+b, and the second one is of a form y = cx + d. Then here are a,b,c,d:
+			double a = l1_coords.first;
+			double b = l1_coords.second;
+			
+			double c = l2_coords.first;
+			double d = l2_coords.second;
+			
+			if ( dbg )
+			{
+				std::cerr << "Here are the formulas for a line: \n";
+				std::cerr << "a : " << a << std::endl;
+				std::cerr << "b : " << b << std::endl;
+				std::cerr << "c : " << c << std::endl;
+				std::cerr << "d : " << d << std::endl;				
+			}
+			
+			//now, to compute the inner product in this interval we need to compute the integral of (ax+b)(cx+d) = acx^2 + (ad+bc)x + bd in the interval from previous_x to current_x:
+			//The integal is ac/3*x^3 + (ac+bd)/2*x^2 + bd*x
+		
+			double added_value = (a*c/3*current_x*current_x*current_x + (a*d+b*c)/2*current_x*current_x + b*d*current_x)-
+			          (a*c/3*previous_x*previous_x*previous_x + (a*d+b*c)/2*previous_x*previous_x + b*d*previous_x);
+
+			if ( dbg )
+			{
+				std::cerr << "Value of the integral on the left end ie : " << previous_x << " is : " <<  a*c/3*previous_x*previous_x*previous_x + (a*d+b*c)/2*previous_x*previous_x + b*d*previous_x << std::endl;
+				std::cerr << "Value of the integral on the right end i.e. : " << current_x << " is " << a*c/3*current_x*current_x*current_x + (a*d+b*c)/2*current_x*current_x + b*d*current_x << std::endl;
+			}
+		
+			result += added_value;
+			
+			if ( dbg )
+			{
+				std::cerr << "added_value : " << added_value << std::endl;
+				std::cerr << "result : " << result << std::endl;
+				getchar();
+			}
+			
+			
+			previous_x = current_x;
+			previous_y_l1 = current_y_l1;
+			previous_y_l2 = current_y_l2;
+			
+		}		
+		return result;
+	}
+
+
+
+
+	//concretization of abstract functions:
+
+	/**
+	 * The number of projections to R is defined to the number of nonzero landscape functions. I-th projection is an integral of i-th landscape function over whole R.
+	**/
+    double project_to_R( int number_of_function )
+    {
+		return this->compute_integral_of_landscape( (size_t)number_of_function );
+	}
+
+    std::vector<double> vectorize( int number_of_function )
+    {
+		//TODO, think of something smarter over here
+		if ( ( number_of_function < 0 ) || ( (size_t)number_of_function >= this->values_of_landscapes.size() ) )
+		{
+			throw "Wrong number of function\n";
+		}
+		std::vector<double> v = this->values_of_landscapes[ number_of_function ];
+		return v;
+	}
+	
+	/**
+	 * A function to compute averaged persistence landscape on a grid, based on vector of persistence landscapes on grid (projected to Abs_Topological_data_with_averages*, since this is a concretization of a virtual method.
+	**/  
+    void compute_average( std::vector< Abs_Topological_data_with_averages* > to_average )
+    {
+		
+		bool dbg = false;
+		//After execution of this procedure, the average is supposed to be in the current object. To make sure that this is the case, we need to do some cleaning first.
+		this->values_of_landscapes.clear();
+		this->grid_min = this->grid_max = 0;
+		
+		//if there is nothing to averate, then the average is a zero landscape. 
+		if ( to_average.size() == 0 )return;
+		
+		//now we need to check if the grids in all objects of to_average are the same:
+		for ( size_t i = 0 ; i != to_average.size() ; ++i )
+		{	
+			if ( !check_if_defined_on_the_same_domain(*((Persistence_landscape_on_grid*)(to_average[0])),*((Persistence_landscape_on_grid*)(to_average[i]))) )throw "Two grids are not compatible"; 				
+		}
+		
+		this->values_of_landscapes = std::vector< std::vector<double> >( ((Persistence_landscape_on_grid*)(to_average[0]))->values_of_landscapes.size() );
+		this->grid_min = ((Persistence_landscape_on_grid*)(to_average[0]))->grid_min;
+		this->grid_max = ((Persistence_landscape_on_grid*)(to_average[0]))->grid_max;
+		
+		if ( dbg )
+		{
+			cerr << "Computations of average. The data from the current landscape have been cleared. We are ready to do the computations. \n";
+		}
+						
+		//for every point in the grid:
+		for ( size_t grid_point = 0 ; grid_point != ((Persistence_landscape_on_grid*)(to_average[0]))->values_of_landscapes.size() ; ++grid_point )
+		{
+			
+			//set up a vector of the correct size:
+			size_t maximal_size_of_vector = 0;
+			for ( size_t land_no = 0 ; land_no != to_average.size() ; ++land_no )
+			{
+				if ( ((Persistence_landscape_on_grid*)(to_average[land_no]))->values_of_landscapes[grid_point].size() > maximal_size_of_vector )
+				maximal_size_of_vector = ((Persistence_landscape_on_grid*)(to_average[land_no]))->values_of_landscapes[grid_point].size();
+			}
+			this->values_of_landscapes[grid_point] = std::vector<double>( maximal_size_of_vector );
+			
+			if ( dbg )
+			{
+				cerr << "We are considering the point : " << grid_point << " of the grid. In this point, there are at most : " << maximal_size_of_vector << " nonzero landscape functions \n";
+			}
+		
+			//and compute an arythmetic average:
+			for ( size_t land_no = 0 ; land_no != to_average.size() ; ++land_no )
+			{
+				//summing:				
+				for ( size_t i = 0 ; i != ((Persistence_landscape_on_grid*)(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] += ((Persistence_landscape_on_grid*)(to_average[land_no]))->values_of_landscapes[grid_point][i];
+				}								
+			}
+			//normalizing:
+			for ( size_t i = 0 ; i != this->values_of_landscapes[grid_point].size() ; ++i )
+			{
+				this->values_of_landscapes[grid_point][i] /= (double)to_average.size();
+			}			 
+		}
+	}//compute_average
+	
+	
+	/**
+	 * Computations of L^{p} distance between two landscapes on a grid. p is the parameter of the procedure.
+	**/
+	friend double compute_discance_of_landscapes_on_grid( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second , int p )
+	{		
+		bool dbg = false;
+		//This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a time:
+
+		if ( dbg )
+		{
+			cerr << "first : " << first << endl;
+			cerr << "second : " << second << endl;
+			getchar();
+		}
+
+		//first-second :
+		Persistence_landscape_on_grid lan = first-second;
+		
+		if ( dbg )
+		{
+			cerr << "Difference : " << lan << endl;
+		}
+
+		//| first-second |:
+		lan.abs();
+		
+		if ( dbg )
+		{
+			cerr << "Abs : " << lan << endl;
+		}
+		
+		if ( p != -1 )
+		{
+			//\int_{- \infty}^{+\infty}| first-second |^p
+			double result;
+			if ( p != 1 )
+			{
+				if (dbg){cerr << "p : " << p << endl; getchar();}
+				result = lan.compute_integral_of_landscape( (double)p );				
+				if (dbg){cerr << "integral : " << result << endl;getchar();}
+			}
+			else
+			{
+				result = lan.compute_integral_of_landscape();			
+				if (dbg){cerr << "integral, wihtout power : " << result << endl;getchar();}
+			}
+			//(\int_{- \infty}^{+\infty}| first-second |^p)^(1/p)
+			return pow( result , 1/(double)p );
+		}
+		else
+		{
+			//p == -1
+			return lan.compute_maximum();
+		}
+	}
+
+	/**
+	* A function to compute distance between persistence landscape on a grid.
+	* The parameter of this functionis a persistence landscapes on grid (projected to Abs_Topological_data_with_averages*, since this is a concretization of a virtual method.
+	**/
+    double distance( const Abs_Topological_data_with_distances* second , double power = 1 )
+    {
+		if ( power != -1 )
+		{
+			return compute_discance_of_landscapes_on_grid( *this , *((Persistence_landscape_on_grid*)second) , power );
+		}
+		else
+		{
+			return compute_max_norm_discance_of_landscapes( *this , *((Persistence_landscape_on_grid*)second) );
+		}
+	}
+
+	/**
+	* A function to compute scalar product of persistence landscape on a grid.
+	* The parameter of this functionis a persistence landscapes on grid (projected to Abs_Topological_data_with_averages*, since this is a concretization of a virtual method.
+	**/
+	double compute_scalar_product( const Abs_Topological_data_with_scalar_product* second )
+	{
+		return compute_inner_product( (*this) , *((Persistence_landscape_on_grid*)second) );
+	}
+
+
+	/**
+	* A function that returns values of landsapes. It can be used for vizualization
+	**/
+	std::vector< std::vector< double > > output_for_visualization()
+	{
+		return this->values_of_landscapes;
+	}
+	
+	/**
+	* function used to create a gnuplot script for visualization of landscapes
+	**/ 
+	void plot( const char* filename , size_t from_ = std::numeric_limits<size_t>::max(), size_t to_ = std::numeric_limits<size_t>::max() );
+
+
+protected:
+	double grid_min;
+	double grid_max;
+	std::vector< std::vector< double > > values_of_landscapes;
+	
+	void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals()
+	{
+		//warning, this function can be only called after filling in the values_of_landscapes vector.
+		this->number_of_functions_for_vectorization = this->values_of_landscapes.size();
+		this->number_of_functions_for_projections_to_reals = this->values_of_landscapes.size();
+	}
+	
+	void set_up_values_of_landscapes( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ );
+	template < typename oper > friend Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 );
+	Persistence_landscape_on_grid multiply_lanscape_by_real_number_not_overwrite( double x )const;
+};
+
+
+void Persistence_landscape_on_grid::set_up_values_of_landscapes( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ )
+{
+	bool dbg = false;
+	if ( dbg )
+	{
+		std::cerr << "Here is the procedure : set_up_values_of_landscapes. The parameters are : grid_min_ : " << grid_min_ << ", grid_max_ : " << grid_max_ << ", number_of_points_ : " << number_of_points_ << endl;
+		//getchar();
+		std::cerr << "Here are the intervals at our disposal : \n";
+		for ( size_t i = 0 ; i != p.size() ; ++i )
+		{
+			std::cerr << p[i].first << " , " << p[i].second << endl;
+		} 
+	}
+	
+	if ( (grid_min_ == std::numeric_limits<double>::max()) || (grid_min_ == std::numeric_limits<double>::max()) )
+	{
+		//in this case, we need to find grid_min_ and grid_min_ based 
+	}
+	
+	this->values_of_landscapes = std::vector< std::vector< double > >( number_of_points_+1 );
+	this->grid_min = grid_min_;
+	this->grid_max = grid_max_;
+	
+	if ( grid_max_ <= grid_min_ )
+	{
+		throw "Wrong parameters of grid_min and grid_max given to the procedure. THe grid have negative, or zero size. The program will now terminate.\n";
+	}
+	
+	double dx = ( grid_max_ - grid_min_ )/(double)(number_of_points_);		
+	//for every interval in the diagram:
+	for ( size_t int_no = 0 ; int_no != p.size() ; ++int_no )
+	{		
+		size_t grid_interval_begin = (p[int_no].first-grid_min_)/dx;
+		size_t grid_interval_end = (p[int_no].second-grid_min_)/dx;
+		size_t grid_interval_midpoint = (size_t)(0.5*(grid_interval_begin+grid_interval_end));
+		
+		if ( dbg )
+		{
+			cerr << "Considering an interval : " << p[int_no].first << "," << p[int_no].second << endl;
+		
+			std::cerr << "grid_interval_begin : " << grid_interval_begin << std::endl;
+			std::cerr << "grid_interval_end : " << grid_interval_end << std::endl;
+			std::cerr << "grid_interval_midpoint : " << grid_interval_midpoint << std::endl;	
+		}
+		
+		double landscape_value = dx;
+		for ( size_t i = grid_interval_begin+1 ; i < grid_interval_midpoint ; ++i )
+		{
+			if ( dbg )
+			{
+				std::cerr << "Adding landscape value (going up) for a point : " << i << " equal : " << landscape_value << std::endl;
+			}
+			this->values_of_landscapes[i].push_back( landscape_value );
+			landscape_value += dx;
+		}
+		for ( size_t i = grid_interval_midpoint ; i <= grid_interval_end ; ++i )
+		{						
+			if ( landscape_value > 0 ) 
+			{ 
+				this->values_of_landscapes[i].push_back( landscape_value );	
+				if ( dbg )
+				{				
+					std::cerr << "AAdding landscape value (going down) for a point : " << i << " equal : " << landscape_value << std::endl;
+				}
+			}
+			landscape_value -= dx;	
+		}  
+	}
+	
+	//and now we need to sort the valuesL
+	for ( size_t pt = 0 ; pt != this->values_of_landscapes.size() ; ++pt )
+	{
+		std::sort( this->values_of_landscapes[pt].begin() , this->values_of_landscapes[pt].end() , greater() );	
+	}
+}//set_up_values_of_landscapes
+
+Persistence_landscape_on_grid::Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ )
+{
+	this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ );
+}//Persistence_landscape_on_grid
+
+Persistence_landscape_on_grid& Persistence_landscape_on_grid::operator=( const Persistence_landscape_on_grid& org )
+{
+	this->grid_min = org.grid_min;
+	this->grid_max = org.grid_max;
+	this->values_of_landscapes = org.values_of_landscapes;
+	return (*this);
+}//operator=
+
+Persistence_landscape_on_grid::Persistence_landscape_on_grid(const Persistence_landscape_on_grid& org)
+{
+	this->grid_min = org.grid_min;
+	this->grid_max = org.grid_max;
+	this->values_of_landscapes = org.values_of_landscapes;
+}//copy constructor
+
+Persistence_landscape_on_grid::Persistence_landscape_on_grid(const char* filename , double grid_min_, double grid_max_ , size_t number_of_points_ , size_t dimension )
+{
+	//standard file with barcode
+    std::vector< std::pair< double , double > > p = read_standard_file( filename );    
+    //gudhi file with barcode
+    //std::vector< std::pair< double , double > > p = read_gudhi_file( filename , dimension );        
+	
+	this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ );
+}
+
+Persistence_landscape_on_grid::Persistence_landscape_on_grid(const char* filename , size_t number_of_points_ )
+{
+	//standard file with barcode
+    std::vector< std::pair< double , double > > p = read_standard_file( filename );    
+    //gudhi file with barcode
+    //std::vector< std::pair< double , double > > p = read_gudhi_file( filename , dimension );     
+    
+    double grid_min_ = std::numeric_limits<double>::max();
+    double grid_max_ = -std::numeric_limits<double>::max();
+    for ( size_t i = 0 ; i != p.size() ; ++i )
+    {
+		if ( p[i].first < grid_min_ )grid_min_ = p[i].first;
+		if ( p[i].second > grid_max_ )grid_max_ = p[i].second;
+	}	
+	this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ );
+}
+
+void Persistence_landscape_on_grid::load_landscape_from_file( const char* filename )
+{
+	//check if the file exist.
+	if ( !( access( filename, F_OK ) != -1 ) )
+	{
+		cerr << "The file : " << filename << " do not exist. The program will now terminate \n";
+		throw "The file from which you are trying to read the persistence landscape do not exist. The program will now terminate \n";
+	}	
+	std::ifstream in;
+	in.open( filename );
+	
+	size_t number_of_points_in_the_grid = 0;
+	in >> this->grid_min >> this->grid_max >> number_of_points_in_the_grid;
+	
+	std::vector< std::vector< double > > v(number_of_points_in_the_grid);
+	std::string line;
+	std::getline(in, line);
+	double number;
+	for ( size_t i = 0 ; i != number_of_points_in_the_grid ; ++i )
+	{
+		//read a line of a file and convert it to a vector.
+		std::vector< double > vv;		
+		std::getline(in, line);		
+		//cerr << "Reading line : " << line << endl;getchar();		
+		std::istringstream stream(line);
+		while (stream >> number)
+		{
+			vv.push_back(number);
+		}		
+		v[i] = vv;
+	}
+	this->values_of_landscapes = v;
+	in.close();
+}
+
+void Persistence_landscape_on_grid::print_to_file( const char* filename )const
+{
+	std::ofstream out;
+	out.open( filename );
+	
+	//first we store the parameters of the grid:
+	out << grid_min << std::endl << grid_max << std::endl << this->values_of_landscapes.size() << std::endl;
+	
+	//and now in the following lines, the values of this->values_of_landscapes for the following arguments:
+	for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+	{
+		for ( size_t j = 0 ; j != this->values_of_landscapes[i].size() ; ++j )
+		{
+			out << this->values_of_landscapes[i][j] << " ";
+		}
+		out << std::endl;
+	}
+	
+	out.close();
+}
+
+void Persistence_landscape_on_grid::plot( const char* filename , size_t from_ , size_t to_ )
+{		
+    //this program create a gnuplot script file that allows to plot persistence diagram.
+    ofstream out;
+
+    std::ostringstream nameSS;
+    nameSS << filename << "_GnuplotScript";
+    std::string nameStr = nameSS.str();
+    out.open( (char*)nameStr.c_str() );
+
+    std::pair<double,double> min_max = compute_minimum_maximum();
+    out << "set xrange [" << this->grid_min << " : " << this->grid_max << "]" << endl;
+    out << "set yrange [" << min_max.first << " : " << min_max.second << "]" << endl;
+    
+    size_t number_of_nonzero_levels = this->number_of_nonzero_levels();
+    double dx = ( this->grid_max - this->grid_min )/((double)this->values_of_landscapes.size()-1);
+    
+    
+    size_t from = 0;
+    if ( from_ != std::numeric_limits<size_t>::max() )
+    {
+		if ( from_ < number_of_nonzero_levels )
+		{
+			from = from_;
+		}
+		else
+		{
+			return;
+		}
+	}	
+	size_t to = number_of_nonzero_levels;
+	if ( to_ != std::numeric_limits<size_t>::max() )
+    {
+		if ( to_ < number_of_nonzero_levels )
+		{
+			to = to_;
+		}
+	}
+    
+    
+    out << "plot ";    
+    for ( size_t lambda= from ; lambda != to ; ++lambda )
+    {
+        //out << "     '-' using 1:2 title 'l" << lambda << "' with lp";
+        out << "     '-' using 1:2 notitle with lp";
+        if ( lambda+1 != to )
+        {
+            out << ", \\";
+        }
+        out << endl;
+    }
+    
+    for ( size_t lambda = from ; lambda != to ; ++lambda )
+    {
+		double point = this->grid_min;
+        for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+        {
+			double value = 0;
+			if ( this->values_of_landscapes[i].size() > lambda )			
+			{
+				value = this->values_of_landscapes[i][lambda];
+			}
+            out << point << " " << value << endl;
+            point += dx;
+        }
+        out << "EOF" << endl;
+    }
+    cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << endl;
+}
+
+template < typename T > 
+Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid ( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 )
+{
+	//first we need to check if the domains are the same:
+	if ( !check_if_defined_on_the_same_domain(land1,land2) )throw "Two grids are not compatible"; 	
+	
+	T oper;
+	Persistence_landscape_on_grid result;
+	result.values_of_landscapes = std::vector< std::vector< double > >( land1.values_of_landscapes.size()  );
+	result.grid_min = land1.grid_min;
+	result.grid_max = land1.grid_max;
+	
+	//now we perorm the operations:
+	for ( size_t grid_point = 0 ; grid_point != land1.values_of_landscapes.size() ; ++grid_point )
+	{
+		result.values_of_landscapes[grid_point] = std::vector< double >( std::max( land1.values_of_landscapes[grid_point].size() , land2.values_of_landscapes[grid_point].size() ) );
+		for ( size_t lambda = 0 ; lambda != std::max( land1.values_of_landscapes[grid_point].size() , land2.values_of_landscapes[grid_point].size() ) ; ++lambda )
+		{
+			double value1 = 0;
+			double value2 = 0;
+			if ( lambda < land1.values_of_landscapes[grid_point].size() )value1 = land1.values_of_landscapes[grid_point][lambda];
+			if ( lambda < land2.values_of_landscapes[grid_point].size() )value2 = land2.values_of_landscapes[grid_point][lambda];
+			result.values_of_landscapes[grid_point][lambda] = oper( value1 , value2 );
+		}
+	}
+	
+	return result;
+}
+
+Persistence_landscape_on_grid Persistence_landscape_on_grid::multiply_lanscape_by_real_number_not_overwrite( double x )const
+{
+	Persistence_landscape_on_grid result;
+	result.values_of_landscapes = std::vector< std::vector< double > >( this->values_of_landscapes.size()  );
+	result.grid_min = this->grid_min;
+	result.grid_max = this->grid_max;
+	
+	for ( size_t grid_point = 0 ; grid_point != this->values_of_landscapes.size() ; ++grid_point )
+	{
+		result.values_of_landscapes[grid_point] = std::vector< double >( this->values_of_landscapes[grid_point].size() );
+		for ( size_t i = 0 ; i != this->values_of_landscapes[grid_point].size() ; ++i )
+		{
+			result.values_of_landscapes[grid_point][i] = x*this->values_of_landscapes[grid_point][i];
+		}
+	}
+	
+	return result;
+}
+
+double compute_max_norm_discance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second )
+{
+	double result = 0;
+	
+	//first we need to check if first and second is defined on the same domain"
+	if ( !check_if_defined_on_the_same_domain(first, second) )throw "Two grids are not compatible"; 
+	
+	for ( size_t i = 0 ; i != first.values_of_landscapes.size() ; ++i )
+	{
+		for ( size_t j = 0 ; j != std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) ; ++j )
+		{
+			if ( result < abs( first.values_of_landscapes[i][j] - second.values_of_landscapes[i][j] ) )
+			{
+				result = abs( first.values_of_landscapes[i][j] - second.values_of_landscapes[i][j] );
+			}
+		}
+		if ( first.values_of_landscapes[i].size() == std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) )
+		{
+			for ( size_t j = first.values_of_landscapes[i].size() ; j != second.values_of_landscapes[i].size() ; ++j )
+			{
+				if ( result < second.values_of_landscapes[i][j] )result = second.values_of_landscapes[i][j];
+			}
+		}
+		if ( second.values_of_landscapes[i].size() == std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) )
+		{
+			for ( size_t j = second.values_of_landscapes[i].size() ; j != first.values_of_landscapes[i].size() ; ++j )
+			{
+				if ( result < first.values_of_landscapes[i][j] )result = first.values_of_landscapes[i][j];
+			}
+		}
+	}  
+	return result;
+}
+
+
+
+}//namespace Gudhi_stat
+}//namespace Gudhi
+
+#endif
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h b/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h
new file mode 100644
index 00000000..d2e6ec3c
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h
@@ -0,0 +1,136 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+#ifndef Read_Persitence_From_File_H
+#define Read_Persitence_From_File_H
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <unistd.h>
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+
+std::vector< std::pair< double , double > > read_standard_file( const char* filename )
+{
+	bool dbg = false;
+	
+	std::ifstream in;
+    in.open( filename );
+    if ( !( access( filename, F_OK ) != -1 ) )
+	{
+		std::cerr << "The file : " << filename << " do not exist. The program will now terminate \n";
+		throw "The file from which you are trying to read the persistence landscape do not exist. The program will now terminate \n";
+	}
+
+    std::string line;
+    std::vector< std::pair<double,double> > barcode;
+
+    while (!in.eof())
+    {
+        getline(in,line);
+        if ( !(line.length() == 0 || line[0] == '#') )
+        {
+            std::stringstream lineSS;
+            lineSS << line;
+            double beginn, endd;
+            lineSS >> beginn;
+            lineSS >> endd;
+            if ( beginn > endd )
+            {
+                double b = beginn;
+                beginn = endd;
+                endd = b;
+            }
+            barcode.push_back( std::make_pair( beginn , endd ) );
+            if (dbg)
+            {
+				std::cerr << beginn << " , " << endd << std::endl;
+			}
+        }
+	}
+	in.close();
+	return barcode;
+}//read_standard_file
+
+std::vector< std::pair< double , double > > read_gudhi_file( const char* filename , size_t dimension = 0 )
+{
+	bool dbg = false;
+	std::ifstream in;
+    in.open( filename );
+
+    std::string line;
+    std::vector< std::pair<double,double> > barcode;
+
+    while (!in.eof())
+    {
+        getline(in,line);
+        if ( !(line.length() == 0 || line[0] == '#') )
+        {
+			if ( line.find("inf") != std::string::npos )
+			{
+				if ( dbg )
+				{
+					std::cerr << "This line: " << line << " contains infinite interval. We will skip it. \n";
+				}
+				continue;
+			}
+            std::stringstream lineSS;
+            lineSS << line;
+            double beginn, endd, field, dim;
+            lineSS >> field;
+            lineSS >> dim;
+            lineSS >> beginn;
+            lineSS >> endd;
+            if ( beginn > endd )
+            {
+                double b = beginn;
+                beginn = endd;
+                endd = b;
+            }
+            if ( dim == dimension )
+            { 
+				barcode.push_back( std::make_pair( beginn , endd ) );
+				if (dbg)
+				{
+					std::cerr << beginn << " , " << endd << std::endl;
+				}
+			}
+        }
+	}
+	in.close();
+	return barcode;
+}//read_gudhi_file
+
+}//namespace Gudhi_stat
+}//namespace Gudhi
+
+
+
+
+#endif