adding a few corrections. Moreover I right now persisgence_heat_map is a template class. A template parameter is a class that gives the scaling of kernel with repsect to diagonal.

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

Former-commit-id: 37e4c01a95f5682f4acb3bba3145014671f95188

4 files changed, 129 insertions, 88 deletions

diff --git a/src/Gudhi_stat/include/gudhi/concretizations/PSSK.h b/src/Gudhi_stat/include/gudhi/concretizations/PSSK.h
index 8168b7a7..e61711f4 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/PSSK.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/PSSK.h
@@ -41,7 +41,7 @@ namespace Gudhi_stat
 * at (q,p)
 **/
 
-class PSSK : public Persistence_heat_maps
+class PSSK : public Persistence_heat_maps<constant_scaling_function>
 {
 public:
 	PSSK():Persistence_heat_maps(){}
@@ -103,7 +103,7 @@ void PSSK::construct( const std::vector< std::pair<double,double> >& intervals_
 
     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
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 d2c1bfc9..f846ac08 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
@@ -122,26 +122,42 @@ std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius ,
  * (4) Scale them with 
 **/ 
 
-double constant_function( const std::pair< double , double >& point_in_diagram )
-{	
-	return 1;
-}
+class constant_scaling_function
+{
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{	
+		return 1;
+	}
+};
 
-double distance_from_diagonal( const std::pair< double , double >& point_in_diagram )
+class distance_from_diagonal_scalling
 {
-	//(point_in_diagram.first+point_in_diagram.second)/2.0
-	return sqrt( pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2)  );
-}
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		//(point_in_diagram.first+point_in_diagram.second)/2.0
+		return sqrt( pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2)  );
+	}
+};
 
-double squared_distance_from_diagonal( const std::pair< double , double >& point_in_diagram )
+class squared_distance_from_diagonal_scaling
 {
-	return pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2);
-}
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		return pow((point_in_diagram.first-(point_in_diagram.first+point_in_diagram.second)/2.0),2) + pow((point_in_diagram.second-(point_in_diagram.first+point_in_diagram.second)/2.0),2);
+	}
+};
 
