adding persistence_representations subfolder (at the moment with all the junk).

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

Former-commit-id: fdf52ca688d992655ff139f42e316d5b7b3f8a12

10 files changed, 4555 insertions, 0 deletions

diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/PSSK.h b/src/Gudhi_stat/include/gudhi/persistence_representations/PSSK.h
new file mode 100644
index 00000000..22513097
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/PSSK.h
@@ -0,0 +1,173 @@
+/*    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 (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/>.
+ */
+
+#pragma once
+#ifndef PSSK_H
+#define PSSK_H
+
+//gudhi include
+#include <gudhi/persistence_representations/Persistence_heat_maps.h>
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+/**
+* This is a version of a representation presented in https://arxiv.org/abs/1412.6821
+* In that paper the authors are using the representation just to compute kernel. Over here, we extend the usability by far.
+* Note that the version presented here is not exact, since we are discretizing the kernel.
+* The only difference with respect to the original class is the method of creation. We have full (square) image, and for every point (p,q), we add a kernel at (p,q) and the negative kernel
+* at (q,p)
+**/
+
+class PSSK : public Persistence_heat_maps<constant_scaling_function>
+{
+public:
+	PSSK():Persistence_heat_maps(){}
+	
+	PSSK(const std::vector< std::pair< double,double > >  & interval , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1 )
+	:Persistence_heat_maps()
+	{
+		this->construct( interval  , filter , number_of_pixels , min_ , max_ );
+	}
+	
+	 
+	 PSSK( const char* filename , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  ):
+	    Persistence_heat_maps()
+	    {
+			std::vector< std::pair< double , double > > intervals_ = read_standard_file( filename );         
+			this->construct( intervals_ , filter , number_of_pixels , min_ , max_ );
+		}
+
+protected:
+	void construct( const std::vector< std::pair<double,double> >& intervals_  , 
+					std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1),
+                    size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1 );
+};
+
+//if min_ == max_, then the program is requested to set up the values itself based on persistence intervals
+void PSSK::construct( const std::vector< std::pair<double,double> >& intervals_  ,  
+									   std::vector< std::vector<double> > filter,
+									   size_t number_of_pixels , double min_ , double max_ )
+{	
+    bool dbg = false;       
+    if ( dbg ){std::cerr << "Entering construct procedure \n";getchar();}
+
+    if ( min_ == max_ )
+    {
+        //in this case, we want the program to set up the min_ and max_ values by itself.
+        min_ = std::numeric_limits<int>::max();
+        max_ = -std::numeric_limits<int>::max();
+        
+        
+        for ( size_t i = 0 ; i != intervals_.size() ; ++i )
+        {
+			if ( intervals_[i].first < min_ )min_ = intervals_[i].first;
+			if ( intervals_[i].second > max_ )max_ = intervals_[i].second;
+		}
+        //now we have the structure filled in, and moreover we know min_ and max_ values of the interval, so we know the range.
+
+        //add some more space:
+        min_ -= fabs(max_ - min_)/100;
+        max_ += fabs(max_ - min_)/100;
+    }
+
+	if ( dbg )
+	{
+		std::cerr << "min_ : " << min_ << std::endl;
+		std::cerr << "max_ : " << max_ << std::endl;
+		std::cerr << "number_of_pixels : " << number_of_pixels << std::endl;
+		getchar();
+	}
+
+    this->min_ = min_;
+    this->max_ = max_;    
+
+
+
+    //initialization of the structure heat_map
+    std::vector< std::vector<double> > heat_map_;
+    for ( size_t i = 0 ; i != number_of_pixels ; ++i )
+    {       
+		std::vector<double> v( number_of_pixels , 0 );
+        heat_map_.push_back( v );
+    }
+    this->heat_map = heat_map_;
+
+    if (dbg)std::cerr << "Done creating of the heat map, now we will fill in the structure \n";
+
+	for ( size_t pt_nr = 0 ; pt_nr != intervals_.size() ; ++pt_nr )
+	{
+		//compute the value of intervals_[pt_nr] in the grid:
+		int x_grid = (int)((intervals_[pt_nr].first - this->min_)/( this->max_-this->min_ )*number_of_pixels);
+		int y_grid = (int)((intervals_[pt_nr].second - this->min_)/( this->max_-this->min_ )*number_of_pixels);
+		
+		if ( dbg )
+		{
+			std::cerr << "point : " << intervals_[pt_nr].first << " , " << intervals_[pt_nr].second << std::endl;
+			std::cerr << "x_grid : " << x_grid << std::endl;
+			std::cerr << "y_grid : " << y_grid << std::endl;
+		}
+		
+		//x_grid and y_grid gives a center of the kernel. We want to have its lower left cordner. To get this, we need to shift x_grid and y_grid by a grid diameter.		
+		x_grid -= filter.size()/2;
+		y_grid -= filter.size()/2;
+		//note that the numbers x_grid and y_grid may be negative. 
+		
+		if ( dbg )
+		{
+			std::cerr << "After shift : \n";;
+			std::cerr << "x_grid : " << x_grid << std::endl;
+			std::cerr << "y_grid : " << y_grid << std::endl;
+			std::cerr << "filter.size() : " << filter.size() << std::endl;
+			getchar();
+		}					
+
+		
+		for ( size_t i = 0 ; i != filter.size() ; ++i )
+		{
+			for ( size_t j = 0 ; j != filter.size() ; ++j )
+			{			
+				//if the point (x_grid+i,y_grid+j) is the correct point in the grid.						
+				if ( 
+					  ((x_grid+i)>=0) && (x_grid+i<this->heat_map.size()) 
+					  &&
+					  ((y_grid+j)>=0) && (y_grid+j<this->heat_map.size()) 
+				   )
+				{
+					if ( dbg ){std::cerr << y_grid+j << " " <<  x_grid+i << std::endl;}
+					this->heat_map[ y_grid+j ][ x_grid+i ] += filter[i][j];
+					this->heat_map[ x_grid+i ][ y_grid+j ] += -filter[i][j];
+				}
+			}
+		}
+		
+	}
+}//construct
+
+
+#endif
+
+}//namespace Gudhi_stat
+}//namespace Gudhi
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_heat_maps.h b/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_heat_maps.h
new file mode 100644
index 00000000..dd95fbec
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_heat_maps.h
@@ -0,0 +1,1003 @@
+/*    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 (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_heat_maps_H
+#define Persistence_heat_maps_H
+ 
+//standard include
+#include <vector>
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <limits>
+#include <vector>
+#include <algorithm>
+
+//gudhi include
+#include <gudhi/read_persitence_from_file.h>
+#include <gudhi/common_gudhi_stat.h>
+
+
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+
+/**
+ * This is a simple procedure to create n by n (or 2*pixel_radius times 2*pixel_radius cubical approximation of a Gaussian kernel.
+**/ 
+std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius , double sigma )
+{
+	bool dbg = false;
+	//we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a size 2*pixel_radius times 2*pixel_radius.
+	
+    double r = 0;
+    double sigma_sqr = sigma * sigma;
+ 
+    // sum is for normalization
+    double sum = 0;
+    
+    //initialization of a kernel:
+    std::vector< std::vector<double> > kernel( 2*pixel_radius +1 );
+    for ( size_t i = 0 ; i != kernel.size() ; ++i )
+    {
+		std::vector<double> v( 2*pixel_radius +1 , 0 );
+		kernel[i] = v;
+	}
+	
+	if ( dbg )
+	{
+		std::cerr << "Kernel initalize \n";	
+		std::cerr << "pixel_radius : " << pixel_radius << std::endl; 
+		std::cerr << "kernel.size() : " << kernel.size() << std::endl;
+		getchar();
+	}
+ 
+    for (int x = -pixel_radius; x <= (int)pixel_radius; x++)
+    {
+        for(int y = -pixel_radius; y <= (int)pixel_radius; y++)
+        {			
+			double real_x = 2*sigma*x/pixel_radius;
+			double real_y = 2*sigma*y/pixel_radius;						
+            r = sqrt(real_x*real_x + real_y*real_y);						
+            kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r*r)/sigma_sqr))/(3.141592 * sigma_sqr);
+            sum += kernel[x + pixel_radius][y + pixel_radius];            
+        }
+    }      
+ 
+    // normalize the kernel
+    for( size_t i = 0; i != kernel.size() ; ++i)
+    {
+        for( size_t j = 0; j != kernel[i].size() ; ++j)
+        {
+            kernel[i][j] /= sum;
+		}
+            
+    }
+    
+    if ( dbg )
+    {
+		std::cerr << "Here is the kernel : \n";
+		for( size_t i = 0; i != kernel.size() ; ++i)
+		{
+			for( size_t j = 0; j != kernel[i].size() ; ++j)
+			{
+				std::cerr << kernel[i][j] << " ";
+			}
+			std::cerr << std::endl;
+		}
+	}
+    return kernel;
+}
+
+
+/**
+ * There are various options to scale the poits depending on their location. One can for instance:
+ * (1) do nothing (scale all of them with the weight 1), as in the function constant_function
+ * (2) Scale them by the distance to the diagonal. This is implemented in function
+ * (3) Scale them with the square of their distance to diagonal. This is implemented in function
+ * (4) Scale them with 
+**/ 
+
+class constant_scaling_function
+{
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{	
+		return 1;
+	}
+};
+
+class distance_from_diagonal_scaling
+{
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		//(point_in_diagram.first+point_in_diagram.second)/2.0
+		return sqrt( pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2)  );
+	}
+};
+
+class squared_distance_from_diagonal_scaling
+{
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		return pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2);
+	}
+};
+
+class arc_tan_of_persistence_of_point
+{
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		return atan( point_in_diagram.second - point_in_diagram.first );
+	}
+};
+
+class weight_by_setting_maximal_interval_to_have_length_one
+{
+public:
+	weight_by_setting_maximal_interval_to_have_length_one( double len ):letngth_of_maximal_interval(len){}
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		return (point_in_diagram.second-point_in_diagram.first)/this->letngth_of_maximal_interval;
+	}
+private:
+	double letngth_of_maximal_interval;
+};
+
+
+/**
+ * This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product
+**/ 
+template <typename Scalling_of_kernels = constant_scaling_function>
+class Persistence_heat_maps
+{
+public:
+    /**
+     * The default constructor. A scaling function from the diagonal is set up to a constant function. The image is not erased below the diagonal. The gaussian have diameter 5. 	 
+    **/ 
+	Persistence_heat_maps()
+	{
+		Scalling_of_kernels f;
+		this->f = f;
+		this->erase_below_diagonal = false;		
+		this->min_ = this->max_ = 0;
+		this->set_up_parameters_for_basic_classes();
+	};
+	
+	/**
+	 * Construction that takes at the input the following parameters:
+	 * (1) A vector of pairs of doubles (representing persistence intervals). All other parameters are optional. They are:
+	 * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), 
+	 * (3) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). 
+	 * (4) a number of pixels in each direction (set to 1000 by default). 
+	 * (5) a min x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data,
+	 * (6) a max x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data.
+	**/  
+    Persistence_heat_maps( const std::vector< std::pair< double,double > >  & interval , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = std::numeric_limits<double>::max() , double max_ = std::numeric_limits<double>::max()  );
+    
+    /**
+	 * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by default), a scaling function (constant by default), a boolean value which determines if the area of image below diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each direction (set to 1000 by default). and min and max values of images (both set to std::numeric_limits<double>::max() by defaulet. If this is the case, the program will pick the right values based on the data).
+	**/  
+	/**
+	 * Construction that takes at the input the following parameters:
+	 * (1) A a name of a file with persistence intervals. The file shold be readable by the function read_standard_file. All other parameters are optional. They are:
+	 * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), 	
+	 * (3) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). 
+	 * (4) a number of pixels in each direction (set to 1000 by default). 
+	 * (5) a min x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data,
+	 * (6) a max x and y value of points that are to be taken into account. By default it is set to std::numeric_limits<double>::max(), in which case the program compute the values based on the data.
+	**/ 	
+    Persistence_heat_maps( const char* filename , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1), bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = std::numeric_limits<double>::max() , double max_ = std::numeric_limits<double>::max()  );
+
+	
+	/**
+	 * Compute a mean value of a collection of heat maps and store it in the current object. Note that all the persistence maps send in a vector to this procedure need to have the same parameters. 
+	 * If this is not the case, the program will throw an exception. 
+	**/ 
+	void compute_mean( const std::vector<Persistence_heat_maps*>& maps );
+	
+	/**
+	 * Compute a median value of a collection of heat maps and store it in the current object. Note that all the persistence maps send in a vector to this procedure need to have the same parameters. 
+	 * If this is not the case, the program will throw an exception. 
+	**/ 
+	void compute_median( const std::vector<Persistence_heat_maps*>& maps );	
+	
+	/**
+	 * Compute a percentage of active (i.e) values above the cutoff of a collection of heat maps. 
+	**/ 
+	void compute_percentage_of_active( const std::vector<Persistence_heat_maps*>& maps , size_t cutoff = 1 );
+    
+    //put to file subroutine
+    /**
+     * The function outputs the perssitence image to a text file. The format as follow:
+     * In the first line, the values min and max of the image are stored
+     * In the next lines, we have the persistence images in a form of a bitmap image. 
+    **/
+    void print_to_file( const char* filename )const;
+    
+    /**
+     * A function that load a heat map from file to the current object (and arase qhatever was stored in the current object before).
+    **/
+    void load_from_file( const char* filename );
+    
+    
+    /**
+     * The procedure checks if min_, max_ and this->heat_maps sizes are the same. 
+    **/ 
+    inline bool check_if_the_same( const Persistence_heat_maps& second  )const
+    {
+		bool dbg = false;
+		if ( this->heat_map.size() != second.heat_map.size() )
+		{
+			if ( dbg )std::cerr << "this->heat_map.size() : " << this->heat_map.size()  << " \n second.heat_map.size() : " << second.heat_map.size() << std::endl;
+			return false;
+		}
+		if ( this->min_ != second.min_ )
+		{
+			if ( dbg )std::cerr << "this->min_ : " << this->min_ << ", second.min_ : " << second.min_ << std::endl;
+			return false;
+		}
+		if ( this->max_ != second.max_ )
+		{
+			if ( dbg )std::cerr << "this->max_ : " << this->max_ << ", second.max_ : " << second.max_ << std::endl;
+			return false;
+		}
+		//in the other case we may assume that the persistence images are defined on the same domain.
+		return true;
+	}
+	
+	
+	/**
+	 * Return minimal range value of persistent image.
+	**/ 
+	inline double get_min()const{return this->min_;}
+
+	/**
+	 * Return maximal range value of persistent image.
+	**/ 
+	inline double get_max()const{return this->max_;}
+	
+	
+	/**
+	 * Operator == to check if to persistence heat maps are the same.
+	**/ 
+	bool operator == ( const Persistence_heat_maps& rhs )const
+	{
+		bool dbg = false;
+		if ( !this->check_if_the_same(rhs) )
+		{
+			if ( dbg )std::cerr << "The domains are not the same \n";
+			return false;//in this case, the domains are not the same, so the maps cannot be the same.
+		}
+		for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+		{
+			for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+			{
+				if ( !almost_equal(this->heat_map[i][j] , rhs.heat_map[i][j]) )
+				{
+					if ( dbg )
+					{
+						std::cerr << "this->heat_map[" << i << "][" << j << "] = " << this->heat_map[i][j] << std::endl;
+						std::cerr << "rhs.heat_map[" << i << "][" << j << "] = " << rhs.heat_map[i][j] << std::endl;						
+					}
+					return false;
+				}
+			}
+		}
+		return true;
+	}
+	
+	/**
+	 * Operator != to check if to persistence heat maps are different.
+	**/ 
+	bool operator != ( const Persistence_heat_maps& rhs )const
+	{
+		return !( (*this) == rhs );
+	}
+    
+    
+    /**
+     * A function to generate a gnuplot script to vizualize the persistent image.  
+    **/
+    void plot( const char* filename )const;
+    
+    
+    /**
+     * Binary (arythmetic) operation on two Persistence_heat_maps.
+    **/     
+    template<typename Operation_type>
+    friend Persistence_heat_maps operation_on_pair_of_heat_maps( const Persistence_heat_maps& first ,  const Persistence_heat_maps& second , Operation_type operation )
+    {
+		//first check if the heat maps are compatible 
+		if ( !first.check_if_the_same( second ) )
+		{
+			std::cerr << "Sizes of the heat maps are not compatible. The program will now terminate \n";
+			throw "Sizes of the heat maps are not compatible. The program will now terminate \n";
+		}
+		Persistence_heat_maps result;
+		result.min_ = first.min_;
+		result.max_ = first.max_;
+		result.heat_map.reserve( first.heat_map.size() );
+		for ( size_t i = 0 ; i != first.heat_map.size() ; ++i )
+		{ 
+			std::vector< double > v;
+			v.reserve( first.heat_map[i].size() );
+			for ( size_t j = 0 ; j != first.heat_map[i].size() ; ++j )
+			{
+				v.push_back( operation( first.heat_map[i][j] , second.heat_map[i][j] ) );
+			}
+			result.heat_map.push_back( v );
+		}
+		return result;
+	}//operation_on_pair_of_heat_maps
+	
+	
+	/**
+	 * Multiplication of Persistence_heat_maps by scalar (so that all values of the heat map gets multiplied by that scalar). 
+	**/ 
+	Persistence_heat_maps multiply_by_scalar( double scalar )const 
+    {
+		Persistence_heat_maps result;
+		result.min_ = this->min_;
+		result.max_ = this->max_;
+		result.heat_map.reserve( this->heat_map.size() );
+		for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+		{ 
+			std::vector< double > v;
+			v.reserve( this->heat_map[i].size() );
+			for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+			{
+				v.push_back( this->heat_map[i][j] * scalar );
+			}
+			result.heat_map.push_back( v );
+		}
+		return result;
+	}
+	
+    /**
+     * This function computes a sum of two objects of a type Persistence_heat_maps.
+    **/    
+    friend Persistence_heat_maps operator+( const Persistence_heat_maps& first , const Persistence_heat_maps& second )
+    {
+		return operation_on_pair_of_heat_maps( first , second , std::plus<double>() );
+	}	
+	/**
+     * This function computes a difference of two objects of a type Persistence_heat_maps.
+    **/
+    friend Persistence_heat_maps operator-( const Persistence_heat_maps& first , const Persistence_heat_maps& second )
+    {
+		return operation_on_pair_of_heat_maps( first , second , std::minus<double>() );
+	}
+	/**
+     * This function computes a product of an object of a type Persistence_heat_maps with real number.
+    **/
+	friend Persistence_heat_maps operator*( double scalar , const Persistence_heat_maps& A )
+	{
+		return A.multiply_by_scalar( scalar );
+	}
+	/**
+     * This function computes a product of an object of a type Persistence_heat_maps with real number.
+    **/
+	friend Persistence_heat_maps operator*( const Persistence_heat_maps& A , double scalar )
+	{
+		return A.multiply_by_scalar( scalar );
+	}
+	/**
+     * This function computes a product of an object of a type Persistence_heat_maps with real number.
+    **/
+	Persistence_heat_maps operator*( double scalar )
+	{
+		return this->multiply_by_scalar( scalar );
+	}
+	/**
+	 * += operator for Persistence_heat_maps.
+	**/ 
+    Persistence_heat_maps operator += ( const Persistence_heat_maps& rhs )
+    {
+        *this = *this + rhs;
+        return *this;
+    }    
+    /**
+	 * -= operator for Persistence_heat_maps.
+	**/ 
+    Persistence_heat_maps operator -= ( const Persistence_heat_maps& rhs )
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+    /**
+	 * *= operator for Persistence_heat_maps.
+	**/ 
+    Persistence_heat_maps operator *= ( double x )
+    {
+        *this = *this*x;
+        return *this;
+    }
+    /**
+	 * /= operator for Persistence_heat_maps.
+	**/ 
+    Persistence_heat_maps operator /= ( double x )
+    {
+        if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." );
+        *this = *this * (1/x);
+        return *this;
+    }
+    
+    
+    //Implementations of functions for various concepts.  
+    
+    /**
+	 * This function produce a vector of doubles based on a persisence heat map. It is required in a concept Vectorized_topological_data
+	*/
+    std::vector<double> vectorize( int number_of_function )const;    
+    /**
+	 * This function return the number of functions that allows vectorization of persistence heat map. It is required in a concept Vectorized_topological_data.
+	 **/ 
+	size_t number_of_vectorize_functions()const
+	{
+		return this->number_of_functions_for_vectorization;	
+	}
+	
+	/**	 
+	 * This function is required by the Real_valued_topological_data concept. It returns various projections od the persistence heat map to a real line.
+	**/ 	
+	double project_to_R( int number_of_function )const;	
+	/**
+	 * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept.
+	**/ 
+	size_t number_of_projections_to_R()const
+	{
+		return this->number_of_functions_for_projections_to_reals;
+	}
+	
+	/**
+	* A function to compute distance between persistence heat maps.
+	* The parameter of this function is a const reference to an object of a class Persistence_heat_maps.
+	* This function is required in Topological_data_with_distances concept.
+    * For max norm distance, set power to std::numeric_limits<double>::max()
+	**/
+	double distance( const Persistence_heat_maps& second_ , double power = 1)const;
+	
+	/**
+	 * A function to compute averaged persistence heat map, based on vector of persistence heat maps.
+	 * This function is required by Topological_data_with_averages concept.
+	**/ 	
+	void compute_average( const std::vector< Persistence_heat_maps* >& to_average );	
+	
+	/**
+	* A function to compute scalar product of persistence heat maps.
+	* The parameter of this functionis a const reference to an object of a class Persistence_heat_maps.
+	* This function is required in Topological_data_with_scalar_product concept.
+	**/
+	double compute_scalar_product( const Persistence_heat_maps& second_ )const;
+	
+	//end of implementation of functions needed for concepts.
+	
+	
+	/**
+	 * The x-range of the persistence heat map. 
+	**/ 
+	std::pair< double , double > get_x_range()const
+	{
+		return std::make_pair( this->min_ , this->max_ );
+	}
+	
+	/**
+	 * The y-range of the persistence heat map. 
+	**/ 
+	std::pair< double , double > get_y_range()const
+	{
+		return this->get_x_range();
+	}
+	
+	
+	
+	
+protected:
+	//private methods
+	std::vector< std::vector<double> > check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps );
+	size_t number_of_functions_for_vectorization;
+	size_t number_of_functions_for_projections_to_reals;
+    void construct( const std::vector< std::pair<double,double> >& intervals_  , 
+					std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1),
+					                    
+                    bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = std::numeric_limits<double>::max() , double max_ = std::numeric_limits<double>::max() );
+                    
+	void set_up_parameters_for_basic_classes()
+	{
+		this->number_of_functions_for_vectorization = 1;
+		this->number_of_functions_for_projections_to_reals = 1;
+	}
+    
+    //data    
+    //double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram );    
+    Scalling_of_kernels f;
+    bool erase_below_diagonal;
+    double min_;
+    double max_;
+    std::vector< std::vector< double > > heat_map;    
+};
+
+
+//if min_ == max_, then the program is requested to set up the values itself based on persistence intervals
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::construct( const std::vector< std::pair<double,double> >& intervals_  ,  
+									   std::vector< std::vector<double> > filter,									   
+									   bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ )
+{
+    bool dbg = false;       
+    if ( dbg )std::cerr << "Entering construct procedure \n";
+    Scalling_of_kernels f;
+    this->f = f;
+
+    if ( min_ == max_ )
+    {
+        //in this case, we want the program to set up the min_ and max_ values by itself.
+        min_ = std::numeric_limits<int>::max();
+        max_ = -std::numeric_limits<int>::max();
+        
+        
+        for ( size_t i = 0 ; i != intervals_.size() ; ++i )
+        {
+			if ( intervals_[i].first < min_ )min_ = intervals_[i].first;
+			if ( intervals_[i].second > max_ )max_ = intervals_[i].second;
+		}
+        //now we have the structure filled in, and moreover we know min_ and max_ values of the interval, so we know the range.
+
+        //add some more space:
+        min_ -= fabs(max_ - min_)/100;
+        max_ += fabs(max_ - min_)/100;
+    }
+
+	if ( dbg )
+	{
+		std::cerr << "min_ : " << min_ << std::endl;
+		std::cerr << "max_ : " << max_ << std::endl;
+		std::cerr << "number_of_pixels : " << number_of_pixels << std::endl;
+		getchar();
+	}
+
+    this->min_ = min_;
+    this->max_ = max_;        
+
+    //initialization of the structure heat_map
+    std::vector< std::vector<double> > heat_map_;
+    for ( size_t i = 0 ; i != number_of_pixels ; ++i )
+    {       
+		std::vector<double> v( number_of_pixels , 0 );
+        heat_map_.push_back( v );
+    }
+    this->heat_map = heat_map_;
+
+    if (dbg)std::cerr << "Done creating of the heat map, now we will fill in the structure \n";
+
+	for ( size_t pt_nr = 0 ; pt_nr != intervals_.size() ; ++pt_nr )
+	{
+		//compute the value of intervals_[pt_nr] in the grid:
+		int x_grid = (int)((intervals_[pt_nr].first - this->min_)/( this->max_-this->min_ )*number_of_pixels);
+		int y_grid = (int)((intervals_[pt_nr].second - this->min_)/( this->max_-this->min_ )*number_of_pixels);
+		
+		if ( dbg )
+		{
+			std::cerr << "point : " << intervals_[pt_nr].first << " , " << intervals_[pt_nr].second << std::endl;
+			std::cerr << "x_grid : " << x_grid << std::endl;
+			std::cerr << "y_grid : " << y_grid << std::endl;
+		}
+		
+		//x_grid and y_grid gives a center of the kernel. We want to have its lower left cordner. To get this, we need to shift x_grid and y_grid by a grid diameter.		
+		x_grid -= filter.size()/2;
+		y_grid -= filter.size()/2;
+		//note that the numbers x_grid and y_grid may be negative. 
+		
+		if ( dbg )
+		{
+			std::cerr << "After shift : \n";;
+			std::cerr << "x_grid : " << x_grid << std::endl;
+			std::cerr << "y_grid : " << y_grid << std::endl;
+		}
+				
+		double scaling_value = this->f(intervals_[pt_nr]);
+
+		
+		for ( size_t i = 0 ; i != filter.size() ; ++i )
+		{
+			for ( size_t j = 0 ; j != filter.size() ; ++j )
+			{			
+				//if the point (x_grid+i,y_grid+j) is the correct point in the grid.						
+				if ( 
+					  ((x_grid+i)>=0) && (x_grid+i<this->heat_map.size()) 
+					  &&
+					  ((y_grid+j)>=0) && (y_grid+j<this->heat_map.size()) 
+				   )
+				{
+					if ( dbg ){std::cerr << y_grid+j << " " <<  x_grid+i << std::endl;}
+					this->heat_map[ y_grid+j ][ x_grid+i ] += scaling_value * filter[i][j];
+					if ( dbg )
+					{
+						std::cerr << "Position : (" << x_grid+i << "," << y_grid+j << ") got increased by the value : " << filter[i][j] << std::endl;
+					}					
+				}
+			}
+		}
+		
+	}
+	
+	//now it remains to cut everything below diagonal if the user wants us to. 
+	if ( erase_below_diagonal )
+	{
+		for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+		{
+			for ( size_t j = i ; j != this->heat_map.size() ; ++j )
+			{
+				this->heat_map[i][j] = 0;
+			}
+		}
+	}
+}//construct
+
+template <typename Scalling_of_kernels>
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( const std::vector< std::pair< double,double > >  & interval , 
+											  std::vector< std::vector<double> > filter,  
+											  bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ )
+{
+    this->construct( interval ,  filter , erase_below_diagonal , number_of_pixels , min_ , max_ );
+    this->set_up_parameters_for_basic_classes();
+}
+
+
+template <typename Scalling_of_kernels>
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( const char* filename , 
+											  std::vector< std::vector<double> > filter, 											  
+											  bool erase_below_diagonal , size_t number_of_pixels  , double min_ , double max_ )
+{    
+    std::vector< std::pair< double , double > > intervals_ = read_standard_file( filename );           
+    this->construct( intervals_ ,  filter, erase_below_diagonal , number_of_pixels , min_ , max_ );
+    this->set_up_parameters_for_basic_classes();
+}
+
+
+template <typename Scalling_of_kernels>
+std::vector< std::vector<double> > Persistence_heat_maps<Scalling_of_kernels>::check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps )
+{
+	//checking if all the heat maps are of the same size:
+	for ( size_t i = 0 ; i != maps.size() ; ++i )
+    {
+		if ( maps[i]->heat_map.size() != maps[0]->heat_map.size() )
+		{
+			std::cerr << "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
+			throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
+		}
+		if ( maps[i]->heat_map[0].size() != maps[0]->heat_map[0].size() )
+		{
+			std::cerr << "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";
+			throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";			
+		}
+	} 	
+	std::vector< std::vector<double> > heat_maps( maps[0]->heat_map.size() );
+	for ( size_t i = 0 ; i != maps[0]->heat_map.size() ; ++i )
+    {
+		 std::vector<double> v( maps[0]->heat_map[0].size() , 0 );
+		 heat_maps[i] = v;
+	}
+	return heat_maps;
+}
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_median( const std::vector<Persistence_heat_maps*>& maps )
+{
+	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
+	
+	std::vector<double> to_compute_median( maps.size() );
+    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
+    {
+        for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j )
+        {
+			for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no )
+			{
+				to_compute_median[map_no] = maps[map_no]->heat_map[i][j];
+			}
+			std::nth_element(to_compute_median.begin(), to_compute_median.begin() + to_compute_median.size()/2, to_compute_median.end());
+            heat_maps[i][j] = to_compute_median[ to_compute_median.size()/2 ];
+        }       
+    }
+	this->heat_map = heat_maps;
+	this->min_= maps[0]->min_;
+	this->max_= maps[0]->max_;
+}
+
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_mean( const std::vector<Persistence_heat_maps*>& maps )
+{
+	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
+		    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
+    {
+        for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j )
+        {
+			double mean = 0;
+			for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no )
+			{	
+				mean += maps[map_no]->heat_map[i][j];
+			}		
+            heat_maps[i][j] =  mean/(double)maps.size();
+        }       
+    }
+	this->heat_map = heat_maps;
+	this->min_ = maps[0]->min_;
+	this->max_ = maps[0]->max_;
+}
+
+
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_percentage_of_active( const std::vector<Persistence_heat_maps*>& maps , size_t cutoff  )
+{
+	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
+
+    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
+    {
+        for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j )
+        {
+            size_t number_of_active_levels = 0;
+            for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no )
+            {
+                if ( maps[map_no]->heat_map[i][j] ) number_of_active_levels++;
+            }
+            if ( number_of_active_levels > cutoff )
+            {
+                heat_maps[i][j] = number_of_active_levels;
+            }
+            else
+            {
+                heat_maps[i][j] = 0;
+            }
+        }       
+    }
+    this->heat_map = heat_maps;
+    this->min_ = maps[0]->min_;
+	this->max_ = maps[0]->max_;
+}
+
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::plot( const char* filename )const
+{
+	std::ofstream out;
+	std::stringstream ss;
+	ss << filename << "_GnuplotScript";
+	
+    out.open( ss.str().c_str() );
+	out << "plot      '-' matrix with image" << std::endl;
+    for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+    {
+        for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+        {
+            out << this->heat_map[i][j] << " ";
+        }
+        out << std::endl;
+    }
+    out.close();
+    std::cout << "Gnuplot script have been created. Open gnuplot and type load \'" << ss.str().c_str() << "\' to see the picture." << std::endl; 
+}
+
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::print_to_file( const char* filename )const
+{
+
+	std::ofstream out;
+	out.open( filename );
+	
+	//First we store this->min_ and this->max_ values:
+	out << this->min_ << " " << this->max_ << std::endl;	
+	for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+    {
+        for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+        {
+            out << this->heat_map[i][j] << " ";
+        }
+        out << std::endl;
+    }
+	out.close();
+}
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::load_from_file( const char* filename )
+{
+	bool dbg = false;	
+	
+	std::ifstream in;
+	in.open( filename );
+	
+	//checking if the file exist / if it was open. 
+	if ( !( 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";
+	}
+	
+	//now we read the file one by one. 
+	
+	 
+	 
+	in >> this->min_ >> this->max_;	    
+	if ( dbg )
+	{
+		std::cerr << "Reading the following values of min and max : " << this->min_ << " , " << this->max_ << std::endl;
+	}
+
+	std::string temp;
+	std::getline(in, temp);
+         
+	while (!in.eof())
+	{
+		std::getline(in, temp);
+        std::stringstream lineSS;
+        lineSS << temp;
+        
+        std::vector<double> line_of_heat_map;    
+        while ( lineSS.good() )
+        {
+			double point;
+			
+			lineSS >> point;
+			line_of_heat_map.push_back( point );
+			if ( dbg )
+			{
+				std::cout << point << " ";
+			}
+		}
+		if ( dbg )
+		{
+			std::cout << std::endl;
+			getchar();
+		}
+
+		if ( in.good() )this->heat_map.push_back( line_of_heat_map );
+	}	
+	in.close();
+	if ( dbg )std::cout << "Done \n";
+}
+
+
+//Concretizations of virtual methods:   
+template <typename Scalling_of_kernels> 
+std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize( int number_of_function )const
+{
+	//convert this->heat_map into one large vector:
+	size_t size_of_result = 0;
+	for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+	{
+		size_of_result += this->heat_map[i].size();
+	}
+	
+	std::vector< double > result;
+	result.reserve( size_of_result );
+	
+	for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+	{
+		for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+		{
+			result.push_back( this->heat_map[i][j] );
+		}
+	}
+	
+	return result;
+}
+
+template <typename Scalling_of_kernels>
+double Persistence_heat_maps<Scalling_of_kernels>::distance( const Persistence_heat_maps& second , double power )const
+{			
+	//first we need to check if (*this) and second are defined on the same domain and have the same dimensions:
+	if ( !this->check_if_the_same(second) )
+	{
+		std::cerr << "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate";
+		throw "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate";
+	}
+	
+	//if we are here, we know that the two persistence iomages are defined on the same domain, so we can start computing their distances:
+	
+	double distance = 0;
+	if ( power < std::numeric_limits<double>::max() )
+	{
+		for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+		{
+			for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+			{
+				distance += pow( fabs(this->heat_map[i][j] - second.heat_map[i][j]) , power );
+			}
+		}
+	}
+	else
+	{		
+		//in this case, we compute max norm distance
+		for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+		{
+			for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+			{
+				if ( distance < fabs(this->heat_map[i][j] - second.heat_map[i][j]) )
+				{
+					distance = fabs(this->heat_map[i][j] - second.heat_map[i][j]);
+				}
+			}
+		}
+	}
+	return distance;
+}
+
+template <typename Scalling_of_kernels>
+double Persistence_heat_maps<Scalling_of_kernels>::project_to_R( int number_of_function )const
+{
+	double result = 0;
+	for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+	{
+		for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+		{
+			result += this->heat_map[i][j];
+		}
+	}
+	return result;
+}
+
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_average( const std::vector< Persistence_heat_maps* >& to_average )
+{				
+	this->compute_mean( to_average );
+}
+
+template <typename Scalling_of_kernels>
+double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product( const Persistence_heat_maps& second )const
+{		
+	//first we need to check if (*this) and second are defined on the same domain and have the same dimensions:
+	if ( !this->check_if_the_same(second) )
+	{
+		std::cerr << "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate";
+		throw "The persistence images are of noncompatible sizes. We cannot therefore compute distance between them. The program will now terminate";
+	}
+	
+	//if we are here, we know that the two persistence iomages are defined on the same domain, so we can start computing their scalar product:
+	double scalar_prod = 0;
+	for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
+	{
+		for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+		{
+			scalar_prod += this->heat_map[i][j]*second.heat_map[i][j];
+		}
+	}
+	return scalar_prod;
+}
+
+
+
+
+}//namespace Gudhi_stat
+}//namespace Gudhi
+
+
+#endif
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_landscape.h b/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_landscape.h
new file mode 100644
index 00000000..c71eb3d2
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_landscape.h
@@ -0,0 +1,1492 @@
+/*    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 (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/read_persitence_from_file.h>
+#include <gudhi/common_gudhi_stat.h>
+
+
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+
+
+
+
+
+/**
+ * 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.
+ * It implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product
+**/
+class Persistence_landscape
+{
+public:
+	/**
+	 * Default constructor.
+	**/
+    Persistence_landscape()
+    {
+		this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+	}
+
+    /**
+	 * Constructor that takes as an input a vector of birth-death pairs.
+	**/
+    Persistence_landscape( const std::vector< std::pair< double , double > >& p );
+
+    /**
+	 * 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_a_level_of_a_landscape( size_t level )const;
+
+
+    /**
+	 * This function compute integral of the landscape p-th power of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions)
+	**/
+	double compute_integral_of_landscape( double p )const;//this function compute integral of p-th power of landscape.
+
+
+    /**
+     * A function that computes the value of a landscape at a given point. The parameters of the function are: unsigned level and double x.
+     * The procedure will compute the value of the level-landscape at the point x.
+    **/
+    double compute_value_at_a_given_point( unsigned level , double x )const;
+
+    /**
+     * Writing landscape into a stream. A i-th level landscape starts with a string "lambda_i". Then the discontinuity points of the landscapes follows.
+     * Shall those points be joined with lines, we will obtain the i-th landscape function.
+    **/
+    friend std::ostream& operator<<(std::ostream& out, Persistence_landscape& land );
+
+
+
+
+
+	/**
+	 *\private A function that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape add_two_landscapes ( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 )
+    {
+        return operation_on_pair_of_landscapes< std::plus<double> >(land1,land2);
+    }
+
+    /**
+	 *\private A function that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape subtract_two_landscapes ( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 )
+    {
+        return operation_on_pair_of_landscapes< std::minus<double> >(land1,land2);
+    }
+
+	/**
+	 * An operator +, that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape operator+( const Persistence_landscape& first , const Persistence_landscape& second )
+    {
+        return add_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator -, that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape operator-( const Persistence_landscape& first , const Persistence_landscape& second )
+    {
+        return subtract_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator * that allows 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 persistence landscape.
+	**/
+    Persistence_landscape operator += ( const Persistence_landscape& rhs )
+    {
+        *this = *this + rhs;
+        return *this;
+    }
+
+	/**
+	 * Operator -=. The second parameter is a persistence landscape.
+	**/
+    Persistence_landscape operator -= ( const Persistence_landscape& rhs )
+    {
+        *this = *this - rhs;
+        return *this;
+    }
+
+
+	/**
+	 * Operator *=. The second parameter is a real number by which the y values of all landscape functions are multiplied. The x-values remain unchanged. 
+	**/
+    Persistence_landscape operator *= ( double x )
+    {
+        *this = *this*x;
+        return *this;
+    }
+
+	/**
+	 * Operator /=. The second parameter is a real number.
+	**/
+    Persistence_landscape operator /= ( double x )
+    {
+        if ( x == 0 )throw( "In operator /=, division by 0. Program terminated." );
+        *this = *this * (1/x);
+        return *this;
+    }
+
+	/**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator == ( const Persistence_landscape& rhs  )const;
+
+
+    /**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator != ( const Persistence_landscape& rhs  )const
+    {
+		return !((*this) == rhs);
+	}
+
+
+	/**
+	 * Computations of maximum (y) value of landscape.
+	**/
+    double compute_maximum()const
+    {
+        double maxValue = 0;
+        if ( this->land.size() )
+        {
+            maxValue = -std::numeric_limits<int>::max();
+            for ( size_t i = 0 ; i != this->land[0].size() ; ++i )
+            {
+                if ( this->land[0][i].second > maxValue )maxValue = this->land[0][i].second;
+            }
+        }
+        return maxValue;
+    }
+    
+    
+    /**
+	 *\private Computations of minimum (y) value of landscape.
+	**/
+    double compute_minimum()const
+    {
+        double minValue = 0;
+        if ( this->land.size() )
+        {
+            minValue = std::numeric_limits<int>::max();
+            for ( size_t i = 0 ; i != this->land[0].size() ; ++i )
+            {
+                if ( this->land[0][i].second < minValue )minValue = this->land[0][i].second;
+            }
+        }
+        return minValue;
+    }
+
+	/**
+	 *\private Computations of a L^i norm of landscape, where i is the input parameter.
+	**/
+    double compute_norm_of_landscape( double i )
+    {
+        Persistence_landscape l;
+        if ( i < std::numeric_limits< double >::max() )
+        {
+            return compute_distance_of_landscapes(*this,l,i);
+        }
+        else
+        {
+            return compute_max_norm_distance_of_landscapes(*this,l);
+        }
+    }
+
+	/**
+ 	 * An operator to compute the value of a landscape in the level 'level' at the argument 'x'.
+	**/
+    double operator()(unsigned level,double x)const{return this->compute_value_at_a_given_point(level,x);}
+
+	/**
+	 *\private Computations of L^{\infty} distance between two landscapes.
+	**/
+    friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second );
+    //friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , unsigned& nrOfLand , double&x , double& y1, double& y2 );
+
+
+	/**
+	 *\private Computations of L^{p} distance between two landscapes. p is the parameter of the procedure.
+	**/
+    friend double compute_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , double p );
+
+
+
+	/**
+	 * Function to compute absolute value of a PL function. The representation of persistence landscapes allow to store general PL-function. When computing distance 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;
+
+	/**
+	 *\private Function to compute inner (scalar) product of two landscapes.
+	**/
+    friend double compute_inner_product( const Persistence_landscape& l1 , const Persistence_landscape& l2 );
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+	//Implementations of functions for various concepts. 
+
+	/**
+	 * The number of projections to R is defined to the number of nonzero landscape functions. I-th projection is an integral of i-th landscape function over whole R.
+	 * This function is required by the Real_valued_topological_data concept.
+	**/
+    double project_to_R( int number_of_function )const
+    {
+		return this->compute_integral_of_a_level_of_a_landscape( (size_t)number_of_function );
+	}
+	
+	/**
+	 * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept.
+	**/ 
+	size_t number_of_projections_to_R()const
+	{
+		return this->number_of_functions_for_projections_to_reals;
+	}
+
+	/**
+	 * This function produce a vector of doubles based on a landscape. It is required in a concept Vectorized_topological_data
+	*/
+    std::vector<double> vectorize( int number_of_function )const
+    {
+		//TODO, think of something smarter over here
+		std::vector<double> v;
+		if ( (size_t)number_of_function > this->land.size() )
+		{
+			return v;
+		}
+		v.reserve( this->land[number_of_function].size() );
+		for ( size_t i = 0 ; i != this->land[number_of_function].size() ; ++i )
+		{
+			v.push_back( this->land[number_of_function][i].second );
+		}
+		return v;
+	}
+	/**
+	 * This function return the number of functions that allows vectorization of persistence laandscape. It is required in a concept Vectorized_topological_data.
+	 **/ 
+	size_t number_of_vectorize_functions()const
+	{
+		return this->number_of_functions_for_vectorization;	
+	}
+	
+	/**
+	 * A function to compute averaged persistence landscape, based on vector of persistence landscapes.
+	 * This function is required by Topological_data_with_averages concept.
+	**/ 
+    void compute_average( const std::vector< Persistence_landscape* >& to_average )
+    {
+		bool dbg = false;		
+		
+		if ( dbg ){std::cerr << "to_average.size() : " << to_average.size() << std::endl;}
+
+		std::vector< Persistence_landscape* > nextLevelMerge( to_average.size() );
+        for ( size_t i = 0 ; i != to_average.size() ; ++i )
+        {
+            nextLevelMerge[i] = to_average[i];            
+        }
+        bool is_this_first_level = true;//in the loop, we will create dynamically a 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 ){std::cerr << "nextLevelMerge.size() : " << nextLevelMerge.size() << std::endl;}
+            std::vector< Persistence_landscape* > nextNextLevelMerge;
+            nextNextLevelMerge.reserve( to_average.size() );
+            for ( size_t i = 0 ; i < nextLevelMerge.size() ; i=i+2 )
+            {
+                if ( dbg ){std::cerr << "i : " << i << std::endl;}
+                Persistence_landscape* l = new Persistence_landscape;
+                if ( i+1 != nextLevelMerge.size() )
+                {					
+                    (*l) = (*nextLevelMerge[i])+(*nextLevelMerge[i+1]);                    
+                }
+                else
+                {
+                    (*l) = *nextLevelMerge[i];
+                }
+                nextNextLevelMerge.push_back( l );
+            }
+            if ( dbg ){std::cerr << "After this iteration \n";getchar();}
+
+            if ( !is_this_first_level )
+            {
+				//deallocate the memory if the vector nextLevelMerge do not consist of the initial landscapes
+				for ( size_t i = 0 ; i != nextLevelMerge.size() ; ++i )
+				{
+					delete nextLevelMerge[i];
+				}
+			}
+			is_this_first_level = false;
+            nextLevelMerge.swap(nextNextLevelMerge);
+        }
+        (*this) = (*nextLevelMerge[0]);
+        (*this) *= 1/( (double)to_average.size() );
+	}
+
+
+	/**
+	* A function to compute distance between persistence landscape.
+	* The parameter of this functionis a Persistence_landscape.
+	* This function is required in Topological_data_with_distances concept.
+	* For max norm distance, set power to std::numeric_limits<double>::max()
+	**/
+    double distance( const Persistence_landscape& second , double power = 1 )const
+    {
+		if ( power < std::numeric_limits<double>::max() )
+		{
+			return compute_distance_of_landscapes( *this , second , power );
+		}
+		else
+		{
+			return compute_max_norm_distance_of_landscapes( *this , second );
+		}
+	}
+
+
+	/**
+	* A function to compute scalar product of persistence landscapes.
+	* The parameter of this functionis a Persistence_landscape.
+	* This function is required in Topological_data_with_scalar_product concept.
+	**/
+	double compute_scalar_product( const Persistence_landscape& second )const
+	{
+		return compute_inner_product( (*this) , second );
+	}	
+	//end of implementation of functions needed for concepts.
+
+
+	/**
+	 * This procedure returns x-range of a given level persistence landscape. If a default value is used, the x-range
+	 * of 0th level landscape is given (and this range contains the ranges of all other landscapes). 
+	**/ 
+	std::pair< double , double > get_x_range( size_t level = 0 )const
+	{
+		std::pair< double , double > result;
+		if ( level < this->land.size() )
+		{
+			result = std::make_pair( this->land[level][1].first , this->land[level][ this->land[level].size() - 2 ].first );
+		}
+		else
+		{
+			result = std::make_pair( 0,0 );
+		}
+		return result;
+	}
+	
+	/**
+	 * This procedure returns y-range of a given level persistence landscape. If a default value is used, the y-range
+	 * of 0th level landscape is given (and this range contains the ranges of all other landscapes). 
+	**/ 
+	std::pair< double , double > get_y_range( size_t level = 0 )const
+	{
+		std::pair< double , double > result;
+		if ( level < this->land.size() )
+		{
+			double maxx = this->compute_maximum();
+			double minn = this->compute_minimum();
+			result = std::make_pair( minn , maxx );
+		}
+		else
+		{
+			result = std::make_pair( 0,0 );
+		}
+		return result;
+	}
+
+	
+	
+	//a function used to create a gnuplot script for visualization of landscapes
+	void plot( const char* filename,  double xRangeBegin = std::numeric_limits<double>::max() , double xRangeEnd = std::numeric_limits<double>::max() , 
+									  double yRangeBegin = std::numeric_limits<double>::max() , double yRangeEnd = std::numeric_limits<double>::max(),
+									  int from = std::numeric_limits<int>::max(), int to = std::numeric_limits<int>::max()  );
+
+
+protected:
+    std::vector< std::vector< std::pair<double,double> > > land;
+    size_t number_of_functions_for_vectorization;
+	size_t number_of_functions_for_projections_to_reals;
+
+    void construct_persistence_landscape_from_barcode( const std::vector< std::pair< double , double > > & p );
+    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 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) ) )
+            {
+				//std::cerr<< this->land[level][i].first << " , " <<  rhs.land[level][i].first << " and " << this->land[level][i].second << " , " << rhs.land[level][i].second << std::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( 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 ){std::cerr << "Persistence_landscape::Persistence_landscape( const std::vector< std::pair< double , double > >& p )" << std::endl;}
+
+	//this is a general algorithm to construct persistence landscapes.
+	std::vector< std::pair<double,double> > bars;
+	bars.insert( bars.begin() , p.begin() , p.end() );
+	std::sort( bars.begin() , bars.end() , compare_points_sorting );
+
+	if (dbg)
+	{
+		std::cerr << "Bars : \n";
+		for ( size_t i = 0 ; i != bars.size() ; ++i )
+		{
+			std::cerr << bars[i].first << " " << bars[i].second << "\n";
+		}
+		getchar();
+	}
+
+	std::vector< std::pair<double,double> > characteristicPoints(p.size());
+	for ( size_t i = 0 ; i != bars.size() ; ++i )
+	{
+		characteristicPoints[i] = std::make_pair((bars[i].first+bars[i].second)/2.0 , (bars[i].second - bars[i].first)/2.0);
+	}
+	std::vector< std::vector< std::pair<double,double> > > Persistence_landscape;
+	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_a_level_of_a_landscape( size_t  level )const
+{
+    double result = 0;
+    if ( level >= this->land.size() )
+    {
+		//this landscape function is constantly equal 0, so is the 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 << "Here \n";
+        std::cerr << "x : " << x << "\n";
+        std::cerr << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n";
+        std::cerr << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n";
+    }
+
+    //in this case x is outside the support of the landscape, therefore the value of the landscape is 0.
+    if ( x <= this->land[level][coordBegin].first )return 0;
+    if ( x >= this->land[level][coordEnd].first )return 0;
+
+    if (compute_value_at_a_given_pointDbg)std::cerr << "Entering to the while loop \n";
+
+    while ( coordBegin+1 != coordEnd )
+    {
+        if (compute_value_at_a_given_pointDbg)
+        {
+            std::cerr << "coordBegin : " << coordBegin << "\n";
+            std::cerr << "coordEnd : " << coordEnd << "\n";
+            std::cerr << "this->land[level][coordBegin].first : " << this->land[level][coordBegin].first << "\n";
+            std::cerr << "this->land[level][coordEnd].first : " << this->land[level][coordEnd].first << "\n";
+        }
+
+
+        unsigned newCord = (unsigned)floor((coordEnd+coordBegin)/2.0);
+
+        if (compute_value_at_a_given_pointDbg)
+        {
+            std::cerr << "newCord : " << newCord << "\n";
+            std::cerr << "this->land[level][newCord].first : " << this->land[level][newCord].first << "\n";
+            std::cin.ignore();
+        }
+
+        if ( this->land[level][newCord].first <= x )
+        {
+            coordBegin = newCord;
+            if ( this->land[level][newCord].first == x )return this->land[level][newCord].second;
+        }
+        else
+        {
+            coordEnd = newCord;
+        }
+    }
+
+    if (compute_value_at_a_given_pointDbg)
+    {
+        std::cout << "x : " << x << " is between : " << this->land[level][coordBegin].first << " a  " << this->land[level][coordEnd].first << "\n";
+        std::cout << "the y coords are : " << this->land[level][coordBegin].second << " a  " << this->land[level][coordEnd].second << "\n";
+        std::cerr << "coordBegin : " << coordBegin << "\n";
+        std::cerr << "coordEnd : " << coordEnd << "\n";
+        std::cin.ignore();
+    }
+    return function_value( this->land[level][coordBegin] , this->land[level][coordEnd] , x );
+}
+
+std::ostream& operator<<(std::ostream& out, Persistence_landscape& land )
+{
+    for ( size_t level = 0 ; level != land.land.size() ; ++level )
+    {
+        out << "Lambda_" << level << ":" << std::endl;
+        for ( size_t i = 0 ; i != land.land[level].size() ; ++i )
+        {
+            if ( land.land[level][i].first == -std::numeric_limits<int>::max() )
+            {
+                out << "-inf";
+            }
+            else
+            {
+                if ( land.land[level][i].first == std::numeric_limits<int>::max() )
+                {
+                    out << "+inf";
+                }
+                else
+                {
+                    out << land.land[level][i].first;
+                }
+            }
+            out << " , " << land.land[level][i].second << std::endl;
+        }
+    }
+    return out;
+}
+
+
+
+
+void Persistence_landscape::multiply_lanscape_by_real_number_overwrite( double x )
+{
+    for ( size_t dim = 0 ; dim != this->land.size() ; ++dim )
+    {
+        for ( size_t i = 0 ; i != this->land[dim].size() ; ++i )
+        {
+             this->land[dim][i].second *= x;
+        }
+    }
+}
+
+bool AbsDbg = false;
+Persistence_landscape Persistence_landscape::abs()
+{
+    Persistence_landscape result;
+    for ( size_t level = 0 ; level != this->land.size() ; ++level )
+    {
+        if ( AbsDbg ){ std::cout << "level: " << level << std::endl; }
+        std::vector< std::pair<double,double> > lambda_n;
+        lambda_n.push_back( std::make_pair( -std::numeric_limits<int>::max() , 0 ) );
+        for ( size_t i = 1 ; i != this->land[level].size() ; ++i )
+        {
+            if ( AbsDbg ){std::cout << "this->land[" << level << "][" << i << "] : " << this->land[level][i].first << " " << this->land[level][i].second << std::endl;}
+            //if a line segment between this->land[level][i-1] and this->land[level][i] crosses the x-axis, then we have to add one landscape point t 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 ) )
+	{
+		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> > 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;
+    
+    if ( operation_on_pair_of_landscapesDBG )
+    {
+		for ( size_t i = 0 ; i != std::min( land1.land.size() , land2.land.size() ) ; ++i )
+		{
+			std::cerr << "land1.land[" << i << "].size() : " << land1.land[i].size() << std::endl;
+            std::cerr << "land2.land[" << i << "].size() : " << land2.land[i].size() << std::endl;             
+		}
+		getchar();
+	}
+
+    for ( size_t i = 0 ; i != std::min( land1.land.size() , land2.land.size() ) ; ++i )
+    {
+        std::vector< std::pair<double,double> > lambda_n;
+        size_t p = 0;
+        size_t q = 0;
+        while ( (p+1 < land1.land[i].size()) && (q+1 < land2.land[i].size()) )
+        {
+            if ( operation_on_pair_of_landscapesDBG )
+            {
+                std::cerr << "p : " << p << "\n";
+                std::cerr << "q : " << q << "\n";
+                std::cerr << "land1.land.size() : " << land1.land.size() << std::endl;
+                std::cerr << "land2.land.size() : " << land2.land.size() << std::endl;
+                std::cerr << "land1.land[" << i << "].size() : " << land1.land[i].size() << std::endl;
+                std::cerr << "land2.land[" << i << "].size() : " << land2.land[i].size() << std::endl;
+                std::cout << "land1.land[i][p].first : " << land1.land[i][p].first << "\n";
+                std::cout << "land2.land[i][q].first : " << land2.land[i][q].first << "\n";
+                //getchar();
+            }
+
+            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";}
+        }
+        while ( (p+1 < land1.land[i].size())&&(q+1 >= land2.land[i].size()) )
+        {
+            if (operation_on_pair_of_landscapesDBG)
+            {
+                std::cout << "New point : " << land1.land[i][p].first << "  oper(land1.land[i][p].second,0) : " <<  oper(land1.land[i][p].second,0) << std::endl;
+            }
+            lambda_n.push_back( std::make_pair(land1.land[i][p].first , oper(land1.land[i][p].second,0) ) );
+            ++p;
+        }
+        while ( (p+1 >= land1.land[i].size())&&(q+1 < land2.land[i].size()) )
+        {
+            if (operation_on_pair_of_landscapesDBG)
+            {
+                std::cout << "New point : " << land2.land[i][q].first << " oper(0,land2.land[i][q].second) : " <<  oper(0,land2.land[i][q].second) << std::endl;
+            }
+            lambda_n.push_back( std::make_pair(land2.land[i][q].first , oper(0,land2.land[i][q].second) ) );
+            ++q;
+        }
+        lambda_n.push_back( std::make_pair( std::numeric_limits<int>::max() , 0 ) );
+        //CHANGE
+        //result.land[i] = lambda_n;
+        result.land[i].swap(lambda_n);
+    }
+    if ( land1.land.size() > std::min( land1.land.size() , land2.land.size() ) )
+    {
+        if (operation_on_pair_of_landscapesDBG){std::cout << "land1.land.size() > std::min( land1.land.size() , land2.land.size() )" << std::endl;}
+        for ( size_t i = std::min( land1.land.size() , land2.land.size() ) ; i != std::max( land1.land.size() , land2.land.size() ) ; ++i )
+        {
+            std::vector< std::pair<double,double> > lambda_n( land1.land[i] );
+            for ( size_t nr = 0 ; nr != land1.land[i].size() ; ++nr )
+            {
+                lambda_n[nr] = std::make_pair( land1.land[i][nr].first , oper( land1.land[i][nr].second , 0 ) );
+            }
+            //CHANGE
+            //result.land[i] = lambda_n;
+            result.land[i].swap(lambda_n);
+        }
+    }
+    if ( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) )
+    {
+        if (operation_on_pair_of_landscapesDBG){std::cout << "( land2.land.size() > std::min( land1.land.size() , land2.land.size() ) ) " << std::endl;}
+        for ( size_t i = std::min( land1.land.size() , land2.land.size() ) ; i != std::max( land1.land.size() , land2.land.size() ) ; ++i )
+        {
+            std::vector< std::pair<double,double> > lambda_n( land2.land[i] );
+            for ( size_t nr = 0 ; nr != land2.land[i].size() ; ++nr )
+            {
+                lambda_n[nr] = std::make_pair( land2.land[i][nr].first , oper( 0 , land2.land[i][nr].second ) );
+            }
+            //CHANGE
+            //result.land[i] = lambda_n;
+            result.land[i].swap(lambda_n);
+        }
+    }
+    if ( operation_on_pair_of_landscapesDBG ){std::cout << "operation_on_pair_of_landscapes END\n";std::cin.ignore();}
+    return result;
+}//operation_on_pair_of_landscapes
+
+
+
+double compute_maximal_distance_non_symmetric( const Persistence_landscape& pl1, const Persistence_landscape& pl2 )
+{
+    bool dbg = false;
+    if (dbg)std::cerr << " compute_maximal_distance_non_symmetric \n";
+    //this distance is not symmetric. It compute ONLY distance between inflection points of pl1 and pl2.
+    double maxDist = 0;
+    size_t minimalNumberOfLevels = std::min( pl1.land.size() , pl2.land.size() );
+    for ( size_t level = 0 ; level != minimalNumberOfLevels ; ++ level )
+    {
+        if (dbg)
+        {
+            std::cerr << "Level : " << level << std::endl;
+            std::cerr << "PL1 : \n";
+            for ( size_t i = 0 ; i  != pl1.land[level].size() ; ++i )
+            {
+                std::cerr << "(" <<pl1.land[level][i].first << "," << pl1.land[level][i].second << ") \n";
+            }
+            std::cerr << "PL2 : \n";
+            for ( size_t i = 0 ; i  != pl2.land[level].size() ; ++i )
+            {
+                std::cerr << "(" <<pl2.land[level][i].first << "," << pl2.land[level][i].second << ") \n";
+            }
+            std::cin.ignore();
+        }
+
+        int p2Count = 0;
+        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_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second , double p )
+{
+	bool dbg = false;
+    //This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a time:
+
+    //first-second :
+    Persistence_landscape lan = first-second;
+
+    //| first-second |:
+    lan = lan.abs();
+    
+    if ( dbg ){std::cerr << "Abs of difference ; " << lan << std::endl;getchar();}
+    
+	if ( p < std::numeric_limits<double>::max() )
+	{
+		//\int_{- \infty}^{+\infty}| first-second |^p
+		double result;
+		if ( p != 1 )
+		{
+			if ( dbg )std::cerr << "Power != 1, compute integral to the power p\n";
+			result = lan.compute_integral_of_landscape( (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 == infty
+		if ( dbg )std::cerr << "Power = infty, compute maximum \n";
+		return lan.compute_maximum();
+	}
+}
+
+double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second )
+{
+    return std::max( compute_maximal_distance_non_symmetric(first,second) , compute_maximal_distance_non_symmetric(second,first) );
+}
+
+
+bool comparePairsForMerging( std::pair< double , unsigned > first , std::pair< double , unsigned > second )
+{
+    return (first.first < second.first);
+}
+
+
+
+
+double compute_inner_product( const Persistence_landscape& l1 , const Persistence_landscape& l2 )
+{
+    bool dbg = false;
+    double result = 0;
+
+    for ( size_t level = 0 ; level != std::min( l1.size() , l2.size() ) ; ++level )
+    {
+        if ( dbg ){std::cerr << "Computing inner product for a level : " << level << std::endl;getchar();}
+        if ( l1.land[level].size() * l2.land[level].size() == 0 )continue;
+
+        //endpoints of the interval on which we will compute the inner product of two locally linear functions:
+        double x1 = -std::numeric_limits<int>::max();
+        double x2;
+        if ( l1.land[level][1].first < l2.land[level][1].first )
+        {
+            x2 = l1.land[level][1].first;
+        }
+        else
+        {
+            x2 = l2.land[level][1].first;
+        }
+
+        //iterators for the landscapes l1 and l2
+        size_t l1It = 0;
+        size_t l2It = 0;
+
+        while ( (l1It < l1.land[level].size()-1) && (l2It < l2.land[level].size()-1) )
+        {
+            //compute the value of a inner product on a interval [x1,x2]
+
+            double a,b,c,d;
+
+			if ( l1.land[level][l1It+1].first != l1.land[level][l1It].first )
+			{
+				a = (l1.land[level][l1It+1].second - l1.land[level][l1It].second)/(l1.land[level][l1It+1].first - l1.land[level][l1It].first);
+			}
+			else
+			{
+				a = 0;
+			}
+            b = l1.land[level][l1It].second - a*l1.land[level][l1It].first;
+            if ( l2.land[level][l2It+1].first != l2.land[level][l2It].first )
+            {
+				c = (l2.land[level][l2It+1].second - l2.land[level][l2It].second)/(l2.land[level][l2It+1].first - l2.land[level][l2It].first);
+			}
+			else
+			{
+				c = 0;
+			}
+            d = l2.land[level][l2It].second - c*l2.land[level][l2It].first;
+
+            double contributionFromThisPart
+            =
+            (a*c*x2*x2*x2/3 + (a*d+b*c)*x2*x2/2 + b*d*x2) - (a*c*x1*x1*x1/3 + (a*d+b*c)*x1*x1/2 + b*d*x1);
+
+            result += contributionFromThisPart;
+
+            if ( dbg )
+            {
+                std::cerr << "[l1.land[level][l1It].first,l1.land[level][l1It+1].first] : " << l1.land[level][l1It].first << " , " << l1.land[level][l1It+1].first << std::endl;
+                std::cerr << "[l2.land[level][l2It].first,l2.land[level][l2It+1].first] : " << l2.land[level][l2It].first << " , " << l2.land[level][l2It+1].first << std::endl;
+                std::cerr << "a : " << a << ", b : " << b << " , c: " << c << ", d : " << d << std::endl;
+                std::cerr << "x1 : " << x1 << " , x2 : " << x2 << std::endl;
+                std::cerr << "contributionFromThisPart : " << contributionFromThisPart << std::endl;
+                std::cerr << "result : " << result << std::endl;
+                getchar();
+            }
+
+            //we have two intervals in which functions are constant:
+            //[l1.land[level][l1It].first , l1.land[level][l1It+1].first]
+            //and
+            //[l2.land[level][l2It].first , l2.land[level][l2It+1].first]
+            //We also have an interval [x1,x2]. Since the intervals in the landscapes cover the whole R, then it is clear that x2
+            //is either l1.land[level][l1It+1].first of l2.land[level][l2It+1].first or both. Lets test it.
+            if ( x2 == l1.land[level][l1It+1].first )
+            {
+                if ( x2 == l2.land[level][l2It+1].first )
+                {
+                    //in this case, we increment both:
+                    ++l2It;
+                    if ( dbg ){std::cerr << "Incrementing both \n";}
+                }
+                else
+                {
+                    if ( dbg ){std::cerr << "Incrementing first \n";}
+                }
+                ++l1It;
+            }
+            else
+            {
+                //in this case we increment l2It
+                ++l2It;
+                if ( dbg ){std::cerr << "Incrementing second \n";}
+            }
+            //Now, we shift x1 and x2:
+            x1 = x2;
+            if ( l1.land[level][l1It+1].first < l2.land[level][l2It+1].first )
+            {
+                x2 = l1.land[level][l1It+1].first;
+            }
+            else
+            {
+                x2 = l2.land[level][l2It+1].first;
+            }
+
+        }
+
+    }
+    return result;
+}
+
+
+void Persistence_landscape::plot( const char* filename,  double xRangeBegin , double xRangeEnd , double yRangeBegin , double yRangeEnd , int from , int to )
+{
+    //this program create a gnuplot script file that allows to plot persistence diagram.
+    std::ofstream out;
+
+    std::ostringstream nameSS;
+    nameSS << filename << "_GnuplotScript";
+    std::string nameStr = nameSS.str();
+    out.open( nameStr );
+
+    if ( (xRangeBegin != std::numeric_limits<double>::max()) || (xRangeEnd != std::numeric_limits<double>::max()) || (yRangeBegin != std::numeric_limits<double>::max()) || (yRangeEnd != std::numeric_limits<double>::max())  )
+    {
+        out << "set xrange [" << xRangeBegin << " : " << xRangeEnd << "]" << std::endl;
+        out << "set yrange [" << yRangeBegin << " : " << yRangeEnd << "]" << std::endl;
+    }
+
+    if ( from == std::numeric_limits<int>::max() ){from = 0;}
+    if ( to == std::numeric_limits<int>::max() ){to = this->land.size();}
+
+    out << "plot ";
+    for ( size_t lambda= std::min((size_t)from,this->land.size()) ; lambda != std::min((size_t)to,this->land.size()) ; ++lambda )
+    {
+        //out << "     '-' using 1:2 title 'l" << lambda << "' with lp";
+        out << "     '-' using 1:2 notitle with lp";
+        if ( lambda+1 != std::min((size_t)to,this->land.size()) )
+        {
+            out << ", \\";
+        }
+        out << std::endl;
+    }
+
+    for ( size_t lambda= std::min((size_t)from,this->land.size()) ; lambda != std::min((size_t)to,this->land.size()) ; ++lambda )
+    {
+        for ( size_t i = 1 ; i != this->land[lambda].size()-1 ; ++i )
+        {
+            out << this->land[lambda][i].first << " " << this->land[lambda][i].second << std::endl;
+        }
+        out << "EOF" << std::endl;
+    }
+    std::cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << std::endl;
+}
+
+
+
+
+}//namespace gudhi stat
+}//namespace gudhi
+
+
+#endif
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_landscape_on_grid.h b/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_landscape_on_grid.h
new file mode 100644
index 00000000..d29a51ab
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/Persistence_landscape_on_grid.h
@@ -0,0 +1,1363 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *h
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2015  INRIA (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_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/read_persitence_from_file.h>
+#include <gudhi/common_gudhi_stat.h>
+
+
+
+namespace Gudhi
+{
+namespace Gudhi_stat
+{
+
+
+//this class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product
+class Persistence_landscape_on_grid  
+{
+public:
+	/**
+	 * Default constructor.
+	**/
+    Persistence_landscape_on_grid()
+    {
+		this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+		this->grid_min = this->grid_max = 0;
+	}
+
+    /**
+	 * Constructor that takes as an input a vector of birth-death pairs.
+	**/
+    Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ );
+    
+    /**
+	 * Constructor that 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 )
+		{			
+			std::cerr << "this->grid_max : " << this->grid_max << std::endl;
+			std::cerr << "this->grid_min : " << this->grid_min << std::endl;
+			std::cerr << "this->values_of_landscapes.size() : " << this->values_of_landscapes.size() << std::endl;
+			getchar();
+		}
+		
+		
+		double previous_x = this->grid_min-dx;
+		double previous_y = 0;
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			double current_x = previous_x + dx;
+			double current_y = 0;
+			if ( this->values_of_landscapes[i].size() > level )current_y = this->values_of_landscapes[i][level];
+		
+			if ( dbg )
+			{
+				std::cerr << "this->values_of_landscapes[i].size() : " << this->values_of_landscapes[i].size() << " , level : " << level << std::endl;	
+				if ( this->values_of_landscapes[i].size() > level )std::cerr << "this->values_of_landscapes[i][level] : " << this->values_of_landscapes[i][level] << std::endl;
+				std::cerr << "previous_y : " << previous_y << std::endl;
+				std::cerr << "current_y : " << current_y << std::endl;
+				std::cerr << "dx : " << dx << std::endl;
+				std::cerr << "0.5*dx*( previous_y + current_y ); " << 0.5*dx*( previous_y + current_y ) << std::endl;
+			}
+			
+			result += 0.5*dx*( previous_y + current_y );			
+			previous_x = current_x;
+			previous_y = current_y;
+		}
+		return result;
+	}
+
+	/**
+	 * This function compute integral of the landscape p-th power of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions)
+	**/
+	double compute_integral_of_landscape( double p )const
+	{
+		size_t maximal_level = this->number_of_nonzero_levels();
+		double result = 0;
+		for ( size_t i = 0 ; i != maximal_level ; ++i )
+		{
+			result += this->compute_integral_of_landscape(p,i);
+		}
+		return result;
+	}
+
+    /**
+	 * This function compute integral of the landscape p-th power of a level of a landscape (defined formally as sum of integrals on R of p-th powers of all landscape functions)
+	**/	
+	double compute_integral_of_landscape( double p , size_t level )const
+	{
+		bool dbg = false;
+		
+		double result = 0;				
+		double dx = (this->grid_max - this->grid_min)/(double)(this->values_of_landscapes.size()-1);		
+		double previous_x = this->grid_min;
+		double previous_y = 0;
+		if ( this->values_of_landscapes[0].size() > level )previous_y = this->values_of_landscapes[0][level];
+		
+		if ( dbg )
+		{
+			std::cerr << "dx : " << dx << std::endl;
+			std::cerr << "previous_x : " << previous_x << std::endl;
+			std::cerr << "previous_y : " << previous_y << std::endl;
+			std::cerr << "power : " << p << std::endl;
+			getchar();
+		}
+		
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			double current_x = previous_x + dx;
+			double current_y = 0;
+			if ( this->values_of_landscapes[i].size() > level )current_y = this->values_of_landscapes[i][level];
+			
+			if ( dbg )std::cerr << "current_y : " << current_y << std::endl;
+			
+			if ( current_y == previous_y )continue; 
+			
+			std::pair<double,double> coef = compute_parameters_of_a_line( std::make_pair( previous_x , previous_y ) , std::make_pair( current_x , current_y ) );
+			double a = coef.first;
+			double b = coef.second;
+			
+			if ( dbg )
+			{
+				std::cerr << "A line passing through points : (" << previous_x << "," << previous_y << ") and (" << current_x << "," << current_y << ") is : " << a << "x+" << b << std::endl;
+			}
+			
+			//In this interval, the landscape has a form f(x) = ax+b. We want to compute integral of (ax+b)^p = 1/a * (ax+b)^{p+1}/(p+1)
+			double value_to_add = 0;
+			if ( a != 0 )
+			{			
+				value_to_add = 1/(a*(p+1)) * ( pow((a*current_x+b),p+1) - pow((a*previous_x+b),p+1));				
+			}
+			else
+			{
+				value_to_add = ( current_x - previous_x )*( pow(b,p) );
+			}				
+			result += value_to_add;
+			if ( dbg )
+			{
+				std::cerr << "Increasing result by : " << value_to_add << std::endl;
+				std::cerr << "restult : " << result << std::endl;
+				getchar();
+			}		
+			previous_x = current_x;
+			previous_y = current_y;
+		}		
+		if ( dbg )std::cerr << "The total result is : " << result << std::endl;
+		return result;
+	}
+	
+	/**
+     * Writing landscape into a stream. A i-th level landscape starts with a string "lambda_i". Then the discontinuity points of the landscapes follows.
+     * Shall those points be joined with lines, we will obtain the i-th landscape function.
+    **/
+    friend std::ostream& operator<<(std::ostream& out, const Persistence_landscape_on_grid& land )
+	{
+		double dx = (land.grid_max - land.grid_min)/(double)(land.values_of_landscapes.size()-1);
+		double x = land.grid_min;
+		for ( size_t i = 0 ; i != land.values_of_landscapes.size() ; ++i )
+		{
+			out << x << " : ";
+			for ( size_t j = 0 ; j != land.values_of_landscapes[i].size() ; ++j )
+			{
+				out << land.values_of_landscapes[i][j] << " ";
+			}
+			out << std::endl;
+			x += dx;
+		}
+		return out;
+	}
+
+
+    /**
+     * 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 << std::endl;
+			std::cerr << "x : " << x << std::endl;
+			std::cerr << "psoition : " << position << std::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;
+	}
+    
+
+	/**
+	 *\private A function that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid add_two_landscapes ( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 )
+    {
+        return operation_on_pair_of_landscapes_on_grid< std::plus<double> >(land1,land2);
+    }
+
+    /**
+	 *\private A function that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid subtract_two_landscapes ( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 )
+    {
+        return operation_on_pair_of_landscapes_on_grid< std::minus<double> >(land1,land2);
+    }
+
+	/**
+	 * An operator +, that compute sum of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid operator+( const Persistence_landscape_on_grid& first , const Persistence_landscape_on_grid& second )
+    {
+        return add_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator -, that compute difference of two landscapes.
+	**/
+    friend Persistence_landscape_on_grid operator-( const Persistence_landscape_on_grid& first , const Persistence_landscape_on_grid& second )
+    {
+        return subtract_two_landscapes( first,second );
+    }
+
+	/**
+	 * An operator * that allows 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) std::cerr << "values_of_landscapes of incompatable sizes\n";
+			return false;
+		}
+		if ( !almost_equal( this->grid_min , rhs.grid_min ) )
+		{
+			if (dbg) std::cerr << "grid_min not equal\n";
+			return false; 
+		}
+		if ( !almost_equal(this->grid_max,rhs.grid_max ) )
+		{	
+			if (dbg) std::cerr << "grid_max not equal\n";
+			return false;		
+		}
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{					
+			for ( size_t aa = 0 ; aa != this->values_of_landscapes[i].size() ; ++aa )
+			{
+				if ( !almost_equal( this->values_of_landscapes[i][aa]  , rhs.values_of_landscapes[i][aa]  ) )
+				{					
+					if (dbg) 
+					{
+						std::cerr << "Problem in the position : " << i << " of values_of_landscapes. \n";
+						std::cerr << this->values_of_landscapes[i][aa]   << " " <<  rhs.values_of_landscapes[i][aa] << std::endl;
+					}
+					return false;
+				}
+			}				
+		}	
+		return true;
+	}
+
+
+    /**
+	 * An operator to compare two persistence landscapes.
+	**/
+    bool operator != ( const Persistence_landscape_on_grid& rhs  )const
+    {
+		return !((*this) == rhs);
+	}
+
+
+	/**
+	 * Computations of maximum (y) value of landscape.
+	**/
+    double compute_maximum()const
+    {
+       //since the function can only be entirely positive or negative, the maximal value will be an extremal value in the arrays:
+		double max_value = -std::numeric_limits<double>::max();
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			if ( this->values_of_landscapes[i].size() )
+			{
+				if ( this->values_of_landscapes[i][0] > max_value )max_value = this->values_of_landscapes[i][0];
+				if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] > max_value )max_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ];
+			}
+		}
+		return max_value;        
+    }
+    
+    /**
+	 * Computations of minimum and maximum value of landscape.
+	**/
+    std::pair<double,double> compute_minimum_maximum()const
+    {        
+       //since the function can only be entirely positive or negative, the maximal value will be an extremal value in the arrays:
+		double max_value = -std::numeric_limits<double>::max();
+		double min_value = 0;
+		for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+		{
+			if ( this->values_of_landscapes[i].size() )
+			{
+				if ( this->values_of_landscapes[i][0] > max_value )max_value = this->values_of_landscapes[i][0];
+				if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] > max_value )max_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ];
+				
+				if ( this->values_of_landscapes[i][0] < min_value )min_value = this->values_of_landscapes[i][0];
+				if ( this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ] < min_value )min_value = this->values_of_landscapes[i][ this->values_of_landscapes[i].size()-1 ];
+			}
+		}
+		return std::make_pair(min_value , max_value);
+    }
+    
+    /**
+	 * This procedure returns x-range of a given level persistence landscape. If a default value is used, the x-range
+	 * of 0th level landscape is given (and this range contains the ranges of all other landscapes). 
+	**/ 
+	std::pair< double , double > get_x_range( size_t level = 0 )const
+	{
+		return std::make_pair( this->grid_min , this->grid_max );
+		//std::pair< double , double > result;
+		//if ( level < this->land.size() )
+		//{					
+		//	double dx = (this->grid_max - this->grid_min)/(double)this->values_of_landscapes.size();
+		//	size_t first_nonzero = 0;
+		//	while ( (first_nonzero != this->values_of_landscapes.size()) && (this->values_of_landscapes[level][first_nonzero] == 0) )++first_nonzero;
+		//							
+		//	if ( first_nonzero == 0 )
+		//	{
+		//		return std::make_pair( 0,0 );//this landscape is empty.
+		//	}						
+		//							
+		//	size_t last_nonzero = 0;
+		//	while ( (last_nonzero != 0) && (this->values_of_landscapes[level][last_nonzero] == 0) )--last_nonzero;
+		//	
+		//	result = std::make_pair( this->grid_min +first_nonzero*dx , this->grid_max - last_nonzero*dx   );
+		//}
+		//else
+		//{
+		//	result = std::make_pair( 0,0 );
+		//}
+		//return result;
+	}
+	
+	/**
+	 * This procedure returns y-range of a persistence landscape. If a default value is used, the y-range
+	 * of 0th level landscape is given (and this range contains the ranges of all other landscapes). 
+	**/ 
+	std::pair< double , double > get_y_range( size_t level = 0 )const
+	{
+		return this->compute_minimum_maximum();
+		//std::pair< double , double > result;
+		//if ( level < this->land.size() )
+		//{
+		//	result = this->compute_minimum_maximum()			
+		//}
+		//else
+		//{
+		//	result = std::make_pair( 0,0 );
+		//}
+		//return result;
+	}
+    
+    /**
+     * 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 )const
+    {
+		std::vector< std::pair< double , double > > p;
+        Persistence_landscape_on_grid l(p,this->grid_min,this->grid_max,this->values_of_landscapes.size()-1);       
+        
+        if ( i < std::numeric_limits<double>::max() )
+        {
+            return compute_distance_of_landscapes_on_grid(*this,l,i);
+        }
+        else
+        {
+            return compute_max_norm_distance_of_landscapes(*this,l);
+        }
+    }
+
+	/**
+ 	 * An operator to compute the value of a landscape in the level 'level' at the argument 'x'.
+	**/
+    double operator()(unsigned level,double x)const{return this->compute_value_at_a_given_point(level,x);}
+
+	/**
+	 * Computations of L^{\infty} distance between two landscapes.
+	**/
+    friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second );
+    //friend double compute_max_norm_distance_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()
+    {		
+		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;
+	}
+	
+	
+	/**
+	 * Computations of L^{p} distance between two landscapes on a grid. p is the parameter of the procedure.
+	 * FIXME: Note that, due to the grid representation, the method below may give non--accurate results in case when the landscape P and Q the difference of which we want to compute
+	 * are interxsecting. This is a consequence of a general way they are computed. In the future, an integral of absolute value of a difference of P and Q will be given as a separated
+	 * function to fix that inaccuracy. 
+	**/
+	friend double compute_distance_of_landscapes_on_grid( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second , double p )
+	{		
+		bool dbg = false;
+		//This is what we want to compute: (\int_{- \infty}^{+\infty}| first-second |^p)^(1/p). We will do it one step at a time:
+
+		if ( dbg )
+		{
+			std::cerr << "first : " << first << std::endl;
+			std::cerr << "second : " << second << std::endl;
+			getchar();
+		}
+
+		//first-second :
+		Persistence_landscape_on_grid lan = first-second;
+		
+		if ( dbg )
+		{
+			std::cerr << "Difference : " << lan << std::endl;
+		}
+
+		//| first-second |:
+		lan.abs();
+		
+		if ( dbg )
+		{
+			std::cerr << "Abs : " << lan << std::endl;
+		}
+		
+		if ( p < std::numeric_limits< double >::max() )
+		{
+			//\int_{- \infty}^{+\infty}| first-second |^p
+			double result;
+			if ( p != 1 )
+			{
+				if (dbg){std::cerr << "p : " << p << std::endl; getchar();}
+				result = lan.compute_integral_of_landscape( (double)p );				
+				if (dbg){std::cerr << "integral : " << result << std::endl;getchar();}
+			}
+			else
+			{
+				result = lan.compute_integral_of_landscape();			
+				if (dbg){std::cerr << "integral, wihtout power : " << result << std::endl;getchar();}
+			}
+			//(\int_{- \infty}^{+\infty}| first-second |^p)^(1/p)
+			return pow( result , 1/(double)p );
+		}
+		else
+		{
+			//p == infty
+			return lan.compute_maximum();
+		}
+	}
+	
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+	//Functions that are needed for that class to implement the concept.
+
+	/**
+	 * The number of projections to R is defined to the number of nonzero landscape functions. I-th projection is an integral of i-th landscape function over whole R.
+	 * This function is required by the Real_valued_topological_data concept. 
+	**/
+    double project_to_R( int number_of_function )const
+    {
+		return this->compute_integral_of_landscape( (size_t)number_of_function );
+	}
+	
+	/**
+	 * The function gives the number of possible projections to R. This function is required by the Real_valued_topological_data concept.
+	**/ 
+	size_t number_of_projections_to_R()const
+	{
+		return number_of_functions_for_projections_to_reals;
+	}
+
+
+
+
+	/**
+	 * This function produce a vector of doubles based on a landscape. It is required in a concept Vectorized_topological_data
+	*/ 
+    std::vector<double> vectorize( int number_of_function )const
+    {
+		//TODO, 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;
+	}
+	
+	/**
+	 * This function return the number of functions that allows vectorization of persistence laandscape. It is required in a concept Vectorized_topological_data.
+	 **/ 
+	size_t number_of_vectorize_functions()const
+	{
+		return number_of_functions_for_vectorization;	
+	}
+
+
+
+
+	
+	/**
+	 * A function to compute averaged persistence landscape on a grid, based on vector of persistence landscapes on grid.
+	 * This function is required by Topological_data_with_averages concept.
+	**/  
+    void compute_average( const std::vector< Persistence_landscape_on_grid* >& to_average )
+    {
+		
+		bool dbg = false;
+		//After execution of this procedure, the average is supposed to be in the current object. To make sure that this is the case, we need to do some cleaning first.
+		this->values_of_landscapes .clear();
+		this->grid_min = this->grid_max = 0;
+		
+		//if there is nothing to 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(*(to_average[0]),*(to_average[i])) )throw "Two grids are not compatible"; 				
+		}
+		
+		this->values_of_landscapes = std::vector< std::vector<double> >( (to_average[0])->values_of_landscapes.size() );
+		this->grid_min = (to_average[0])->grid_min;
+		this->grid_max = (to_average[0])->grid_max;
+		
+		if ( dbg )
+		{
+			std::cerr << "Computations of average. The data from the current landscape have been cleared. We are ready to do the computations. \n";
+		}
+						
+		//for every point in the grid:
+		for ( size_t grid_point = 0 ; grid_point != (to_average[0])->values_of_landscapes.size() ; ++grid_point )
+		{
+			
+			//set up a vector of the correct size:
+			size_t maximal_size_of_vector = 0;
+			for ( size_t land_no = 0 ; land_no != to_average.size() ; ++land_no )
+			{
+				if ( (to_average[land_no])->values_of_landscapes[grid_point].size() > maximal_size_of_vector )
+				maximal_size_of_vector = (to_average[land_no])->values_of_landscapes[grid_point].size();
+			}
+			this->values_of_landscapes[grid_point] = std::vector<double>( maximal_size_of_vector );
+			
+			if ( dbg )
+			{
+				std::cerr << "We are considering the point : " << grid_point << " of the grid. In this point, there are at most : " << maximal_size_of_vector << " nonzero landscape functions \n";
+			}
+		
+			//and compute an arythmetic average:
+			for ( size_t land_no = 0 ; land_no != to_average.size() ; ++land_no )
+			{
+				//summing:				
+				for ( size_t i = 0 ; i != (to_average[land_no])->values_of_landscapes[grid_point].size() ; ++i )
+				{
+					//compute the average in a smarter way.
+					this->values_of_landscapes[grid_point][i] += (to_average[land_no])->values_of_landscapes[grid_point][i];
+				}								
+			}
+			//normalizing:
+			for ( size_t i = 0 ; i != this->values_of_landscapes[grid_point].size() ; ++i )
+			{
+				this->values_of_landscapes[grid_point][i] /= (double)to_average.size();
+			}			 
+		}
+	}//compute_average
+	
+	
+	/**
+	* A function to compute distance between persistence landscape on a grid.
+	* The parameter of this functionis a Persistence_landscape_on_grid.
+	* This function is required in Topological_data_with_distances concept.
+	* For max norm distance, set power to std::numeric_limits<double>::max()
+	**/
+    double distance( const Persistence_landscape_on_grid& second , double power = 1 )const
+    {
+		if ( power < std::numeric_limits<double>::max() )
+		{
+			return compute_distance_of_landscapes_on_grid( *this , second , power );
+		}
+		else
+		{
+			return compute_max_norm_distance_of_landscapes( *this , second );
+		}
+	}
+
+	/**
+	* A function to compute scalar product of persistence landscape on a grid.
+	* The parameter of this functionis a Persistence_landscape_on_grid.
+	* This function is required in Topological_data_with_scalar_product concept.
+	**/
+	double compute_scalar_product( const Persistence_landscape_on_grid& second )
+	{
+		return compute_inner_product( (*this) , second );
+	}
+
+	//end of implementation of functions needed for concepts. 
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+	
+
+
+	/**
+	* A function that returns values of landsapes. It can be used for vizualization
+	**/
+	std::vector< std::vector< double > > output_for_visualization()const
+	{
+		return this->values_of_landscapes;
+	}
+	
+	/**
+	* function used to create a gnuplot script for visualization of landscapes. Over here we need to specify which landscapes do we want to plot.
+	* In addition, the user may specify the range (min and max) where landscape is plot. The fefault values for min and max are std::numeric_limits<double>::max(). If the procedure detect those
+	* values, it will determine the range so that the whole landscape is supported there. If at least one min or max value is different from std::numeric_limits<double>::max(), then the values
+	* provided by the user will be used. 
+	**/ 
+	void plot( const char* filename , size_t from_ , size_t to_ )const
+	{
+		this->plot( filename , std::numeric_limits<double>::max() , std::numeric_limits<double>::max(), std::numeric_limits<double>::max() , std::numeric_limits<double>::max() , from_ , to_ );
+	}
+	
+	/**
+	* function used to create a gnuplot script for visualization of landscapes. Over here we can restrict also x and y range of the landscape. 
+	**/ 
+	void plot( const char* filename, double min_x = std::numeric_limits<double>::max() , double max_x = std::numeric_limits<double>::max() , double min_y = std::numeric_limits<double>::max() , double max_y = std::numeric_limits<double>::max() , size_t from_ = std::numeric_limits<size_t>::max(), size_t to_= std::numeric_limits<size_t>::max()  )const;
+
+
+protected:
+	double grid_min;
+	double grid_max;
+	std::vector< std::vector< double > > values_of_landscapes;
+	size_t number_of_functions_for_vectorization;
+	size_t number_of_functions_for_projections_to_reals;
+	
+	void set_up_numbers_of_functions_for_vectorization_and_projections_to_reals()
+	{
+		//warning, this function can be only called after filling in the values_of_landscapes vector.
+		this->number_of_functions_for_vectorization = this->values_of_landscapes.size();
+		this->number_of_functions_for_projections_to_reals = this->values_of_landscapes.size();
+	}
+	
+	void set_up_values_of_landscapes( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ );
+	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_ << std::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 << std::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 )
+		{
+			std::cerr << "Considering an interval : " << p[int_no].first << "," << p[int_no].second << std::endl;
+		
+			std::cerr << "grid_interval_begin : " << grid_interval_begin << std::endl;
+			std::cerr << "grid_interval_end : " << grid_interval_end << std::endl;
+			std::cerr << "grid_interval_midpoint : " << grid_interval_midpoint << std::endl;	
+		}
+		
+		double landscape_value = dx;
+		for ( size_t i = grid_interval_begin+1 ; i < grid_interval_midpoint ; ++i )
+		{
+			if ( dbg )
+			{
+				std::cerr << "Adding landscape value (going up) for a point : " << i << " equal : " << landscape_value << std::endl;
+			}
+			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 values:
+	for ( size_t pt = 0 ; pt != this->values_of_landscapes.size() ; ++pt )
+	{
+		std::sort( this->values_of_landscapes[pt].begin() , this->values_of_landscapes[pt].end() , std::greater<double>() );	
+	}
+}//set_up_values_of_landscapes
+
+Persistence_landscape_on_grid::Persistence_landscape_on_grid( const std::vector< std::pair< double , double > >& p , double grid_min_ , double grid_max_ , size_t number_of_points_ )
+{
+	this->set_up_values_of_landscapes( p , grid_min_ , grid_max_ , number_of_points_ );
+}//Persistence_landscape_on_grid
+
+
+Persistence_landscape_on_grid::Persistence_landscape_on_grid(const 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 ) )
+	{
+		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::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);		
+		//std::cerr << "Reading line : " << line << std::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, double min_x , double max_x , double min_y , double max_y, size_t from_ , size_t to_  )const
+{		
+    //this program create a gnuplot script file that allows to plot persistence diagram.
+    std::ofstream out;
+
+    std::ostringstream nameSS;
+    nameSS << filename << "_GnuplotScript";
+    std::string nameStr = nameSS.str();
+    out.open( nameStr );
+
+	if ( min_x == max_x )
+	{
+		std::pair<double,double> min_max = compute_minimum_maximum();
+		out << "set xrange [" << this->grid_min << " : " << this->grid_max << "]" << std::endl;
+		out << "set yrange [" << min_max.first << " : " << min_max.second << "]" << std::endl;
+	}
+	else
+	{
+		out << "set xrange [" << min_x << " : " << max_x << "]" << std::endl;
+		out << "set yrange [" << min_y << " : " << max_y << "]" << std::endl;
+	}
+    
+    size_t number_of_nonzero_levels = this->number_of_nonzero_levels();
+    double dx = ( this->grid_max - this->grid_min )/((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 << std::endl;
+    }
+    
+    for ( size_t lambda = from ; lambda != to ; ++lambda )
+    {
+		double point = this->grid_min;
+        for ( size_t i = 0 ; i != this->values_of_landscapes.size() ; ++i )
+        {
+			double value = 0;
+			if ( this->values_of_landscapes[i].size() > lambda )			
+			{
+				value = this->values_of_landscapes[i][lambda];
+			}
+            out << point << " " << value << std::endl;
+            point += dx;
+        }
+        out << "EOF" << std::endl;
+    }
+    std::cout << "Gnuplot script to visualize persistence diagram written to the file: " << nameStr << ". Type load '" << nameStr << "' in gnuplot to visualize." << std::endl;
+}
+
+template < typename T > 
+Persistence_landscape_on_grid operation_on_pair_of_landscapes_on_grid ( const Persistence_landscape_on_grid& land1 ,  const Persistence_landscape_on_grid& land2 )
+{
+	//first we need to check if the domains are the same:
+	if ( !check_if_defined_on_the_same_domain(land1,land2) )throw "Two grids are not compatible"; 	
+	
+	T oper;
+	Persistence_landscape_on_grid result;
+	result.values_of_landscapes = std::vector< std::vector< double > >( land1.values_of_landscapes.size()  );
+	result.grid_min = land1.grid_min;
+	result.grid_max = land1.grid_max;
+	
+	//now we 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_distance_of_landscapes( const Persistence_landscape_on_grid& first, const Persistence_landscape_on_grid& second )
+{
+	double result = 0;
+	
+	//first we need to check if first and second is defined on the same domain"
+	if ( !check_if_defined_on_the_same_domain(first, second) )throw "Two grids are not compatible"; 
+	
+	for ( size_t i = 0 ; i != first.values_of_landscapes.size() ; ++i )
+	{
+		for ( size_t j = 0 ; j != std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) ; ++j )
+		{
+			if ( result < abs( first.values_of_landscapes[i][j] - second.values_of_landscapes[i][j] ) )
+			{
+				result = abs( first.values_of_landscapes[i][j] - second.values_of_landscapes[i][j] );
+			}
+		}
+		if ( first.values_of_landscapes[i].size() == std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) )
+		{
+			for ( size_t j = first.values_of_landscapes[i].size() ; j != second.values_of_landscapes[i].size() ; ++j )
+			{
+				if ( result < second.values_of_landscapes[i][j] )result = second.values_of_landscapes[i][j];
+			}
+		}
+		if ( second.values_of_landscapes[i].size() == std::min( first.values_of_landscapes[i].size() , second.values_of_landscapes[i].size() ) )
+		{
+			for ( size_t j = second.values_of_landscapes[i].size() ; j != first.values_of_landscapes[i].size() ; ++j )
+			{
+				if ( result < first.values_of_landscapes[i][j] )result = first.values_of_landscapes[i][j];
+			}
+		}
+	}  
+	return result;
+}
+
+
+
+}//namespace Gudhi_stat
+}//namespace Gudhi
+
+#endif
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/coo_kernel.h b/src/Gudhi_stat/include/gudhi/persistence_representations/coo_kernel.h
new file mode 100644
index 00000000..02c9527a
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/coo_kernel.h
@@ -0,0 +1,79 @@
+//this is an implementation of a kernel by Mathieu Carrière,
+//Steve Y. Oudot, and Maksim Ovsjanikov presended in Stable Topological
+//Signatures for Points on 3D Shapes. Proc. Sympos. on Geometry Processing, 2015
+
+#pragma once
+#include <cmath>
+#include <climits>
+
+#include "persistence_intervals.h"
+#include "kernel.h"
+#include "Distances_of_points_in_diagram.h"
+
+template <typename T , typename F>
+class coo_kernel : public Persistence_kernel<T>
+{
+public:
+    //Distances_of_points_in_diagram( std::vector< std::pair< T , T > > intervals , F f , size_t where_to_cut );
+    //Distances_of_points_in_diagram( char* filename , F f , size_t where_to_cut );
+
+    coo_kernel( std::vector< Persistence_intervals<T> >& intervals , F f , size_t where_to_cut = 100 );
+    coo_kernel( char* filename , F f , size_t where_to_cut = 100 );
+
+    //i tutaj pojawia sie problem, bo reprezentacja jest liczona wielokrotnie, strata czasu.
+    double compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode );
+private:
+    //in this case, we will construct the representation of diagrams in the vector space. That is why, we store here some more stuff compared to
+    //the Persistence_kernel class.
+    std::vector< Distances_of_points_in_diagram<T,F> > representation;
+    void constuct_representation( std::vector< Persistence_intervals<T> >& intervals , F f , size_t where_to_cut);
+};
+
+template <typename T , typename F>
+coo_kernel<T,F>::coo_kernel( std::vector< Persistence_intervals<T> >& intervals , F f , size_t where_to_cut ):Persistence_kernel<T>(intervals)
+{
+    this->constuct_representation(intervals, f , where_to_cut);
+}
+
+template <typename T , typename F>
+coo_kernel<T,F>::coo_kernel( char* filename , F f , size_t where_to_cut )
+{
+    cerr << "Constuct representation \n";
+    std::vector< Persistence_intervals<T> > intervals;
+    std::vector<std::string> names = readFileNames( filename );
+    for ( size_t file_no = 0 ; file_no != names.size() ; ++file_no )
+    {
+        cout << "Reading file : " << names[file_no] << endl;
+        Persistence_intervals<T> interval( (char*)names[file_no].c_str() );
+        intervals.push_back( interval );
+    }
+    this->constuct_representation(intervals, f , where_to_cut);
+}
+
+template <typename T , typename F>
+void coo_kernel<T,F>::constuct_representation( std::vector< Persistence_intervals<T> >& intervals , F f , size_t where_to_cut)
+{
+    std::vector< Distances_of_points_in_diagram<T,F> > representation( intervals.size() );
+    this->representation = representation;
+    #pragma omp parallel for
+    for ( size_t i = 0 ; i < intervals.size() ; ++i )
+    {
+        this->representation[i] = Distances_of_points_in_diagram<T,F>( intervals[i].intervals , f , where_to_cut );
+
+    }
+    this->number_of_intervals = intervals.size();
+}
+
+template <typename T , typename F>
+double coo_kernel<T,F>::compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode )
+{
+    double result = 0;
+    size_t range = std::min( this->representation[number_of_first_barcode].size() , this->representation[number_of_second_barcode].size() );
+    for ( size_t i = 0 ; i != range ; ++i )
+    {
+        result += this->representation[number_of_first_barcode].vector_in_position(i)*this->representation[number_of_second_barcode].vector_in_position(i);
+    }
+    return result;
+}
+
+
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/kernel.h b/src/Gudhi_stat/include/gudhi/persistence_representations/kernel.h
new file mode 100644
index 00000000..5c226303
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/kernel.h
@@ -0,0 +1,76 @@
+#pragma once
+#include "persistence_intervals.h"
+#include "read_files_names.h"
+
+template <typename T>
+class Persistence_kernel
+{
+public:
+    Persistence_kernel(){}
+    //Persistence_kernel( std::vector< Persistence_intervals<T> >& inter );
+    //Persistence_kernel( const char* filename );
+    Persistence_kernel(size_t number_of_intervals_):number_of_intervals(number_of_intervals_){}
+	std::vector< std::vector<double> > compute_Gram_matrix();
+	virtual double compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode )=0;
+	size_t& number(){return this->number_of_intervals;}
+protected:
+    //std::vector< Persistence_intervals<T> > intervals;
+    size_t number_of_intervals;
+};
+/*
+template <typename T>
+Persistence_kernel<T>::Persistence_kernel( std::vector< Persistence_intervals<T> >& inter )
+{
+    this->intervals(inter);
+}
+
+template <typename T>
+Persistence_kernel<T>::Persistence_kernel( const char* filename )
+{
+    std::vector<std::string> names = readFileNames( filename );
+    for ( size_t file_no = 0 ; file_no != names.size() ; ++file_no )
+    {
+        cout << "Reading file : " << names[file_no] << endl;
+        Persistence_intervals<double> interval( (char*)names[file_no].c_str() );
+        this->intervals.push_back( interval );
+    }
+}
+*/
+
+template <typename T>
+std::vector< std::vector<double> > Persistence_kernel<T>::compute_Gram_matrix()
+{
+    cerr << "compute_Gram_matrix \n";
+    bool dbg = true;
+    std::vector< std::vector<double> >  result( this->number_of_intervals );
+    for ( size_t i = 0 ; i != this->number_of_intervals ; ++i )
+    {
+        result[i] = std::vector<double>( this->number_of_intervals );
+    }
+    for ( size_t i = 0 ; i != this->number_of_intervals ; ++i )
+    {
+        cerr << "Computing row number : " << i << endl;
+        #pragma omp parallel for
+        for ( size_t j = i ; j < this->number_of_intervals ; ++j )
+        {
+            cerr << j << " ";
+            result[i][j] = result [j][i] = this->compute_scalar_product( i,j );
+        }
+    }
+    return result;
+}//compute_Gram_matrix
+
+void pint_Gram_matrix_to_file( std::vector< std::vector<double> >& matrix , char* filename )
+{
+    ofstream out;
+    out.open(filename);
+    for ( size_t i = 0 ; i != matrix.size() ; ++i )
+    {
+        for ( size_t j = 0 ; j != matrix[i].size() ; ++j )
+        {
+            out << matrix[i][j] << " ";
+        }
+        out << std::endl;
+    }
+    out.close();
+}
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/one_d_gaussians.h b/src/Gudhi_stat/include/gudhi/persistence_representations/one_d_gaussians.h
new file mode 100644
index 00000000..96a3215b
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/one_d_gaussians.h
@@ -0,0 +1,186 @@
+#pragma once
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <vector>
+#include <sstream>
+#include <iostream>
+
+
+
+#include "persistence_intervals.h"
+#include "Distances_of_points_in_diagram.h"
+#include "read_files_names.h"
+
+
+
+
+using namespace std;
+
+
+class normal_dist
+{
+public:
+    normal_dist(double mean , double stdev):stdev(stdev),mean(mean){}
+    double operator() (double x)
+    {
+        //compute value of the normal distribution at x:
+        return 1.0/(this->stdev*2.506628)*exp(-(x-this->mean)*(x-this->mean)/(2*this->stdev*this->stdev));
+    }
+private:
+    double mean;
+    double stdev;
+};
+
+
+template <typename T>
+class one_d_gaussians
+{
+public:
+    one_d_gaussians( std::vector< Persistence_intervals<T> > intervals , size_t number_of_points , double min_ = -1 , double max_ = -1 , size_t range = -1 );
+    void print_to_file( char* filename  );
+    size_t size()const{return this->gaussians.size();}
+    inline T operator [](size_t i)
+    {
+        if ( i >= this->gaussians.size() )throw("Index out of range! Operator [], one_d_gaussians class\n");
+        return this->gaussians[i].second;
+    }
+protected:
+    std::vector< std::pair<T,T> > gaussians;
+    size_t number_of_points;
+    T min_;
+    T max_;
+    size_t range;
+};
+
+
+template <typename T>
+one_d_gaussians<T>::one_d_gaussians( std::vector< Persistence_intervals<T> > intervals , size_t number_of_points , double min_ , double max_ , size_t range )
+{
+    bool dbg = true;
+
+    this->number_of_points = number_of_points;
+    //setting up min max and range:
+    if ( range == -1 )
+    {
+        this->range = number_of_points/10.00;
+    }
+    else
+    {
+        this->range = range;
+    }
+
+    if ( (min_ == -1) && (max_ == -1) )
+    {
+        if ( dbg )cerr << "Setting up min and max \n";
+        //in this case, we need to find min and max:
+        min_ = INT_MAX;
+        max_ = -INT_MAX;
+        for ( size_t i = 0 ; i != intervals.size() ; ++i )
+        {
+            for ( size_t j = 0 ; j != intervals[i].intervals.size() ; ++j )
+            {
+                if ( intervals[i].intervals[j].first < min_ )min_ = intervals[i].intervals[j].first;
+                if (intervals[i].intervals[j].first > max_ )max_ = intervals[i].intervals[j].first;
+            }
+        }
+    }
+
+    this->min_= min_;
+    this->max_ = max_;
+
+    if ( this->min_ == this->max_ )return;
+
+    if ( dbg )
+    {
+        cerr << "min_ : " << this->min_ << endl;
+        cerr << "max_ : " << this->max_ << endl;
+        cerr << "number_of_points : " << this->number_of_points << endl;
+        cerr << "range : " << this->range << endl;
+
+    }
+
+    //initialization of the output variable:
+    std::vector< std::pair<double,double> > result( this->number_of_points+1 );
+    T dx = (this->max_-this->min_)/(T)this->number_of_points;
+    T current = this->min_;
+    for ( size_t i = 0 ; i != this->number_of_points ; ++i )
+    {
+        result[i] = std::make_pair( current , 0 );
+        current += dx;
+    }
+
+    //for every diagram
+    for ( size_t diag_no = 0 ; diag_no != intervals.size() ; ++diag_no )
+    {
+        if ( dbg )
+        {
+            cerr << "We are considering diagram number : " << diag_no << endl;
+        }
+        //for every point in this diagram:
+        for ( size_t point_no = 0 ; point_no != intervals[diag_no].intervals.size() ; ++point_no )
+        {
+            if (dbg){cerr << "(intervals[diag_no].intervals[point_no].first-min_) : " << (intervals[diag_no].intervals[point_no].first-min_) << "\n";}
+
+            size_t position = ((intervals[diag_no].intervals[point_no].first-this->min_)/( this->max_ - this->min_ )) * this->number_of_points;
+            //now we need to draw a gaussian distribution 'range' value to left, and to right from this point
+
+            if ( dbg )
+            {
+                cerr << "point : (" << intervals[diag_no].intervals[point_no].first << "," << intervals[diag_no].intervals[point_no].second << "), its position : " << position << endl;
+            }
+            normal_dist d(intervals[diag_no].intervals[point_no].first,dx*this->range);
+            double length = ( intervals[diag_no].intervals[point_no].second - intervals[diag_no].intervals[point_no].first );
+
+            if ( position > this->range )
+            {
+                cerr << "Position of this element extends the current range. This point will not be taken into consideration \n";
+                continue;
+            }
+
+            for ( size_t i = 0 ; i != this->range ; ++i )
+            {
+                //cerr << "i : " << i << " range : " << this->range << endl;
+                if ( position+i < number_of_points )
+                {
+                    result[position+i].second += length*d( min_+(position+i)*dx );
+                    //cerr << "position+i : " << position+i << endl;
+                }
+                if ( position >= i )
+                {
+                    result[position-i].second += length*d( min_+(position-i)*dx );
+                    //cerr << "position-i : " << position-i << endl;
+                }
+            }
+
+            //if ( dbg )cerr << "Done \n";
+        }
+    }
+
+    //weighting, normalization. This is needed, othervise the result depends on the number of diagram involved, which should not be the case for averages.
+    for ( size_t i = 0 ; i != result.size() ; ++i )
+    {
+        //cerr << "result[i].second  : " << result[i].second  << endl;
+        std::pair<double,double> newResult = std::make_pair( result[i].first , result[i].second / (double)intervals.size() );
+        result[i] = newResult;
+        this->gaussians.push_back(newResult);
+        //cerr << "result[i].second  : " << result[i].second  << endl;getchar();
+    }
+
+}//plot_gaussians_centered_at_birth_times
+
+
+
+template <typename T>
+void one_d_gaussians<T>::print_to_file( char* filename  )
+{
+    cerr << "print_gaussians_centered_at_birth_times_to_file\n";
+    ofstream out;
+    out.open( filename );
+    for ( size_t i = 0 ; i != gaussians.size() ; ++i )
+    {
+        out << gaussians[i].first << " " << gaussians[i].second << endl;
+    }
+    out.close();
+}//print_to_file
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/one_d_gaussians_kernel.h b/src/Gudhi_stat/include/gudhi/persistence_representations/one_d_gaussians_kernel.h
new file mode 100644
index 00000000..671e7ad1
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/one_d_gaussians_kernel.h
@@ -0,0 +1,82 @@
+#pragma once
+
+#include <string>
+
+#include "kernel.h"
+#include "one_d_gaussians.h"
+
+
+template <typename T>
+class one_d_gaussians_kernel : public Persistence_kernel<T>
+{
+public:
+    one_d_gaussians_kernel( std::vector< Persistence_intervals<T> >&  intervals_ , size_t number_of_points , double min_ = -1 , double max_ = -1 , size_t range = -1 );
+    one_d_gaussians_kernel( char* filename , size_t number_of_points , double min_ = -1 , double max_ = -1 , size_t range = -1 );
+
+    double compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode );
+private:
+    std::vector< one_d_gaussians<T> > gaussians;
+    void construct( std::vector< Persistence_intervals<T> >&  intervals , size_t number_of_points , double min_ = -1 , double max_ = -1 , size_t range = -1 );
+};
+
+
+template <typename T>
+void one_d_gaussians_kernel<T>::construct( std::vector< Persistence_intervals<T> >&  intervals_ , size_t number_of_points , double min_ , double max_ , size_t range )
+{
+    for ( size_t nr = 0 ; nr != intervals_.size() ; ++nr )
+    {
+        std::vector< Persistence_intervals<T> > aa;
+        aa.push_back( intervals_[nr] );
+        cerr << "intervals_[nr].size() : " << intervals_[nr].size() << endl;
+        cerr << "A \n";
+        this->gaussians.push_back( one_d_gaussians<T>( aa , number_of_points , min_ , max_ , range ) );
+        cerr << "B \n";
+    }
+}
+
+template <typename T>
+one_d_gaussians_kernel<T>::one_d_gaussians_kernel( std::vector< Persistence_intervals<T> >& intervals , size_t number_of_points , double min_ , double max_ , size_t range ):Persistence_kernel<T>(intervals.size())
+{
+    this->construct( intervals , number_of_points , min_ , max_ , range );
+}
+
+template <typename T>
+one_d_gaussians_kernel<T>::one_d_gaussians_kernel( char* filename , size_t number_of_points , double min_ , double max_ , size_t range )
+{
+    std::vector<std::string> names = readFileNames( filename );
+    std::vector< Persistence_intervals<double> > intervals;
+    for ( size_t file_no = 0 ; file_no != names.size() ; ++file_no )
+    {
+        cout << "Reading file : " << names[file_no] << endl;
+        Persistence_intervals<double> interval( (char*)names[file_no].c_str() );
+        intervals.push_back( interval );
+    }
+
+    //in this case, if  min_ and max_ are not specified by the user, we should compute it, so that they are the same for all the considered cases:
+    if ( min_ == max_ == -1 )
+    {
+        min_ = INT_MAX;
+        max_ = INT_MIN;
+        //in this case, the values of min_ and max_ are not specified by the user, and we should specify them here.
+        for ( size_t diagNo = 0 ; diagNo != intervals.size() ; ++diagNo )
+        {
+            std::pair<T,T> mM = intervals[diagNo].min_max();
+            if ( min_ > mM.first )min_ = mM.first;
+            if ( max_ < mM.second )max_ = mM.second;
+        }
+    }
+    this->construct( intervals , number_of_points , min_ , max_ , range );
+    this->number_of_intervals = names.size();
+}
+
+template <typename T>
+double one_d_gaussians_kernel<T>::compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode )
+{
+    double result = 0;
+    for ( size_t i = 0 ; i != this->gaussians[number_of_first_barcode].size() ; ++i )
+    {
+        result += this->gaussians[number_of_first_barcode][i] * this->gaussians[number_of_second_barcode][i];
+    }
+    //cerr << "Scalar product : " << result << endl; getchar();
+    return result;
+}
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/pssk_old.h b/src/Gudhi_stat/include/gudhi/persistence_representations/pssk_old.h
new file mode 100644
index 00000000..48dee045
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/pssk_old.h
@@ -0,0 +1,101 @@
+//This is an implementation of a kernel described in //http://arxiv.org/pdf/1412.6821v1.pdf
+//by Reininghaus, Huber, Bauer and Kwitt.
+
+#pragma once
+#include <cmath>
+#include <climits>
+
+#include "persistence_intervals.h"
+#include "kernel.h"
+
+inline double norm_squared( double x , double y )
+{
+    return x*x + y*y;
+};
+
+template <typename T>
+class pssk_kernel : public Persistence_kernel<T>
+{
+public:
+    pssk_kernel( std::vector< Persistence_intervals<T> >& inter , double stdev_ ):Persistence_kernel<T>(inter.size()),stdev(stdev_)
+    {
+        this->intervals(inter);
+    }
+    pssk_kernel( const char* filename , double stdev_ = -1 )
+    {
+        //first, read the intervals:
+        std::vector<std::string> names = readFileNames( filename );
+        this->number_of_intervals = names.size();
+        std::vector< Persistence_intervals<T> > intervals;
+        this->intervals = intervals;
+        for ( size_t file_no = 0 ; file_no != names.size() ; ++file_no )
+        {
+            cout << "Reading file : " << names[file_no] << endl;
+            Persistence_intervals<T> interval( (char*)names[file_no].c_str() );
+            this->intervals.push_back(interval);
+            cerr << "Number of intervals : " << interval.size() << endl;
+
+        }
+        if ( stdev_ != -1 )
+        {
+            this->stdev = stdev_;
+        }
+        else
+        {
+            T min_ = INT_MAX;
+            T max_ = INT_MIN;
+            //compute range of the considered diagrams:
+            for ( size_t i = 0 ; i != this->intervals.size() ; ++i )
+            {
+                std::pair< T , T > min_max = this->intervals[i].min_max();
+                if ( min_ > min_max.first )min_ = min_max.first;
+                if ( max_ < min_max.second )max_ = min_max.second;
+            }
+            //and arbitrary choice based on min and max:
+            this->stdev = (max_ - min_)/100;
+        }
+
+    }
+    double compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode );
+private:
+    //since we do not construct explicit representation of the diagrams in the feature space, we do not store it here.
+    double stdev;
+    std::vector< Persistence_intervals<T> > intervals;
+};
+
+//First version, a naive implementation:
+template <typename T>
+double pssk_kernel<T>::compute_scalar_product( size_t number_of_first_barcode , size_t number_of_second_barcode )
+{
+    bool dbg = false;
+    //cerr << "Entering \n";getchar();
+    double result = 0;
+    for ( size_t i = 0 ; i != this->intervals[number_of_first_barcode].intervals.size() ; ++i )
+    {
+        //cerr << "i :" << i <<endl;
+        #pragma omp parallel for
+        for ( size_t j = 0 ; j != this->intervals[number_of_second_barcode].intervals.size() ; ++j )
+        {
+            if ( dbg )
+            {
+                cerr << "( " << this->intervals[number_of_first_barcode].intervals[i].first << "," << this->intervals[number_of_first_barcode].intervals[i].second << ") \n";
+                cerr << "( " << this->intervals[number_of_second_barcode].intervals[j].first << "," << this->intervals[number_of_second_barcode].intervals[j].second << ") \n";
+                getchar();
+            }
+            T toAdd =
+            exp( -( norm_squared( this->intervals[number_of_first_barcode].intervals[i].first-this->intervals[number_of_second_barcode].intervals[j].first, this->intervals[number_of_first_barcode].intervals[i].second- this->intervals[number_of_second_barcode].intervals[j].second ) ) / (8*this->stdev) )
+            +
+            exp( -( norm_squared( this->intervals[number_of_first_barcode].intervals[i].first-this->intervals[number_of_second_barcode].intervals[j].second, this->intervals[number_of_first_barcode].intervals[i].second- this->intervals[number_of_second_barcode].intervals[j].first ) ) / (8*this->stdev) );
+            result += toAdd;
+            if ( dbg )
+            {
+                cerr << "toAdd : " << toAdd << endl;
+                cerr << result << " : " << result << endl;
+                getchar();
+            }
+        }
+    }
+    result *= 1/( 8*this->stdev*3.141592 );
+    //cerr << "result : " << result << endl; getchar();
+    return result;
+}
diff --git a/src/Gudhi_stat/include/gudhi/persistence_representations/temp b/src/Gudhi_stat/include/gudhi/persistence_representations/temp
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/src/Gudhi_stat/include/gudhi/persistence_representations/temp