Adding tests to the implementation of persitence heat maps.

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

Former-commit-id: 065fdfed0c7000ff00bc4b1752239fedb2cdc9e4

3 files changed, 432 insertions, 112 deletions

diff --git a/src/Gudhi_stat/example/persistence_heat_maps.cpp b/src/Gudhi_stat/example/persistence_heat_maps.cpp
index 475569c8..fc9e31a0 100644
--- a/src/Gudhi_stat/example/persistence_heat_maps.cpp
+++ b/src/Gudhi_stat/example/persistence_heat_maps.cpp
@@ -89,10 +89,98 @@ int main( int argc , char** argv )
 	q.load_from_file( "persistence_heat_map_from_file_with_diagram" );	
 	
 	cerr << (p == q) << endl;
+*/
+/*
+	//test of computations of a mean:
+	std::vector< std::pair< double,double > > intervals;
+	intervals.push_back( std::make_pair(5,5) );	
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1);		
+	Persistence_heat_maps p( intervals , filter ,  constant_function, false , 100 , 0 , 10 );
+	p.plot( "heat_map_1" );
+
+	
+	std::vector< std::pair< double,double > > intervals2;
+	intervals2.push_back( std::make_pair(7,7) );		
+	Persistence_heat_maps q( intervals2 , filter ,  constant_function, false , 100 , 0 , 10 );
+	q.plot( "heat_map_2" );
+	
+	
+	std::vector< Abs_Topological_data_with_averages* > to_average;
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&p) );
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&q) );
+	Persistence_heat_maps av;	
+	av.compute_average( to_average );
+	av.plot( "average" );
+*/
+
+/*
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 10 );
+	Persistence_heat_maps q( "file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 10 );
+	Persistence_heat_maps r( "file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 10 );
+	
+	std::vector< Abs_Topological_data_with_averages* > to_average;
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&p) );
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&q) );
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&r) );
+	Persistence_heat_maps av;	
+	av.compute_average( to_average );
+	
+	av.write_to_file( "template_average_of_heat_maps" );
 */	
+
+/*
+	std::vector< std::pair< double,double > > intervals;
+	intervals.push_back( std::make_pair(5,5) );	
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1);		
+	Persistence_heat_maps p( intervals , filter ,  constant_function, false , 10 , 0 , 10 );
+	p.plot( "heat_map_1" );
+	
+	std::vector< std::pair< double,double > > intervals2;
+	intervals2.push_back( std::make_pair(7,7) );		
+	Persistence_heat_maps q( intervals2 , filter ,  constant_function, false , 10 , 0 , 10 );
+	q.plot( "heat_map_2" );
+	
+	std::vector< Persistence_heat_maps* > to_compute_median;
+	to_compute_median.push_back( &p );
+	to_compute_median.push_back( &q );	
+	Persistence_heat_maps median;	
+	median.compute_median( to_compute_median );
+	median.plot( "median" );
+*/
+
+/*
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
 	
+	std::vector< Persistence_heat_maps* > to_compute_median;
+	to_compute_median.push_back( &p );
+	to_compute_median.push_back( &q );
+	to_compute_median.push_back( &r );
+	Persistence_heat_maps median;	
+	median.compute_median( to_compute_median );
 	
+	median.write_to_file( "template_median_of_heat_maps" );
+*/
+
 
