/*    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) 2017  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_H_
#define PERSISTENCE_LANDSCAPE_H_

//standard include
#include <cmath>
#include <iostream>
#include <vector>
#include <limits>
#include <fstream>
#include <sstream>
#include <algorithm>


//gudhi include
#include <gudhi/read_persistence_from_file.h>
#include <gudhi/common_persistence_representations.h>




namespace Gudhi
{
namespace Persistence_representations
{



// pre declaration
class Persistence_landscape;
template < typename operation > 
Persistence_landscape operation_on_pair_of_landscapes( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 );



/**
 * \class Persistence_landscape Persistence_landscape.h gudhi/Persistence_landscape.h
 * \brief A class implementing persistence landscapes data structures.
 *
 * \ingroup Persistence_representations
 *
 * \details
 * For theoretical description, please consult <i>Statistical topological data analysis using persistence
 * landscapes</i>\cite bubenik_landscapes_2015 , and for details of algorithms,
 * <i>A persistence landscapes toolbox for topological statistics</i>\cite bubenik_dlotko_landscapes_2016.
 *
 * Persistence landscapes allow vectorization, 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
 *
 * Note that at the moment, due to rounding errors during the construction of persistence landscapes, elements which
 * are different by 0.000005 are considered the same. If the scale in your persistence diagrams is comparable to this
 * value, please rescale them before use this code.
 *
**/
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 = std::numeric_limits<unsigned>::max() );



    /**
     * 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 );
    
    template < typename operation > 
	friend Persistence_landscape operation_on_pair_of_landscapes( const Persistence_landscape& land1 ,  const Persistence_landscape& land2 );




	/**
	 *\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 multiplication 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 multiplication 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 \f$L^i\f$ 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 \f$L^{\infty}\f$ distance between two landscapes.
	**/
    friend double compute_max_norm_distance_of_landscapes( const Persistence_landscape& first, const Persistence_landscape& second );


	/**
	 *\private Computations of \f$L^{p}\f$ 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 between 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(); }

	/**
	 *  Compute 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.
         * At the moment this function is not tested, since it is quite likely to be changed in the future. Given this, when using it, keep in mind that it
         * will be most likely changed in the next versions.
	**/
    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 landscape. 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 number of intermediate complexes. We have to clean that up, but we cannot erase the initial landscapes 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 using 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 function is 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 function is 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 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 );    
    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)
{    
    std::vector< std::pair< double , double > > barcode;    
    if ( dimension < std::numeric_limits<double>::max() )
    {
		barcode = read_persistence_intervals_in_one_dimension_from_file( filename , dimension );	
	}
    else
    {
		barcode = read_persistence_intervals_in_one_dimension_from_file( filename );
    }    
    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) ) )
            {
                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 up 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 integral.
		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 up 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;
            }
        }
    }
    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 o result
            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 before head
        std::ifstream in;
        in.open( filename );
	if ( !in.good() )
	{
		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 point : " << 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";
            }

            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";
                }
                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 ) // In this case, I consider points at the infinity
        {
            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 Persistence_representations
}  // namespace gudhi


#endif  // PERSISTENCE_LANDSCAPE_H_