A few more modifications to the Gudhi stat. Now all the utilities are tested, and some corrections to the tests are added.

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

Former-commit-id: 5032a8c66a13764ee929a76ac301667f4d67de45

8 files changed, 382 insertions, 276 deletions

diff --git a/src/Gudhi_stat/example/utilities/persistence_vectors/average_persistence_vectors.cpp b/src/Gudhi_stat/example/utilities/persistence_vectors/average_persistence_vectors.cpp
index f0c0ba61..1b269e9c 100644
--- a/src/Gudhi_stat/example/utilities/persistence_vectors/average_persistence_vectors.cpp
+++ b/src/Gudhi_stat/example/utilities/persistence_vectors/average_persistence_vectors.cpp
@@ -35,7 +35,8 @@ using namespace Gudhi::Gudhi_stat;
 
 int main( int argc , char** argv )
 {
-	std::cout << "This program computes average persistence landscape of persistence landscapes created based on persistence diagrams provided as an input. Please call this program with the names of files with persistence diagrams \n";
+	std::cout << "This program computes average persistence vector of persistence vectors created based on persistence diagrams provided as an input. \n";
+	std::cout << "Please call this program with the names of files with persistence diagrams \n";
 	std::vector< const char* > filenames;
 	
 	if ( argc == 1 )
@@ -49,10 +50,10 @@ int main( int argc , char** argv )
 		filenames.push_back( argv[i] );
 	}
 	
-	std::cout << "Creating persistence landscapes...\n";
+	std::cout << "Reading persistence vectors...\n";
 	std::vector< Abs_Topological_data_with_averages* > lands;
 	for ( size_t i = 0 ; i != filenames.size() ; ++i )