+
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	
+	std::vector< Persistence_heat_maps* > to_compute_percentage_of_active;
+	to_compute_percentage_of_active.push_back( &p );
+	to_compute_percentage_of_active.push_back( &q );
+	to_compute_percentage_of_active.push_back( &r );
+	Persistence_heat_maps percentage_of_active;	
+	percentage_of_active.compute_percentage_of_active( to_compute_percentage_of_active , 0.1 );
+	
+	percentage_of_active.write_to_file( "template_percentage_of_active_of_heat_maps" );
+	//percentage_of_active.plot( "template_percentage_of_active_of_heat_maps" );
 	
 	
 
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
index f6b17915..0910ae71 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
@@ -189,24 +189,62 @@ public:
     **/ 
     inline bool check_if_the_same( const Persistence_heat_maps& second  )
     {
-		if ( this->heat_map.size() != second.heat_map.size() )return false;
-		if ( this->min_ != second.min_ )return false;
-		if ( this->max_ != second.max_ )return false;
+		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 gimme_min(){return this->min_;}
+
+	/**
+	 * Return maximal range value of persistent image.
+	**/ 
+	inline double gimme_max(){return this->max_;}
+	
+	
 	/**
 	 * Operator == to check if to persistence heat maps are the same.
 	**/ 
 	bool operator == ( const Persistence_heat_maps& rhs )
 	{
-		if ( !this->check_if_the_same(rhs) ) return false;//in this case, the domains are not the same, so the maps cannot be the same.
+		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]) )return false;
+				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;
@@ -228,99 +266,11 @@ public:
     
     
     //Concretizations of virtual methods:    
-    std::vector<double> vectorize( int number_of_function )
-    {
-		//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;
-	}
-	
-	double distance( const Abs_Topological_data_with_distances* second_ , double power = 1)
-	{		
-		Persistence_heat_maps* second = (Persistence_heat_maps*)second_;
-		
-		//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;
-		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( abs(this->heat_map[i][j] - second->heat_map[i][j]) , power );
-			}
-		}
-		return distance;
-	}
-	
-	double project_to_R( int number_of_function )
-	{
-		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;
-	}
-	
-	void compute_average( std::vector< Abs_Topological_data_with_averages* > to_average )
-	{
-		std::vector<Persistence_heat_maps*> maps;
-		maps.reserve( to_average.size() );
-		for ( size_t i = 0 ; i != to_average.size() ; ++i )
-		{
-			maps.push_back( (Persistence_heat_maps*)to_average[i] );
-		}		
-		this->compute_mean( maps );
-	}
-	
-	double compute_scalar_product( const Abs_Topological_data_with_scalar_product* second_ )
-	{
-		Persistence_heat_maps* second = (Persistence_heat_maps*)second_;
-		
-		//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;
-	}
+    std::vector<double> vectorize( int number_of_function );
+	double distance( const Abs_Topological_data_with_distances* second_ , double power = 1);
+	double project_to_R( int number_of_function );	
+	void compute_average( std::vector< Abs_Topological_data_with_averages* > to_average );	
+	double compute_scalar_product( const Abs_Topological_data_with_scalar_product* second_ );
 	
 private:
 	//private methods
@@ -352,7 +302,7 @@ void Persistence_heat_maps::construct( const std::vector< std::pair<double,doubl
 									   bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ )
 {
     bool dbg = false;       
-    if ( dbg )std::cerr << "Eneterring construct procedure \n";
+    if ( dbg )std::cerr << "Entering construct procedure \n";
 
     if ( min_ == max_ )
     {
@@ -504,9 +454,9 @@ void Persistence_heat_maps::compute_median( const std::vector<Persistence_heat_m
 	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_map.size() ; ++i )
+    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
     {
-        for ( size_t j = 0 ; j != heat_map[i].size() ; ++j )
+        for ( size_t j = 0 ; j != heat_maps[i].size() ; ++j )
         {
 			for ( size_t map_no = 0 ; map_no != maps.size() ; ++map_no )
 			{
@@ -517,6 +467,8 @@ void Persistence_heat_maps::compute_median( const std::vector<Persistence_heat_m
         }       
     }
 	this->heat_map = heat_maps;
+	this->min_= maps[0]->min_;
+	this->max_= maps[0]->max_;
 }
 
 
@@ -524,34 +476,33 @@ void Persistence_heat_maps::compute_median( const std::vector<Persistence_heat_m
 void Persistence_heat_maps::compute_mean( 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_map.size() ; ++i )
+		    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
     {
-        for ( size_t j = 0 ; j != heat_map[i].size() ; ++j )
+        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_;
 }
 
 
 
 
-
 void Persistence_heat_maps::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 != this->heat_map.size() ; ++i )
+    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
     {
-        for ( size_t j = 0 ; j != this->heat_map[i].size() ; ++j )
+        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 )
@@ -569,6 +520,8 @@ void Persistence_heat_maps::compute_percentage_of_active( const std::vector<Pers
         }       
     }
     this->heat_map = heat_maps;
+    this->min_ = maps[0]->min_;
+	this->max_ = maps[0]->max_;
 }
 
 
@@ -669,6 +622,103 @@ void Persistence_heat_maps::load_from_file( const char* filename )
 }
 
 