-double arc_tan_of_persistence_of_point( const std::pair< double , double >& point_in_diagram )
+class arc_tan_of_persistence_of_point
 {
-	return atan( point_in_diagram.second - point_in_diagram.first );
-}
+public:
+	double operator()( const std::pair< double , double >& point_in_diagram )
+	{
+		return atan( point_in_diagram.second - point_in_diagram.first );
+	}
+};
 
 class weight_by_setting_maximal_interval_to_have_length_one
 {
@@ -159,6 +175,7 @@ private:
 /**
  * This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product
 **/ 
+template <typename Scalling_of_kernels = constant_scaling_function>
 class Persistence_heat_maps
 {
 public:
@@ -167,7 +184,8 @@ public:
     **/ 
 	Persistence_heat_maps()
 	{
-		this->scalling_function_with_respect_to_distance_from_diagonal = constant_function;
+		Scalling_of_kernels f;
+		this->f = f;
 		this->erase_below_diagonal = false;		
 		this->min_ = this->max_ = 0;
 		this->set_up_parameters_for_basic_classes();
@@ -177,13 +195,12 @@ public:
 	 * Construction that takes at the input the following parameters:
 	 * (1) A vector of pairs of doubles (representing persistence intervals). All other parameters are optional. They are:
 	 * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), 
-	 * (3) a scaling function (constant by default). This function allows to put a different weights to points depending on their position, in particular, depending on their distance from diagonal, 
-	 * (4) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). 
-	 * (5) a number of pixels in each direction (set to 1000 by default). 
-	 * (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,
+	 * (3) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). 
+	 * (4) a number of pixels in each direction (set to 1000 by default). 
+	 * (5) a min x and y value of points that are to be taken into account. By default it is set to -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 std::vector< std::pair< double,double > >  & 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 std::vector< std::pair< double,double > >  & interval , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1) , bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  );
     
     /**
 	 * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by default), a scaling function (constant by default), a boolean value which determines if the area of image below diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each direction (set to 1000 by default). and min and max values of images (both set to -1 by defaulet. If this is the case, the program will pick the right values based on the data).
@@ -191,14 +208,13 @@ public:
 	/**
 	 * Construction that takes at the input the following parameters:
 	 * (1) A a name of a file with persistence intervals. The file shold be readable by the function read_standard_file. All other parameters are optional. They are:
-	 * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), 
-	 * (3) a scaling function (constant by default). This function allows to put a different weights to points depending on their position, in particular, depending on their distance from diagonal, 
-	 * (4) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). 
-	 * (5) a number of pixels in each direction (set to 1000 by default). 
-	 * (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,
+	 * (2) a Gausian filter generated by create_Gaussian_filter filter (the default value of this vaiable is a Gaussian filter of a radius 5), 	
+	 * (3) a boolean value which determines if the area of image below diagonal should, or should not be erased (it will be erased by default). 
+	 * (4) a number of pixels in each direction (set to 1000 by default). 
+	 * (5) a min x and y value of points that are to be taken into account. By default it is set to -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* 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  );
+    Persistence_heat_maps( const char* filename , std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1), bool erase_below_diagonal = false , size_t number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  );
 
 	
 	/**
@@ -368,7 +384,7 @@ protected:
 	size_t number_of_functions_for_projections_to_reals;
     void construct( const std::vector< std::pair<double,double> >& intervals_  , 
 					std::vector< std::vector<double> > filter = create_Gaussian_filter(5,1),
-                    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 );
                     
 	void set_up_parameters_for_basic_classes()
@@ -378,7 +394,8 @@ protected:
 	}
     
     //data    
-    double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram );    
+    //double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram );    
+    Scalling_of_kernels f;
     bool erase_below_diagonal;
     double min_;
     double max_;
@@ -387,13 +404,15 @@ protected:
 
 
 //if min_ == max_, then the program is requested to set up the values itself based on persistence intervals
-void Persistence_heat_maps::construct( const std::vector< std::pair<double,double> >& intervals_  ,  
-									   std::vector< std::vector<double> > filter,
-									   double (*scalling_function_with_respect_to_distance_from_diagonal)( const std::pair< double , double >& point_in_diagram ), 
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::construct( const std::vector< std::pair<double,double> >& intervals_  ,  
+									   std::vector< std::vector<double> > filter,									   
 									   bool erase_below_diagonal , size_t number_of_pixels , double min_ , double max_ )
 {
     bool dbg = false;       
     if ( dbg )std::cerr << "Entering construct procedure \n";
+    Scalling_of_kernels f;
+    this->f = f;
 
     if ( min_ == max_ )
     {
@@ -423,9 +442,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;
-
+    this->max_ = max_;        
 
     //initialization of the structure heat_map
     std::vector< std::vector<double> > heat_map_;
@@ -463,7 +480,7 @@ void Persistence_heat_maps::construct( const std::vector< std::pair<double,doubl
 			std::cerr << "y_grid : " << y_grid << endl;
 		}
 				
-		double scaling_value = this->scalling_function_with_respect_to_distance_from_diagonal(intervals_[pt_nr]);
+		double scaling_value = this->f(intervals_[pt_nr]);
 
 		
 		for ( size_t i = 0 ; i != filter.size() ; ++i )
@@ -502,30 +519,29 @@ void Persistence_heat_maps::construct( const std::vector< std::pair<double,doubl
 	}
 }//construct
 
-
-Persistence_heat_maps::Persistence_heat_maps( const std::vector< std::pair< double,double > >  & interval , 
+template <typename Scalling_of_kernels>
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( const std::vector< std::pair< double,double > >  & interval , 
 											  std::vector< std::vector<double> > filter,  
-											  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_ )
 {
-    this->construct( interval ,  filter ,  scalling_function_with_respect_to_distance_from_diagonal, erase_below_diagonal , number_of_pixels , min_ , max_ );
+    this->construct( interval ,  filter , erase_below_diagonal , number_of_pixels , min_ , max_ );
     this->set_up_parameters_for_basic_classes();
 }
 
 
-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 ), 
+template <typename Scalling_of_kernels>
+Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( const char* filename , 
+											  std::vector< std::vector<double> > filter, 											  
 											  bool erase_below_diagonal , size_t number_of_pixels  , double min_ , double max_ )
 {    
     std::vector< std::pair< double , double > > intervals_ = read_standard_file( filename );           
-    this->construct( intervals_ ,  filter ,  scalling_function_with_respect_to_distance_from_diagonal, erase_below_diagonal , number_of_pixels , min_ , max_ );
+    this->construct( intervals_ ,  filter, erase_below_diagonal , number_of_pixels , min_ , max_ );
     this->set_up_parameters_for_basic_classes();
 }
 
 
-
-std::vector< std::vector<double> > Persistence_heat_maps::check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps )
+template <typename Scalling_of_kernels>
+std::vector< std::vector<double> > Persistence_heat_maps<Scalling_of_kernels>::check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps )
 {
 	//checking if all the heat maps are of the same size:
 	for ( size_t i = 0 ; i != maps.size() ; ++i )
@@ -550,8 +566,8 @@ std::vector< std::vector<double> > Persistence_heat_maps::check_and_initialize_m
 	return heat_maps;
 }
 
-
-void Persistence_heat_maps::compute_median( const std::vector<Persistence_heat_maps*>& maps )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_median( const std::vector<Persistence_heat_maps*>& maps )
 {
 	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
 	
@@ -574,8 +590,8 @@ 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 )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_mean( const std::vector<Persistence_heat_maps*>& maps )
 {
 	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
 		    for ( size_t i = 0 ; i != heat_maps.size() ; ++i )
@@ -597,8 +613,8 @@ void Persistence_heat_maps::compute_mean( const std::vector<Persistence_heat_map
 
 
 
-
-void Persistence_heat_maps::compute_percentage_of_active( const std::vector<Persistence_heat_maps*>& maps , size_t cutoff  )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_percentage_of_active( const std::vector<Persistence_heat_maps*>& maps , size_t cutoff  )
 {
 	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
 
@@ -627,8 +643,8 @@ void Persistence_heat_maps::compute_percentage_of_active( const std::vector<Pers
 }
 
 
-
-void Persistence_heat_maps::plot( const char* filename )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::plot( const char* filename )
 {
 	ofstream out;
 	std::stringstream ss;
@@ -649,8 +665,8 @@ void Persistence_heat_maps::plot( const char* filename )
 }
 
 
-
-void Persistence_heat_maps::print_to_file( const char* filename )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::print_to_file( const char* filename )
 {
 
 	ofstream out;
@@ -669,8 +685,8 @@ void Persistence_heat_maps::print_to_file( const char* filename )
 	out.close();
 }
 
-
-void Persistence_heat_maps::load_from_file( const char* filename )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::load_from_file( const char* filename )
 {
 	bool dbg = false;	
 	
@@ -728,8 +744,9 @@ 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 )
+//Concretizations of virtual methods:   
+template <typename Scalling_of_kernels> 
+std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize( int number_of_function )
 {
 	//convert this->heat_map into one large vector:
 	size_t size_of_result = 0;
@@ -752,7 +769,8 @@ std::vector<double> Persistence_heat_maps::vectorize( int number_of_function )
 	return result;
 }
 
-double Persistence_heat_maps::distance( const Persistence_heat_maps& second , double power )
+template <typename Scalling_of_kernels>
+double Persistence_heat_maps<Scalling_of_kernels>::distance( const Persistence_heat_maps& second , double power )
 {			
 	//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) )
@@ -774,7 +792,8 @@ double Persistence_heat_maps::distance( const Persistence_heat_maps& second , do
 	return distance;
 }
 
-double Persistence_heat_maps::project_to_R( int number_of_function )
+template <typename Scalling_of_kernels>
+double Persistence_heat_maps<Scalling_of_kernels>::project_to_R( int number_of_function )
 {
 	double result = 0;
 	for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
@@ -787,12 +806,14 @@ double Persistence_heat_maps::project_to_R( int number_of_function )
 	return result;
 }
 
-void Persistence_heat_maps::compute_average( const std::vector< Persistence_heat_maps* >& to_average )
+template <typename Scalling_of_kernels>
+void Persistence_heat_maps<Scalling_of_kernels>::compute_average( const std::vector< Persistence_heat_maps* >& to_average )
 {				
 	this->compute_mean( to_average );
 }
 
-double Persistence_heat_maps::compute_scalar_product( const Persistence_heat_maps& second )
+template <typename Scalling_of_kernels>
+double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product( const 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) )
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h b/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h
index 0a740996..953547c7 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/read_persitence_from_file.h
@@ -105,7 +105,7 @@ std::vector< std::pair< double,double > > read_file_names_that_may_contain_inf_s
 /**
  * This procedure reads names of files which are stored in a file. 
 **/ 