-	{
+	{		
 		Vector_distances_in_diagram< euclidean_distance<double> >* l = new Vector_distances_in_diagram< euclidean_distance<double> >;
 		l->load_from_file( filenames[i] );
 		lands.push_back( (Abs_Topological_data_with_averages*)l );
@@ -61,7 +62,7 @@ int main( int argc , char** argv )
 	Vector_distances_in_diagram< euclidean_distance<double> > av;
 	av.compute_average( lands );
 	
-	av.print_to_file( "average" );
+	av.print_to_file( "average.vect" );
 	
 	for ( size_t i = 0 ; i != filenames.size() ; ++i )
 	{
diff --git a/src/Gudhi_stat/example/utilities/persistence_vectors/compute_distance_of_persistence_vectors.cpp b/src/Gudhi_stat/example/utilities/persistence_vectors/compute_distance_of_persistence_vectors.cpp
index f93c910a..72026497 100644
--- a/src/Gudhi_stat/example/utilities/persistence_vectors/compute_distance_of_persistence_vectors.cpp
+++ b/src/Gudhi_stat/example/utilities/persistence_vectors/compute_distance_of_persistence_vectors.cpp
@@ -36,9 +36,9 @@ using namespace Gudhi::Gudhi_stat;
 
 int main( int argc , char** argv )
 {
-	std::cout << "This program compute dsitance of persistence vectors stored in a file (the file needs to be created beforehand). \n";	
-	std::cout << "The first parameter of a program is an interger p. The program compute l^p distance of the vectors. For L^infty distance choose p = -1. \n";
-	std::cout << "The remaining parameters of this programs are names of files with persistence vectors.";
+	std::cout << "This program compute distance of persistence vectors stored in a file (the file needs to be created beforehand). \n";	
+	std::cout << "The first parameter of a program is an interger p. The program compute l^p distance of the vectors. For l^infty distance choose p = -1. \n";
+	std::cout << "The remaining parameters of this programs are names of files with persistence vectors.\n";
 	
 	if ( argc < 3 )
 	{
@@ -46,7 +46,7 @@ int main( int argc , char** argv )
 		return 1;
 	}
 	
-	int p = atoi( argv[3] );
+	int p = atoi( argv[1] );
 
 	std::vector< const char* > filenames;
 	for ( int i = 2 ; i < argc ; ++i )
@@ -57,6 +57,7 @@ int main( int argc , char** argv )
 	vectors.reserve( filenames.size() );
 	for ( size_t file_no = 0 ; file_no != filenames.size() ; ++file_no )
 	{
+		//cerr << filenames[file_no] << endl;
 		Vector_distances_in_diagram< euclidean_distance<double> >* l = new Vector_distances_in_diagram< euclidean_distance<double> >;
 		l->load_from_file( filenames[file_no] );
 		vectors.push_back( l );
@@ -71,11 +72,11 @@ int main( int argc , char** argv )
 		std::vector< double > v( filenames.size() , 0 );
 		distance[i] = v;
 	}
-	
+		
 	//and now we can compute the distances:
 	for ( size_t i = 0 ; i != vectors.size() ; ++i )
 	{
-		for ( size_t j = i ; j != vectors.size() ; ++j )
+		for ( size_t j = i+1 ; j != vectors.size() ; ++j )
 		{
 			distance[i][j] = distance[j][i] = ((Vector_distances_in_diagram< euclidean_distance<double> >*)vectors[i])->distance( vectors[j] , p ) ;
 		}
diff --git a/src/Gudhi_stat/example/utilities/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp b/src/Gudhi_stat/example/utilities/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
index c886d076..16192d52 100644
--- a/src/Gudhi_stat/example/utilities/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
+++ b/src/Gudhi_stat/example/utilities/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
@@ -37,7 +37,7 @@ using namespace Gudhi::Gudhi_stat;
 int main( int argc , char** argv )
 {
 	std::cout << "This program compute scalar product of persistence vectors stored in a file (the file needs to be created beforehand). \n";	
-	std::cout << "The parameters of this programs are names of files with persistence vectors.";
+	std::cout << "The parameters of this programs are names of files with persistence vectors.\n";
 
 	std::vector< const char* > filenames;
 	for ( int i = 1 ; i < argc ; ++i )
diff --git a/src/Gudhi_stat/example/utilities/persistence_vectors/create_persistence_vectors.cpp b/src/Gudhi_stat/example/utilities/persistence_vectors/create_persistence_vectors.cpp
index b800c872..99015e41 100644
--- a/src/Gudhi_stat/example/utilities/persistence_vectors/create_persistence_vectors.cpp
+++ b/src/Gudhi_stat/example/utilities/persistence_vectors/create_persistence_vectors.cpp
@@ -43,10 +43,10 @@ int main( int argc , char** argv )
 		filenames.push_back( argv[i] );
 	}
 	
-	std::cout << "Creating persistence vectors...\n";	
 	for ( size_t i = 0 ; i != filenames.size() ; ++i )
 	{
-		Vector_distances_in_diagram< euclidean_distance<double> > l( filenames[i] , -1 );
+		std::cerr << "Creatign persistence vectors based on a file : " << filenames[i] << std::endl;
+		Vector_distances_in_diagram< euclidean_distance<double> > l( filenames[i] , -1 );				
 		std::stringstream ss;
 		ss << filenames[i] << ".vect";
 		l.print_to_file( ss.str().c_str() );
@@ -54,3 +54,4 @@ int main( int argc , char** argv )
 	std::cout << "Done \n";
 	return 0;
 }
+
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 e758b30a..6d06ac8e 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
@@ -88,7 +88,7 @@ std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius ,
             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)
@@ -99,6 +99,19 @@ std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius ,
 		}
             
     }
+    
+    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;
 }
 
@@ -112,7 +125,7 @@ std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius ,
 **/ 
 
 double constant_function( const std::pair< double , double >& point_in_diagram )
-{
+{	
 	return 1;
 }
 
@@ -179,7 +192,7 @@ public:
 	 * (6) a min x and y value of points that are to be taken into account. By default it is set to -1, 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 -1, in which case the program compute the values based on the data.
 	**/ 	
-    Persistence_heat_maps( const char* name_of_file_with_names_of_files_with_interval , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram ) = constant_function, bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  );
+    Persistence_heat_maps( const char* filename , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram ) = constant_function, bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  );
 
 	
 	/**
@@ -362,6 +375,7 @@ void Persistence_heat_maps::construct( const std::vector< std::pair<double,doubl
 
     this->min_ = min_;
     this->max_ = max_;    
+    this->scalling_function_with_respect_to_distance_from_diagonal = scalling_function_with_respect_to_distance_from_diagonal;
 
 
     //initialization of the structure heat_map
@@ -393,12 +407,20 @@ void Persistence_heat_maps::construct( const std::vector< std::pair<double,doubl
 		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 << endl;
+			std::cerr << "y_grid : " << y_grid << endl;
+		}
+				
 		double scaling_value = this->scalling_function_with_respect_to_distance_from_diagonal(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()) 
@@ -406,6 +428,7 @@ void Persistence_heat_maps::construct( const std::vector< std::pair<double,doubl
 					  ((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 )
 					{
@@ -441,12 +464,12 @@ Persistence_heat_maps::Persistence_heat_maps( const std::vector< std::pair< doub
 }
 
 
-Persistence_heat_maps::Persistence_heat_maps( const char* name_of_file_with_names_of_files_with_interval , 
+Persistence_heat_maps::Persistence_heat_maps( const char* filename , 
 											  std::vector< std::vector<double> > filter, 
 											  double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram ), 
 											  bool erase_below_diagonal , size_t number_of_pixels  , double min_ , double max_ )
 {    
-    std::vector< std::pair< double , double > > intervals_ = read_standard_file( name_of_file_with_names_of_files_with_interval );   
+    std::vector< std::pair< double , double > > intervals_ = read_standard_file( filename );           
     this->construct( intervals_ ,  filter ,  constant_function, erase_below_diagonal , number_of_pixels , min_ , max_ );
     this->set_up_parameters_for_basic_classes();
 }
@@ -559,7 +582,10 @@ void Persistence_heat_maps::compute_percentage_of_active( const std::vector<Pers
 void Persistence_heat_maps::plot( const char* filename )
 {
 	ofstream out;
-    out.open( filename );
+	std::stringstream ss;
+	ss << filename << "_Gnuplot_script";
+	
+    out.open( ss.str().c_str() );
 	out << "plot      '-' matrix with image" << std::endl;
     for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
     {
@@ -569,14 +595,15 @@ void Persistence_heat_maps::plot( const char* filename )
         }
         out << 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; 
 }
 
 
 
 void Persistence_heat_maps::print_to_file( const char* filename )
 {
+
 	ofstream out;
 	out.open( filename );
 	
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h b/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h
index 56fad188..8c3d0f37 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Vector_distances_in_diagram.h
@@ -38,6 +38,7 @@
 #include <gudhi/abstract_classes/Abs_Real_valued_topological_data.h>
 #include <gudhi/abstract_classes/Abs_Topological_data_with_scalar_product.h>
 #include <gudhi/concretizations/read_persitence_from_file.h>
+#include <gudhi/common.h>
 using namespace std;
 
 namespace Gudhi 
@@ -137,19 +138,37 @@ public:
 	**/
     void write_to_file( const char* filename );
     