+//Concretizations of virtual methods:    
+std::vector<double> Persistence_heat_maps::vectorize( int number_of_function )
+{
+	//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;
+}
+
+double Persistence_heat_maps::distance( const Abs_Topological_data_with_distances* second_ , double power )
+{		
+	Persistence_heat_maps* second = (Persistence_heat_maps*)second_;
+	
+	//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;
+	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( abs(this->heat_map[i][j] - second->heat_map[i][j]) , power );
+		}
+	}
+	return distance;
+}
+
+double Persistence_heat_maps::project_to_R( int number_of_function )
+{
+	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;
+}
+
+void Persistence_heat_maps::compute_average( std::vector< Abs_Topological_data_with_averages* > to_average )
+{
+	std::vector<Persistence_heat_maps*> maps;
+	maps.reserve( to_average.size() );
+	for ( size_t i = 0 ; i != to_average.size() ; ++i )
+	{
+		maps.push_back( (Persistence_heat_maps*)to_average[i] );
+	}					
+	this->compute_mean( maps );
+}
+
+double Persistence_heat_maps::compute_scalar_product( const Abs_Topological_data_with_scalar_product* second_ )
+{
+	Persistence_heat_maps* second = (Persistence_heat_maps*)second_;
+	
+	//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
diff --git a/src/Gudhi_stat/test/persistence_heat_maps_test.cpp b/src/Gudhi_stat/test/persistence_heat_maps_test.cpp
index c05aff35..78297cea 100644
--- a/src/Gudhi_stat/test/persistence_heat_maps_test.cpp
+++ b/src/Gudhi_stat/test/persistence_heat_maps_test.cpp
@@ -41,7 +41,7 @@ using namespace std;
 
 	
 
-BOOST_AUTO_TEST_CASE(check_construction_of_landscape) 
+BOOST_AUTO_TEST_CASE(check_construction_of_heat_maps) 
 {	
 	std::vector< std::vector<double> > filter = create_Gaussian_filter(100,1);		
 	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
@@ -49,8 +49,190 @@ BOOST_AUTO_TEST_CASE(check_construction_of_landscape)
 	
 	Persistence_heat_maps q;
 	q.load_from_file( "data/persistence_heat_map_from_file_with_diagram" );	
+		
+	BOOST_CHECK( p == q );
+}
+
+
+BOOST_AUTO_TEST_CASE(check_averages_of_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 10 );
+	Persistence_heat_maps q( "data/file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 10 );
+	Persistence_heat_maps r( "data/file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 10 );
 	
+	std::vector< Abs_Topological_data_with_averages* > to_average;
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&p) );
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&q) );
+	to_average.push_back( (Abs_Topological_data_with_averages*)(&r) );
+	Persistence_heat_maps av;	
+	av.compute_average( to_average );
 	