-std::vector< std::string > readFileNames( char* filenameWithFilenames )
+std::vector< std::string > readFileNames( const char* filenameWithFilenames )
 {
     bool dbg = false;
     
diff --git a/src/Gudhi_stat/include/gudhi/topological_process.h b/src/Gudhi_stat/include/gudhi/topological_process.h
index d8d8e696..729ef2f8 100644
--- a/src/Gudhi_stat/include/gudhi/topological_process.h
+++ b/src/Gudhi_stat/include/gudhi/topological_process.h
@@ -20,12 +20,8 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef PERMUTATION_TEST_H
-#define PERMUTATION_TEST_H
-//abstract classes:
-#include <gudhi/abstract_classes/Abs_Topological_data.h>
-#include <gudhi/abstract_classes/Abs_Topological_data_with_averages.h>
-#include <gudhi/abstract_classes/Abs_Topological_data_with_distances.h>
+#ifndef TOPOLOGICAL_PROCESS_H
+#define TOPOLOGICAL_PROCESS_H
 
 
 //concretizations
@@ -41,14 +37,48 @@ namespace Gudhi
 {
 namespace Gudhi_stat 
 {
+	
+
+template <typename Representation>
+std::vector< Representation* > construct_representation_from_file( const char* filename )	
+{
+	bool dbg = false;
+	std::vector< std::string > files = readFileNames( filename );	
+	
+	std::cout << "Here are the filenames in the file : " << filename << std::endl;
+	for ( size_t i = 0 ; i != files.size() ; ++i )
+	{
+		std::cout << files[i] << std::endl;
+	}
+	
+	std::vector< Representation* > result( files.size() );
+	for ( size_t i = 0 ; i != files.size() ; ++i )
+	{
+		std::vector< std::pair< double , double > > diag = read_standard_file( files[i].c_str() );	
+		
+		if ( dbg )
+		{
+			std::cerr << "Here is a diagram from a file : " << files[i].c_str() << std::endl;
+			for ( size_t aa = 0  ; aa != diag.size() ; ++aa )
+			{
+				std::cout << diag[aa].first << " " << diag[aa].second << std::endl;
+			}
+			getchar();
+		}
+						
+		Representation* l = new Representation( diag );			
+		result[i] = l;		
+	}
+	return result;
+}
+	
 
 template <typename Representation>
 class Topological_process
 {
 public:
 	Topological_process();
-	Topological_process( const std::vetor< Representation* >& data_ ):data(data_){}
-	Topological_process( const char* filename );//we assume that the filename contains a list of files with the diagrams. 
+	Topological_process( const std::vector< Representation* >& data_ ):data(data_){}
 	double distance( const Topological_process& second )
 	{
 		if ( this->data.size() != second.data.size() )
@@ -58,12 +88,12 @@ public:
 		double result = 0;
 		for ( size_t i = 0 ; i != this->data.size() ; ++i )
 		{
-			result += this->data[i].distance( secod.data[i] );
+			result += this->data[i]->distance( *second.data[i] );
 		}
 		return result;
 	}	
 	
-	compute_average( const std::vector< Representation* >& to_average )
+	void compute_average( const std::vector< Representation* >& to_average )
 	{
 		//since we will substitute whatever data we have in this object with an average, we clear the data in this object first:
 		this->data.clear();
@@ -80,20 +110,10 @@ public:
 	//scalar products?
 	//confidence bounds?
 private:
-	std::vetor< Representation* > data;
+	std::vector< Representation* > data;
 };
 
-template <typename Representation>
-Topological_process<Representation>::Topological_process( const char* filename )
-{
-	std::vector< std::string > filenames = readFileNames( filename );
-	this->data.reserve( filenames.size() );
-	for ( size_t file = 0 ; file != filenames.size() ; ++file )
-	{
-		std::vector< std::pair<double,double> > intervals_in_this_file = read_file_names_that_may_contain_inf_string( filenames[i] );
-		this-data[i] = new Representation( intervals_in_this_file );
-	}	
-}
+
 
 }//Gudhi_stat
 }//Gudhi