-    /**
+     /**
 	 * Write a vector to a file.
 	**/
     void print_to_file( const char* filename )
     {
-		this->write_to_file(filename);
-	}
+		this->write_to_file( filename );
+	}  
     
     /**
 	 * Loading a vector to a file.
 	**/
     void load_from_file( const char* filename );
     
+    /**
+     * Comparision operators:
+    **/ 
+    bool operator == ( const Vector_distances_in_diagram& second )
+    {
+		if ( this->sorted_vector_of_distnaces.size() != second.sorted_vector_of_distnaces.size() )return false;
+		for ( size_t i = 0 ; i != this->sorted_vector_of_distnaces.size() ; ++i )
+		{
+			if ( !almost_equal(this->sorted_vector_of_distnaces[i] , second.sorted_vector_of_distnaces[i]) )return false;
+		}
+		return true;
+	}
+	
+	bool operator != ( const Vector_distances_in_diagram& second )
+    {
+		return !( *this == second );
+	}
+    
     //concretization of abstract methods:
     double project_to_R( int number_of_function );
     std::vector<double> vectorize( int number_of_function );
@@ -190,7 +209,8 @@ Vector_distances_in_diagram<F>::Vector_distances_in_diagram( std::vector< std::p
 {    
     std::vector< std::pair< double,double > > i( intervals_ );
     this->intervals = i;
-    this->compute_sorted_vector_of_distnaces_via_heap( where_to_cut );
+    //this->compute_sorted_vector_of_distnaces_via_heap( where_to_cut );
+    this->compute_sorted_vector_of_distnaces_via_vector_sorting(where_to_cut);
     this->set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
 }
 
@@ -218,15 +238,16 @@ Vector_distances_in_diagram<F>& Vector_distances_in_diagram<F>::operator =( cons
 
 template <typename F>
 Vector_distances_in_diagram<F>::Vector_distances_in_diagram( const char* filename , size_t where_to_cut  ):Abs_Vectorized_topological_data(where_to_cut)
-{   
+{   	
     //standard file with barcode
     std::vector< std::pair< double , double > > intervals = read_standard_file( filename );    
     //gudhi file with barcode
-    //std::vector< std::pair< double , double > > intervals = read_gudhi_file( filename , dimension );   
-    
+    //std::vector< std::pair< double , double > > intervals = read_gudhi_file( filename , dimension );    
+     
     this->intervals = intervals;
     this->compute_sorted_vector_of_distnaces_via_heap( where_to_cut );
-    set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();
+    //this->compute_sorted_vector_of_distnaces_via_vector_sorting( where_to_cut );
+    set_up_numbers_of_functions_for_vectorization_and_projections_to_reals();  
 }
 
 template < typename F>
@@ -241,10 +262,10 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_heap
         {
             cerr << this->intervals[i].first << " , " << this->intervals[i].second << endl;
         }
-    }    
+    }       
+	where_to_cut = std::min(where_to_cut , (size_t)(0.5 * this->intervals.size() * ( this->intervals.size() - 1 ) + this->intervals.size()));
 
-    std::vector< double > heap( where_to_cut );
-    std::fill( heap.begin() , heap.end() , std::numeric_limits<int>::max() );
+    std::vector< double > heap( where_to_cut , std::numeric_limits<int>::max() );
     std::make_heap (heap.begin(),heap.end());
     F f;
 
@@ -252,7 +273,7 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_heap
     for ( size_t i = 0 ; i < this->intervals.size() ; ++i )
     {
         for ( size_t j = i+1 ; j < this->intervals.size() ; ++j )
-        {
+        {			
             double value = std::min(
                                 f( this->intervals[i] , this->intervals[j] ),                               
                                 std::min(
@@ -261,13 +282,16 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_heap
                                         )    
                                 );
 
-                            //cerr << "Value : " << value << endl;
-                            //cerr << "heap.front() : " << heap.front() << endl;
-                            //getchar();
+							if ( dbg )
+							{
+								cerr << "Value : " << value << endl;
+								cerr << "heap.front() : " << heap.front() << endl;
+								getchar();
+							}
 
             if ( -value < heap.front() )
             {
-                //cerr << "Replacing : " << heap.front() << " with : " << -value << endl;getchar();
+                if ( dbg ){cerr << "Replacing : " << heap.front() << " with : " << -value << endl;getchar();}
                 //remove the first element from the heap
                 std::pop_heap (heap.begin(),heap.end());
                 //heap.pop_back();
@@ -309,6 +333,16 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_heap
             heap[i] *= -1;
         }
     }
+    
+    if ( dbg )
+    {
+		std::cerr << "This is the heap after all the operations :\n";
+		for ( size_t i = 0 ; i != heap.size() ; ++i )
+		{
+			std::cout << heap[i] << " ";
+		} 
+		std::cout << endl;
+	} 
 
     this->sorted_vector_of_distnaces = heap;
 }
@@ -318,11 +352,9 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_heap
 
 template < typename F>
 void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_vector_sorting( size_t where_to_cut )
-{
-    bool dbg = false;
- 
-	std::vector< double > distances;
-	distances.reserve( 0.5 * this->intervals.size() * ( this->intervals.size() - 1 ) + this->intervals.size() );
+{		
+	std::vector< double > distances;	
+	distances.reserve( (size_t)(0.5 * this->intervals.size() * ( this->intervals.size() - 1 ) + this->intervals.size()) );
     F f;
 
 	//for every pair of points in the diagram, compute the minimum of their distance, and distance of those points from diagonal	
@@ -343,8 +375,8 @@ void Vector_distances_in_diagram<F>::compute_sorted_vector_of_distnaces_via_vect
             
         }
     }
-    std::sort( distances.begin() , distances.end() , std::greater<double>() );
-    distances.resize( where_to_cut );
+    std::sort( distances.begin() , distances.end() , std::greater<double>() );    
+    if ( distances.size() > where_to_cut )distances.resize( where_to_cut );       
     
     this->sorted_vector_of_distnaces =  distances;
 }
@@ -410,28 +442,58 @@ template <typename F>
 double Vector_distances_in_diagram<F>::distance( const Abs_Topological_data_with_distances* second , double power )
 {
 	bool dbg = false;
-	if ( this->sorted_vector_of_distnaces.size() != ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces.size() )throw "Incompatible sizes of sorted_vector_of_distnaces in the method Vector_distances_in_diagram<F>::distance";
+	
+	if ( dbg )
+	{
+		std::cerr << "Entering double Vector_distances_in_diagram<F>::distance( const Abs_Topological_data_with_distances* second , double power ) procedure \n";
+		std::cerr << "Power : " << power << std::endl;
+		std::cerr << "This : " << *this << std::endl;
+		std::cerr << "second : " << *((Vector_distances_in_diagram<F>*)second) << std::endl;
+	}
+	
+	Vector_distances_in_diagram<F>* second_ = ((Vector_distances_in_diagram<F>*)second);
+	
 	
 	double result = 0;
-	for ( size_t i = 0 ; i != this->sorted_vector_of_distnaces.size() ; ++i )
+	for ( size_t i = 0 ; i != std::min(this->sorted_vector_of_distnaces.size(), second_->sorted_vector_of_distnaces.size())  ; ++i )
 	{
 		if ( power == 1 )
 		{
 			if ( dbg )
 			{
-				cerr << "|" << this->sorted_vector_of_distnaces[i] << " -  " << ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces[i] << " |  : " << fabs( this->sorted_vector_of_distnaces[i] - ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces[i] ) << endl;
+				cerr << "|" << this->sorted_vector_of_distnaces[i] << " -  " << second_->sorted_vector_of_distnaces[i] << " |  : " << fabs( this->sorted_vector_of_distnaces[i] - second_->sorted_vector_of_distnaces[i] ) << endl;
 			}
-			result += fabs( this->sorted_vector_of_distnaces[i] - ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces[i] );
+			result += fabs( this->sorted_vector_of_distnaces[i] - second_->sorted_vector_of_distnaces[i] );
 		}
 		else
 		{
-			result += std::pow( fabs( this->sorted_vector_of_distnaces[i] - ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces[i] ) , power );
+			result += std::pow( fabs( this->sorted_vector_of_distnaces[i] - second_->sorted_vector_of_distnaces[i] ) , power );
 			if ( dbg )
 			{
-				cerr << "| " << this->sorted_vector_of_distnaces[i] << " - " << ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces[i]  << " : " << fabs( this->sorted_vector_of_distnaces[i] - ((Vector_distances_in_diagram<F>*)second)->sorted_vector_of_distnaces[i] ) << endl;
+				cerr << "| " << this->sorted_vector_of_distnaces[i] << " - " << second_->sorted_vector_of_distnaces[i]  << " : " << fabs( this->sorted_vector_of_distnaces[i] - second_->sorted_vector_of_distnaces[i] ) << endl;
 			}
 		}
 	}
+	if ( this->sorted_vector_of_distnaces.size() != second_->sorted_vector_of_distnaces.size() )
+	{
+		if ( this->sorted_vector_of_distnaces.size() > second_->sorted_vector_of_distnaces.size() )
+		{
+			for ( size_t i = second_->sorted_vector_of_distnaces.size() ; i != this->sorted_vector_of_distnaces.size() ; ++i )			
+			{
+				result += fabs( this->sorted_vector_of_distnaces[i] );
+			}
+		}
+		else
+		{
+			//this->sorted_vector_of_distnaces.size() < second_->sorted_vector_of_distnaces.size()
+			for ( size_t i = this->sorted_vector_of_distnaces.size() ; i != second_->sorted_vector_of_distnaces.size() ; ++i )			
+			{
+				result += fabs( second_->sorted_vector_of_distnaces[i] );
+			}			
+		}
+	}
+	
+	
 	if ( power != 1 )
 	{
 		result = std::pow( result , (1.0/power) );
@@ -494,13 +556,9 @@ template < typename F>
 double Vector_distances_in_diagram<F>::compute_scalar_product( const Abs_Topological_data_with_scalar_product* second )
 {
 	Vector_distances_in_diagram<F>* second_vector = (Vector_distances_in_diagram<F>*)second;
-	if ( this->sorted_vector_of_distnaces.size() != second_vector->sorted_vector_of_distnaces.size() )
-	{
-		std::cerr << "Non compatible sizes of sorted vectors of distances. The program will now terminate \n";
-		throw "Non compatible sizes of sorted vectors of distances. The program will now terminate";
-	}	
+		
 	double result = 0;
-	for ( size_t i = 0 ; i != this->sorted_vector_of_distnaces.size() ; ++i )
+	for ( size_t i = 0 ; i != std::min(this->sorted_vector_of_distnaces.size(),second_vector->sorted_vector_of_distnaces.size()) ; ++i )
 	{
 		result += this->sorted_vector_of_distnaces[i] * second_vector->sorted_vector_of_distnaces[i];
 	}
diff --git a/src/Gudhi_stat/test/data/average_of_persistence_vectors b/src/Gudhi_stat/test/data/average_of_persistence_vectors
new file mode 100644
index 00000000..c7078160
--- /dev/null
+++ b/src/Gudhi_stat/test/data/average_of_persistence_vectors
@@ -0,0 +1 @@
+3.53553 3.53553 2.82843 0.707107 0.707107 0.707107 
\ No newline at end of file
diff --git a/src/Gudhi_stat/test/vector_representation_test.cpp b/src/Gudhi_stat/test/vector_representation_test.cpp
index bbbca93b..e87d1b23 100644
--- a/src/Gudhi_stat/test/vector_representation_test.cpp
+++ b/src/Gudhi_stat/test/vector_representation_test.cpp
@@ -21,8 +21,8 @@
  */
 
 
-
 #include <gudhi/concretizations/Vector_distances_in_diagram.h>
+#include <gudhi/common.h>
 #include <iostream>
 
 #define BOOST_TEST_DYN_LINK
@@ -33,139 +33,155 @@
 #include <cmath>
 #include <iomanip>
 
+
 using namespace Gudhi;
 using namespace Gudhi::Gudhi_stat;
 
 using namespace std;
 
 
-double epsilon = 0.000005;
 
 
 
 
+BOOST_AUTO_TEST_CASE(check_read_write_to_files) 
+{
+	std::vector< std::pair<double,double> > intervals;
+	intervals.push_back( std::make_pair(2,3) );
+	intervals.push_back( std::make_pair(4,7) );
+	intervals.push_back( std::make_pair(9,10) );
+	intervals.push_back( std::make_pair(3,11) );
+	Vector_distances_in_diagram< euclidean_distance<double> > p( intervals , -1 );
+	p.write_to_file( "test_vector_representation_write_read" );
+	
+	Vector_distances_in_diagram< euclidean_distance<double> > q;
+	q.load_from_file( "test_vector_representation_write_read" );
+	
+	BOOST_CHECK( p == q );
+	
+}
+
+
 BOOST_AUTO_TEST_CASE(check_sortev_vector_distances_template) 
 {	
 	Vector_distances_in_diagram< euclidean_distance<double> > p( "data/file_with_diagram" , 100 );
 	std::vector< double > sortev_vector_distances_template;
-	sortev_vector_distances_template.push_back(0.609968);
-	sortev_vector_distances_template.push_back(0.566317);
-	sortev_vector_distances_template.push_back(0.538858);
-	sortev_vector_distances_template.push_back(0.534927);
-	sortev_vector_distances_template.push_back(0.534927);
-	sortev_vector_distances_template.push_back(0.534927);
-	sortev_vector_distances_template.push_back(0.534927);
-	sortev_vector_distances_template.push_back(0.515741);
-	sortev_vector_distances_template.push_back(0.515741);
-	sortev_vector_distances_template.push_back(0.515741);
-	sortev_vector_distances_template.push_back(0.515741);
-	sortev_vector_distances_template.push_back(0.507828);
-	sortev_vector_distances_template.push_back(0.507828);
-	sortev_vector_distances_template.push_back(0.507828);
-	sortev_vector_distances_template.push_back(0.500911);
-	sortev_vector_distances_template.push_back(0.500911);
-	sortev_vector_distances_template.push_back(0.500911);
-	sortev_vector_distances_template.push_back(0.500911);
-	sortev_vector_distances_template.push_back(0.500911);
-	sortev_vector_distances_template.push_back(0.496986);
-	sortev_vector_distances_template.push_back(0.496986);
-	sortev_vector_distances_template.push_back(0.496986);
-	sortev_vector_distances_template.push_back(0.496986);
-	sortev_vector_distances_template.push_back(0.495306);
-	sortev_vector_distances_template.push_back(0.495306);
-	sortev_vector_distances_template.push_back(0.495306);
-	sortev_vector_distances_template.push_back(0.495306);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.439945);
-	sortev_vector_distances_template.push_back(0.424097);
-	sortev_vector_distances_template.push_back(0.424097);
-	sortev_vector_distances_template.push_back(0.424097);
-	sortev_vector_distances_template.push_back(0.424097);
-	sortev_vector_distances_template.push_back(0.424097);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.413891);
-	sortev_vector_distances_template.push_back(0.410613);
-	sortev_vector_distances_template.push_back(0.410613);
-	sortev_vector_distances_template.push_back(0.410613);
-	sortev_vector_distances_template.push_back(0.410613);
-	sortev_vector_distances_template.push_back(0.410613);
-	sortev_vector_distances_template.push_back(0.407853);
-	sortev_vector_distances_template.push_back(0.407853);
-	sortev_vector_distances_template.push_back(0.407853);
-	sortev_vector_distances_template.push_back(0.407853);
-	sortev_vector_distances_template.push_back(0.407853);
-	sortev_vector_distances_template.push_back(0.407853);
-	sortev_vector_distances_template.push_back(0.402306);
-	sortev_vector_distances_template.push_back(0.402306);
-	sortev_vector_distances_template.push_back(0.402306);
-	sortev_vector_distances_template.push_back(0.402306);
-	sortev_vector_distances_template.push_back(0.402306);
-	sortev_vector_distances_template.push_back(0.401937);
-	sortev_vector_distances_template.push_back(0.377605);
-	sortev_vector_distances_template.push_back(0.377605);
-	sortev_vector_distances_template.push_back(0.377605);
-	sortev_vector_distances_template.push_back(0.377605);
-	sortev_vector_distances_template.push_back(0.377605);
-	sortev_vector_distances_template.push_back(0.363859);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.357765);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.345124);
-	sortev_vector_distances_template.push_back(0.34469);
-	sortev_vector_distances_template.push_back(0.34469);  
+
+	sortev_vector_distances_template.push_back(	0.609968	);
+	sortev_vector_distances_template.push_back(	0.566317	);
+	sortev_vector_distances_template.push_back(	0.538858	);
+	sortev_vector_distances_template.push_back(	0.534927	);
+	sortev_vector_distances_template.push_back(	0.515741	);
+	sortev_vector_distances_template.push_back(	0.507828	);
+	sortev_vector_distances_template.push_back(	0.500911	);
+	sortev_vector_distances_template.push_back(	0.496986	);
+	sortev_vector_distances_template.push_back(	0.495306	);
+	sortev_vector_distances_template.push_back(	0.439945	);
+	sortev_vector_distances_template.push_back(	0.424097	);
+	sortev_vector_distances_template.push_back(	0.413891	);
+	sortev_vector_distances_template.push_back(	0.413891	);
+	sortev_vector_distances_template.push_back(	0.413891	);
+	sortev_vector_distances_template.push_back(	0.412621	);
+	sortev_vector_distances_template.push_back(	0.410613	);
+	sortev_vector_distances_template.push_back(	0.407853	);
+	sortev_vector_distances_template.push_back(	0.407853	);
+	sortev_vector_distances_template.push_back(	0.402306	);
+	sortev_vector_distances_template.push_back(	0.401937	);
+	sortev_vector_distances_template.push_back(	0.377605	);
+	sortev_vector_distances_template.push_back(	0.363859	);
+	sortev_vector_distances_template.push_back(	0.357765	);
+	sortev_vector_distances_template.push_back(	0.357765	);
+	sortev_vector_distances_template.push_back(	0.357765	);
+	sortev_vector_distances_template.push_back(	0.357765	);
+	sortev_vector_distances_template.push_back(	0.357765	);
+	sortev_vector_distances_template.push_back(	0.353401	);
+	sortev_vector_distances_template.push_back(	0.348004	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.345124	);
+	sortev_vector_distances_template.push_back(	0.34469	);
+	sortev_vector_distances_template.push_back(	0.339466	);
+	sortev_vector_distances_template.push_back(	0.33935	);
+	sortev_vector_distances_template.push_back(	0.32834	);
+	sortev_vector_distances_template.push_back(	0.327276	);
+	sortev_vector_distances_template.push_back(	0.318626	);
+	sortev_vector_distances_template.push_back(	0.318082	);
+	sortev_vector_distances_template.push_back(	0.30603	);
+	sortev_vector_distances_template.push_back(	0.30525	);
+	sortev_vector_distances_template.push_back(	0.297308	);
+	sortev_vector_distances_template.push_back(	0.296333	);
+	sortev_vector_distances_template.push_back(	0.296333	);
+	sortev_vector_distances_template.push_back(	0.296333	);
+	sortev_vector_distances_template.push_back(	0.296333	);
+	sortev_vector_distances_template.push_back(	0.293372	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.292666	);
+	sortev_vector_distances_template.push_back(	0.29029	);
+	sortev_vector_distances_template.push_back(	0.290218	);
+	sortev_vector_distances_template.push_back(	0.289782	);
+	sortev_vector_distances_template.push_back(	0.288128	);
+	sortev_vector_distances_template.push_back(	0.286416	);
+	sortev_vector_distances_template.push_back(	0.285969	);
+	sortev_vector_distances_template.push_back(	0.282046	);
+	sortev_vector_distances_template.push_back(	0.28154	);
+	sortev_vector_distances_template.push_back(	0.281085	);
+	sortev_vector_distances_template.push_back(	0.280227	);
+	sortev_vector_distances_template.push_back(	0.279273	);
+	sortev_vector_distances_template.push_back(	0.278936	);
+	sortev_vector_distances_template.push_back(	0.278706	);
+	sortev_vector_distances_template.push_back(	0.278507	);
+	sortev_vector_distances_template.push_back(	0.278097	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276293	);
+	sortev_vector_distances_template.push_back(	0.276169	);
+	sortev_vector_distances_template.push_back(	0.270563	);
+	sortev_vector_distances_template.push_back(	0.264009	);
+
+
+ 
 	size_t proj_no = p.number_of_vectorize_functions();
 	std::vector< double > aa = p.vectorize(proj_no);
 	
 	for ( size_t i = 0 ; i != aa.size() ; ++i )	
-	{
-		BOOST_CHECK( fabs ( sortev_vector_distances_template[i] - aa[i] ) <= epsilon );
+	{		
+		BOOST_CHECK( almost_equal( sortev_vector_distances_template[i] , aa[i] ) );
 	}
 }	
 
 
 
-	
-	
-	
-	
-	
-	
-	
+
 BOOST_AUTO_TEST_CASE(check_projections_to_R) 
 {
 	Vector_distances_in_diagram< euclidean_distance<double> > p( "data/file_with_diagram" , 100 );
@@ -175,108 +191,106 @@ BOOST_AUTO_TEST_CASE(check_projections_to_R)
 	proj.push_back(	1.176284775	);
 	proj.push_back(	1.715142954	);
 	proj.push_back(	2.25006986	);
-	proj.push_back(	2.784996767	);
-	proj.push_back(	3.319923673	);
-	proj.push_back(	3.854850579	);
-	proj.push_back(	4.370591225	);
-	proj.push_back(	4.886331872	);
-	proj.push_back(	5.402072518	);
-	proj.push_back(	5.917813164	);
-	proj.push_back(	6.425641089	);
-	proj.push_back(	6.933469014	);
-	proj.push_back(	7.441296939	);
-	proj.push_back(	7.942207817	);
-	proj.push_back(	8.443118694	);
-	proj.push_back(	8.944029572	);
-	proj.push_back(	9.44494045	);
-	proj.push_back(	9.945851328	);
-	proj.push_back(	10.44283706	);
-	proj.push_back(	10.93982279	);
-	proj.push_back(	11.43680853	);
-	proj.push_back(	11.93379426	);
-	proj.push_back(	12.42910009	);
-	proj.push_back(	12.92440592	);
-	proj.push_back(	13.41971176	);
-	proj.push_back(	13.91501759	);
-	proj.push_back(	14.35496286	);
-	proj.push_back(	14.79490814	);
-	proj.push_back(	15.23485341	);
-	proj.push_back(	15.67479869	);
-	proj.push_back(	16.11474396	);
-	proj.push_back(	16.55468923	);
-	proj.push_back(	16.99463451	);
-	proj.push_back(	17.41873131	);
-	proj.push_back(	17.84282811	);
-	proj.push_back(	18.26692491	);
-	proj.push_back(	18.69102171	);
-	proj.push_back(	19.11511851	);
-	proj.push_back(	19.52900989	);
-	proj.push_back(	19.94290127	);
-	proj.push_back(	20.35679265	);
-	proj.push_back(	20.77068403	);
-	proj.push_back(	21.18457541	);
-	proj.push_back(	21.59846679	);
-	proj.push_back(	22.01235817	);
-	proj.push_back(	22.42297162	);
-	proj.push_back(	22.83358506	);
-	proj.push_back(	23.2441985	);
-	proj.push_back(	23.65481194	);
-	proj.push_back(	24.06542538	);
-	proj.push_back(	24.47327795	);
-	proj.push_back(	24.88113051	);
-	proj.push_back(	25.28898308	);
-	proj.push_back(	25.69683564	);
-	proj.push_back(	26.1046882	);
-	proj.push_back(	26.51254077	);
-	proj.push_back(	26.91484703	);
-	proj.push_back(	27.3171533	);
-	proj.push_back(	27.71945957	);
-	proj.push_back(	28.12176583	);
-	proj.push_back(	28.5240721	);
-	proj.push_back(	28.92600955	);
-	proj.push_back(	29.30361425	);
-	proj.push_back(	29.68121895	);
-	proj.push_back(	30.05882365	);
-	proj.push_back(	30.43642834	);
-	proj.push_back(	30.81403304	);
-	proj.push_back(	31.17789162	);
-	proj.push_back(	31.53565711	);
-	proj.push_back(	31.89342259	);
-	proj.push_back(	32.25118808	);
-	proj.push_back(	32.60895356	);
-	proj.push_back(	32.96671905	);
-	proj.push_back(	33.32448454	);
-	proj.push_back(	33.68225002	);
-	proj.push_back(	34.04001551	);
-	proj.push_back(	34.397781	);
-	proj.push_back(	34.75554648	);
-	proj.push_back(	35.11331197	);
-	proj.push_back(	35.47107745	);
-	proj.push_back(	35.82884294	);
-	proj.push_back(	36.18660843	);
-	proj.push_back(	36.53173245	);
-	proj.push_back(	36.87685646	);
-	proj.push_back(	37.22198048	);
-	proj.push_back(	37.5671045	);
-	proj.push_back(	37.91222852	);
-	proj.push_back(	38.25735254	);
-	proj.push_back(	38.60247656	);
-	proj.push_back(	38.94760057	);
-	proj.push_back(	39.29272459	);
-	proj.push_back(	39.63784861	);
-	proj.push_back(	39.98297263	);
-	proj.push_back(	40.32809665	);
-	proj.push_back(	40.67322067	);
-	proj.push_back(	41.01834468	);
-	proj.push_back(	41.3634687	);
-	proj.push_back(	41.70815824	);
-
-
+	proj.push_back(	2.765810506	);
+	proj.push_back(	3.273638431	);
+	proj.push_back(	3.774549309	);
+	proj.push_back(	4.271535042	);
+	proj.push_back(	4.766840875	);
+	proj.push_back(	5.206786149	);
+	proj.push_back(	5.63088295	);
+	proj.push_back(	6.04477433	);
+	proj.push_back(	6.45866571	);
+	proj.push_back(	6.87255709	);
+	proj.push_back(	7.285177939	);
+	proj.push_back(	7.695791381	);
+	proj.push_back(	8.103643945	);
+	proj.push_back(	8.511496508	);
+	proj.push_back(	8.913802775	);
+	proj.push_back(	9.315740229	);
+	proj.push_back(	9.693344927	);
+	proj.push_back(	10.0572035	);
+	proj.push_back(	10.41496899	);
+	proj.push_back(	10.77273448	);
+	proj.push_back(	11.13049996	);
+	proj.push_back(	11.48826545	);
+	proj.push_back(	11.84603094	);
+	proj.push_back(	12.19943233	);
+	proj.push_back(	12.5474364	);
+	proj.push_back(	12.89256042	);
+	proj.push_back(	13.23768444	);
+	proj.push_back(	13.58280846	);
+	proj.push_back(	13.92793248	);
+	proj.push_back(	14.2730565	);
+	proj.push_back(	14.61818051	);
+	proj.push_back(	14.96330453	);
+	proj.push_back(	15.30842855	);
+	proj.push_back(	15.65355257	);
+	proj.push_back(	15.99867659	);
+	proj.push_back(	16.34380061	);
+	proj.push_back(	16.68892462	);
+	proj.push_back(	17.03404864	);
+	proj.push_back(	17.37917266	);
+	proj.push_back(	17.7238622	);
+	proj.push_back(	18.06332781	);
+	proj.push_back(	18.40267754	);
+	proj.push_back(	18.73101759	);
+	proj.push_back(	19.05829313	);
+	proj.push_back(	19.3769189	);
+	proj.push_back(	19.69500045	);
+	proj.push_back(	20.0010306	);
+	proj.push_back(	20.30628026	);
+	proj.push_back(	20.60358868	);
+	proj.push_back(	20.89992192	);
+	proj.push_back(	21.19625516	);
+	proj.push_back(	21.4925884	);
+	proj.push_back(	21.78892164	);
+	proj.push_back(	22.08229394	);
+	proj.push_back(	22.37495987	);
+	proj.push_back(	22.66762581	);
+	proj.push_back(	22.96029174	);
+	proj.push_back(	23.25295768	);
+	proj.push_back(	23.54562361	);
+	proj.push_back(	23.83828955	);
+	proj.push_back(	24.13095549	);
+	proj.push_back(	24.42362142	);
+	proj.push_back(	24.71628736	);
+	proj.push_back(	25.00895329	);
+	proj.push_back(	25.30161923	);
+	proj.push_back(	25.59428516	);
+	proj.push_back(	25.8869511	);
+	proj.push_back(	26.17961703	);
+	proj.push_back(	26.47228297	);
+	proj.push_back(	26.76257262	);
+	proj.push_back(	27.05279049	);
+	proj.push_back(	27.34257265	);
+	proj.push_back(	27.63070097	);
+	proj.push_back(	27.91711687	);
+	proj.push_back(	28.20308566	);
+	proj.push_back(	28.48513176	);
+	proj.push_back(	28.76667161	);
+	proj.push_back(	29.04775635	);
+	proj.push_back(	29.32798359	);
+	proj.push_back(	29.60725702	);
+	proj.push_back(	29.88619335	);
+	proj.push_back(	30.16489915	);
+	proj.push_back(	30.44340655	);
+	proj.push_back(	30.72150329	);
+	proj.push_back(	30.99779604	);
+	proj.push_back(	31.27408878	);
+	proj.push_back(	31.55038153	);
+	proj.push_back(	31.82667427	);
+	proj.push_back(	32.10296702	);
+	proj.push_back(	32.37925976	);
+	proj.push_back(	32.6555525	);
+	proj.push_back(	32.93184525	);
+	proj.push_back(	33.20813799	);
+	proj.push_back(	33.48430662	);
+	proj.push_back(	33.7548692	);
 
 	for ( size_t proj_no = 0 ; proj_no != p.number_of_projections_to_R() ; ++proj_no )	 
 	{		
 		//cout << std::setprecision(10)  << 	p.project_to_R(proj_no) << endl;
-		BOOST_CHECK( fabs ( p.project_to_R(proj_no) - proj[proj_no] ) <= epsilon );
+		BOOST_CHECK( almost_equal( p.project_to_R(proj_no) , proj[proj_no] ) );
 	}
 }
 
@@ -299,8 +313,9 @@ BOOST_AUTO_TEST_CASE(check_distance_computations)
 	intervals[7] = std::pair<double,double>( 15,16 );
 	intervals[8] = std::pair<double,double>( 17,18 );
 	intervals[9] = std::pair<double,double>( 19,20 );		
-	Vector_distances_in_diagram< euclidean_distance<double> > p_bis( intervals , 10 );		
-	BOOST_CHECK( fabs ( p_prime.distance( (Abs_Topological_data_with_distances*)(&p_bis) , 1 ) - 1.669 ) <= epsilon );		
+	Vector_distances_in_diagram< euclidean_distance<double> > p_bis( intervals , 10 );	
+	//cerr << "p_prime.distance( (Abs_Topological_data_with_distances*)(&p_bis) , 1 ) : " << p_prime.distance( (Abs_Topological_data_with_distances*)(&p_bis) , 1 ) << endl;	
+	BOOST_CHECK( almost_equal ( p_prime.distance( (Abs_Topological_data_with_distances*)(&p_bis) , 1 ) , 1.86428 ) );		
 }
 
 