+	Persistence_heat_maps template_average;
+	template_average.load_from_file( "data/template_average_of_heat_maps" );
+	
+	BOOST_CHECK( av == template_average );
+}
+
+
+BOOST_AUTO_TEST_CASE(check_median_of_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "data/file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "data/file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	
+	std::vector< Persistence_heat_maps* > to_compute_median;
+	to_compute_median.push_back( &p );
+	to_compute_median.push_back( &q );
+	to_compute_median.push_back( &r );
+	Persistence_heat_maps median;	
+	median.compute_median( to_compute_median );
+	
+	Persistence_heat_maps template_median;
+	template_median.load_from_file( "data/template_median_of_heat_maps" );
+	
+	BOOST_CHECK( median == template_median );	
+}
+
+
+BOOST_AUTO_TEST_CASE(check_compute_percentage_of_active_of_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "data/file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "data/file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	
+	std::vector< Persistence_heat_maps* > to_compute_percentage_of_active;
+	to_compute_percentage_of_active.push_back( &p );
+	to_compute_percentage_of_active.push_back( &q );
+	to_compute_percentage_of_active.push_back( &r );
+	Persistence_heat_maps percentage_of_active;	
+	percentage_of_active.compute_percentage_of_active( to_compute_percentage_of_active , 0.1 );
+	
+	Persistence_heat_maps template_percentage_of_active;
+	template_percentage_of_active.load_from_file( "data/template_percentage_of_active_of_heat_maps" );
+	
+	BOOST_CHECK( percentage_of_active == template_percentage_of_active );	
+}
+
+BOOST_AUTO_TEST_CASE(check_vectorize_for_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 5 , 0 , 1 );
+	
+	std::vector<double> p_vect_template;
+	
+	p_vect_template.push_back(0.0606728);
+	p_vect_template.push_back(0.0610023); 
+	p_vect_template.push_back(0.0607978 );
+	p_vect_template.push_back(0.0600647 );
+	p_vect_template.push_back(0.0588224 );
+	p_vect_template.push_back(0.0619829 );
+	p_vect_template.push_back(0.0623218 );
+	p_vect_template.push_back(0.0621152 );
+	p_vect_template.push_back(0.0613686 );
+	p_vect_template.push_back(0.0601016 );
+	p_vect_template.push_back(0.0627679 );
+	p_vect_template.push_back(0.0631134 );
+	p_vect_template.push_back(0.0629066 );
+	p_vect_template.push_back(0.0621528 );
+	p_vect_template.push_back(0.0608719 );
+	p_vect_template.push_back(0.0630073 );
+	p_vect_template.push_back(0.0633564 );
+	p_vect_template.push_back(0.0631511 );
+	p_vect_template.push_back(0.0623968 );
+	p_vect_template.push_back(0.0611132 );
+	p_vect_template.push_back(0.0626947 );
+	p_vect_template.push_back(0.0630445 );
+	p_vect_template.push_back(0.0628425 );
+	p_vect_template.push_back(0.0620941 );
+	p_vect_template.push_back(0.060819);
+	
+	std::vector<double> p_vect = p.vectorize(0);
+	for ( size_t i = 0 ; i != p_vect.size() ; ++i )
+	{
+		BOOST_CHECK( almost_equal( p_vect_template[i] , p_vect[i] ) );
+	}	
+}
+
+
+BOOST_AUTO_TEST_CASE(check_distance_for_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "data/file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "data/file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
+
+	//cerr << p.distance( (Abs_Topological_data_with_distances*)&p ) << endl;
+	//cerr << p.distance( (Abs_Topological_data_with_distances*)&q ) << endl;
+	//cerr << p.distance( (Abs_Topological_data_with_distances*)&r ) << endl;	
+	//cerr << q.distance( (Abs_Topological_data_with_distances*)&p ) << endl;
+	//cerr << q.distance( (Abs_Topological_data_with_distances*)&q ) << endl;
+	//cerr << q.distance( (Abs_Topological_data_with_distances*)&r ) << endl;	
+	//cerr << r.distance( (Abs_Topological_data_with_distances*)&p ) << endl;
+	//cerr << r.distance( (Abs_Topological_data_with_distances*)&q ) << endl;
+	//cerr << r.distance( (Abs_Topological_data_with_distances*)&r ) << endl;	
+	//0         624.183   415.815
+	//624.183   0         528.06Z
+	//415.815   528.066   0
+
+	BOOST_CHECK( fabs( p.distance( (Abs_Topological_data_with_distances*)&p ) - 0) < 0.0005);
+	BOOST_CHECK( fabs( p.distance( (Abs_Topological_data_with_distances*)&q ) - 624.183)< 0.0005);
+	BOOST_CHECK( fabs( p.distance( (Abs_Topological_data_with_distances*)&r ) - 415.815)< 0.0005);
+	BOOST_CHECK( fabs( q.distance( (Abs_Topological_data_with_distances*)&p ) - 624.183)< 0.0005);
+	BOOST_CHECK( fabs( q.distance( (Abs_Topological_data_with_distances*)&q ) - 0)< 0.0005);
+	BOOST_CHECK( fabs( q.distance( (Abs_Topological_data_with_distances*)&r ) - 528.066)< 0.0005);
+	BOOST_CHECK( fabs( r.distance( (Abs_Topological_data_with_distances*)&p ) - 415.815)< 0.0005);
+	BOOST_CHECK( fabs( r.distance( (Abs_Topological_data_with_distances*)&q ) - 528.066)< 0.0005);
+	BOOST_CHECK( fabs( r.distance( (Abs_Topological_data_with_distances*)&r ) - 0)< 0.0005);
+}
+
+
+
+BOOST_AUTO_TEST_CASE(check_projections_to_R_for_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "data/file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "data/file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
+
+	//cerr << p.project_to_R(0) << endl;
+	//cerr << q.project_to_R(0) << endl;
+	//cerr << r.project_to_R(0) << endl;
+
+	BOOST_CHECK( fabs( p.project_to_R(0) - 44.3308 )< 0.0005);
+	BOOST_CHECK( fabs( q.project_to_R(0) - 650.568 )< 0.0005);
+	BOOST_CHECK( fabs( r.project_to_R(0) - 429.287 )< 0.0005);
+}
+
+
+BOOST_AUTO_TEST_CASE(check_scalar_products_for_heat_maps) 
+{
+	std::vector< std::vector<double> > filter = create_Gaussian_filter(30,1);	
+	Persistence_heat_maps p( "data/file_with_diagram" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps q( "data/file_with_diagram_1" , filter ,  constant_function, false , 1000 , 0 , 1 );
+	Persistence_heat_maps r( "data/file_with_diagram_2" , filter ,  constant_function, false , 1000 , 0 , 1 );
+
+	//cerr << p.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&p ) << endl;
+	//cerr << p.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&q ) << endl;
+	//cerr << p.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&r ) << endl;
+	//cerr << q.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&p ) << endl;
+	//cerr << q.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&q ) << endl;
+	//cerr << q.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&r ) << endl;	
+	//cerr << r.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&p ) << endl;
+	//cerr << r.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&q ) << endl;
+	//cerr << r.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&r ) << endl;
+	
+
+	BOOST_CHECK( fabs( p.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&p ) - 0.0345687 )< 0.0005);
+	BOOST_CHECK( fabs( p.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&q ) - 0.0509357 )< 0.0005);
+	BOOST_CHECK( fabs( p.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&r ) - 0.0375608 )< 0.0005);
+	BOOST_CHECK( fabs( q.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&p ) - 0.0509357 )< 0.0005);
+	BOOST_CHECK( fabs( q.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&q ) - 1.31293 )< 0.0005);
+	BOOST_CHECK( fabs( q.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&r ) - 0.536799)< 0.0005);
+	BOOST_CHECK( fabs( r.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&p ) - 0.0375608)< 0.0005);
+	BOOST_CHECK( fabs( r.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&q ) - 0.536799)< 0.0005);
+	BOOST_CHECK( fabs( r.compute_scalar_product( (Abs_Topological_data_with_scalar_product*)&r ) - 0.672907)< 0.0005);
 	
-	BOOST_CHECK( p == q );
 }