@@ -311,16 +326,16 @@ BOOST_AUTO_TEST_CASE(check_compute_average)
 	//compute average 
 	std::vector< std::pair<double,double> > i1(3);
 	i1[0] = std::pair<double,double>( 1,2 );
-	i1[1] = std::pair<double,double>( 3,4 );
-	i1[2] = std::pair<double,double>( 5,6 );
+	i1[1] = std::pair<double,double>( 3,8 );
+	i1[2] = std::pair<double,double>( 1,6 );
 	
 	std::vector< std::pair<double,double> > i2(3);
-	i2[0] = std::pair<double,double>( 2,3 );
-	i2[1] = std::pair<double,double>( 4,5);
-	i2[2] = std::pair<double,double>( 6,7 );
+	i2[0] = std::pair<double,double>( 2,9 );
+	i2[1] = std::pair<double,double>( 2,15);
+	i2[2] = std::pair<double,double>( 6,17 );
 	
-	Vector_distances_in_diagram< euclidean_distance<double> > A( i1 , 3 );
-	Vector_distances_in_diagram< euclidean_distance<double> > B( i1 , 3 );
+	Vector_distances_in_diagram< euclidean_distance<double> > A( i1 , -1 );
+	Vector_distances_in_diagram< euclidean_distance<double> > B( i1 , -1 );
 	
 	std::vector< Abs_Topological_data_with_averages* > to_average;
 	to_average.push_back( (Abs_Topological_data_with_averages*)(&A) );
@@ -329,8 +344,10 @@ BOOST_AUTO_TEST_CASE(check_compute_average)
 	Vector_distances_in_diagram< euclidean_distance<double> > average; 
 	average.compute_average( to_average );
 	
-	
-	std::vector< double > vect = average.output_for_visualization();
+	Vector_distances_in_diagram< euclidean_distance<double> > template_average; 
+	template_average.load_from_file( "data/average_of_persistence_vectors" );
+		
+	BOOST_CHECK( template_average == average );
 	
 }