a few modification to the code to compute persistence landscapes on a grid

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

Former-commit-id: 6adb9cecf2465330c121e44236ad1d7c4bce2d07

4 files changed, 146 insertions, 167 deletions

diff --git a/src/Gudhi_stat/example/CMakeLists.txt b/src/Gudhi_stat/example/CMakeLists.txt
index aadfc83b..0229ba7c 100644
--- a/src/Gudhi_stat/example/CMakeLists.txt
+++ b/src/Gudhi_stat/example/CMakeLists.txt
@@ -1,6 +1,9 @@
 cmake_minimum_required(VERSION 2.6)
 project(GUDHI_STAT)
 
+add_executable ( persistence_landscape_on_grid persistence_landscape_on_grid.cpp )
+target_link_libraries(persistence_landscape_on_grid ${Boost_SYSTEM_LIBRARY})
+
 add_executable ( persistence_landscape persistence_landscape.cpp )
 target_link_libraries(persistence_landscape ${Boost_SYSTEM_LIBRARY})
 
@@ -13,7 +16,6 @@ target_link_libraries(create_landscapes ${Boost_SYSTEM_LIBRARY})
 add_executable ( plot_landscapes plot_landscapes.cpp )
 target_link_libraries(plot_landscapes ${Boost_SYSTEM_LIBRARY})
 
-
 add_executable ( persistence_intervals persistence_intervals.cpp )
 target_link_libraries(persistence_intervals ${Boost_SYSTEM_LIBRARY})
 
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 a0a74eba..c3c665b2 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_heat_maps.h
@@ -20,35 +20,6 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
  
- 
- /*
-I would not say that homology is rovust when you speak about invariance with respect to continuous deformation. 
-* 
-Definition of simplicial compelx is redundant. There is no need to include vertices separatell.
-* 
-term 'codimension' used without being defined.
-* 
-"because boundary of boundary is always empty' -- in an intro survey one should not assume it.
-* 
-*what you call 'simple homotopies' in the literature is called 'free face collapses'. They are based on a simple homotopy theory, but this is not how they are called.
-* 
-* "hovewer when 1 \neq -1 ... -- not clear for people who do not know what you are talking about.
-* 
-* It is a urban legend that Cech complexes are difficult. They are not, one use exactly the same algorithm as in Rips complex construction and use minibal to check the intersection. 
-* 
-* correlations or measures of dissimilarity are typically not distances.
-* 
-* topological inference means soemthing else than you say in section 4 ?????????????????????????????????
-* 
-* Definition of persistence is wrong. Authors in the next paragraph are trying to provide a 'refined and ultimately more useful' definition which takes into account that classes can merge.
-* This is a very clumzy and extremally confusing and in my opinion useless. What is the aim of providing wrong definition and then trying to fix it in informal text? Now it is not clear what do 
-* they mean by birth and death later in the paper.
-* In the previous iteration the feedback on this was given by both reviewers. Still, the authors failed to copy-pasted the correct definition 
-* from vast literature. The fact that despite critical comments in the previous review none of many authors of the paper have found a problem in definition clearly indicate that authors are not 
-* experts in persistent homology. In my opinion such a type of survey should be written by scientists who are more experiences in a field. 
-* 
-* 
-  */
 
 
 //standard include
@@ -57,20 +28,73 @@ term 'codimension' used without being defined.
 #include <iostream>
 
 //gudhi include
-#include "persistence_intervals.h"
+#include "read_persitence_from_file.h"
 #include "Distances_of_points_in_diagram.h"
-#include "read_files_names.h"
 
 
 using namespace std;
 
-template <typename T>
+std::vector< std::vector<double> > create_Gaussian_filter( size_t pixel_radius , double sigma )
+{
+	//we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a size 2*pixel_radius times 2*pixel_radius.
+	
+    double r = 0;
+    double sigma_sqr = sigma * sigma;
+ 
+    // sum is for normalization
+    double sum = 0;
+    
+    //initialization of a kernel:
+    std::vector< std::vector<double> > kernel( 2*pixel_radius  );
+    for ( size_t i = 0 ; i != 2*pixel_radius ; ++i )
+    {
+		std::vector<double> v( 2*pixel_radius , 0 );
+		kernel[i] = v;
+	}
+ 
+    for (int x = -pixel_radius; x <= pixel_radius; x++)
+    {
+        for(int y = -pixel_radius; y <= pixel_radius; y++)
+        {
+			double real_x = 2*sigma*x/pixel_radius;
+			double real_y = 2*sigma*y/pixel_radius;						
+            r = sqrt(real_x*real_x + real_y*real_y);
+            kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r*r)/sigma_sqr))/(3.141592 * sigma_sqr);
+            sum += gKernel[x + pixel_radius][y + pixel_radius];
+        }
+    }
+ 
+    // normalize the kernel
+    for(int i = 0; i != kernel.size() ; ++i)
+    {
+        for(int j = 0; j != kernel[i].size() ; ++j)
+        {
+            gKernel[i][j] /= sum;
+		}
+            
+    }
+    return kernel;
+}
+
+double constant_function( std::pair< double , double >& point_in_diagram )
+{
+	return 1;
+}
+
+
+
 class Persistence_heat_maps
 {
 public:
-	Persistence_heat_maps(){};
-    Persistence_heat_maps( Persistence_intervals<T>& interval , int number_of_pixels = 1000 , int diameter_of_gaussian = 4 , T min_ = -1 , T max_ = -1  );
-    Persistence_heat_maps( char* name_of_file_with_names_of_files_with_intervals , int number_of_pixels = 1000 , int diameter_of_gaussian = 4 , T min_ = -1 , T max_ = -1  );
+	Persistence_heat_maps()
+	{
+		this->scalling_function_with_respect_to_distance_from_diagonal = constant_function;
+		this->erase_below_diagonal = false;
+		this->filter = create_Gaussian_filter(5,1);
+		this->min_ = this->max_ = 0;
+	};
+    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)( std::pair< double , double >& point_in_diagram ) = constant_function, bool erase_below_diagonal = false , int number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  );
+    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)( std::pair< double , double >& point_in_diagram ) = constant_function, bool erase_below_diagonal = false , int number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  );
 
 	
 	void load_mean( const std::vector<Persistence_heat_maps*>& maps );
@@ -85,44 +109,37 @@ public:
     
     void plot( const char* filename );
 private:
-	std::vector< std::vector<T> > check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps );
-    void construct( const Persistence_intervals<T>& intervals_ , T min_ = -1 , T max_ = -1);
-    const Persistence_intervals<T>& intervals;
-    int number_of_pixels;
-    //int diameter_of_gaussian;
-    //we should keep here the Gaussian filter
-    
-    //here we should have the scalling function
+	//private methods
+	std::vector< std::vector<double> > check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps );
+    void construct( const std::vector< std::pair<double,double> >& intervals_ , double min_ = -1 , double max_ = -1);
     
-    //here we should have the info if we cut out everything below diagonal or not
-
-    T min_;
-    T max_;
-    std::vector< std::vector<T> > heat_map;    
+    //data
+    std::vector< std::vector<double> > filter;      
+    double (*scalling_function_with_respect_to_distance_from_diagonal)( std::pair< double , double >& point_in_diagram );    
+    bool erase_below_diagonal;
+    double min_;
+    double max_;
+    std::vector< std::vector< double > > heat_map;    
 };
 
 
-template <typename T>
-void Persistence_heat_maps<T>::construct( std::vector< Persistence_intervals<T>* > intervals_ , T min_ , T max_ )
+//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,doubl> >& intervals_ , double min_ , double max_ )
 {
     bool dbg = false;
 
-    for ( size_t i = 0 ; i != intervals_.size() ; ++i )
-    {
-        this->intervals.push_back( intervals_[i] );
-    }
-
     if ( min_ == max_ )
     {
         //in this case, we want the program to set up the min_ and max_ values by itself.
         min_ = INT_MAX;
         max_ = -INT_MAX;
+        
+        
         for ( size_t i = 0 ; i != intervals_.size() ; ++i )
         {
-            std::pair<T,T> min_max = intervals_[i]->min_max();
-            if ( min_max.first < min_ )min_ = min_max.first;
-            if ( min_max.second > max_ )max_ = min_max.second;
-        }
+			if ( intervals_[i].first < min_ )min_ = min_max.first;
+			if ( intervals_[i].second > max_ )max_ = min_max.second;
+		}
         //now we have the structure filled in, and moreover we know min_ and max_ values of the interval, so we know the range.
 
         //add some more space:
@@ -130,104 +147,53 @@ void Persistence_heat_maps<T>::construct( std::vector< Persistence_intervals<T>*
         max_ += fabs(max_ - min_)/100;
     }
 
-
-    cerr << "min_ : " << min_ << endl;
-    cerr << "max_ : " << max_ << endl;
+	if ( dng )
+	{
+		cerr << "min_ : " << min_ << endl;
+		cerr << "max_ : " << max_ << endl;
+	}
 
     this->min_ = min_;
     this->max_ = max_;
 
 
     //initialization of the structure heat_map
-    std::vector< std::vector< std::vector<T> > > heat_map_;
+    std::vector< std::vector<double> > heat_map_;
     for ( size_t i = 0 ; i != this->number_of_pixels ; ++i )
-    {
-        std::vector< std::vector<T> > row;
-        for ( size_t j = 0 ; j != this->number_of_pixels ; ++j )
-        {
-            std::vector<T> v( this->intervals.size() );
-            std::fill( v.begin() , v.end() , 0 );
-            row.push_back( v );
-        }
-        heat_map_.push_back( row );
+    {       
+		std::vector<double> v( this->number_of_pixels , 0 );
+        heat_map_.push_back( v );
     }
     this->heat_map = heat_map_;
 
     if (dbg)cerr << "Done creating of the heat map, now we will fill in the structure \n";
 
-    //and filling-in the structure:
-    //for every persistence diagram:
-    for ( size_t diag_no = 0 ; diag_no != this->intervals.size() ; ++diag_no )
-    {
-        //for every point in this diagram:
-        for ( size_t point_no = 0 ; point_no != this->intervals[diag_no]->intervals.size() ; ++point_no )
-        {
-            //take the coordinates of this->intervals[diag_no][point_no]
-            //this->heat_map[diag_no]
-            size_t x_position_matrix = ((this->intervals[diag_no]->intervals[point_no].first-min_)/( max_ - min_ )) * this->number_of_pixels;
-            size_t y_position_matrix = ((this->intervals[diag_no]->intervals[point_no].second-min_)/( max_ - min_ )) * this->number_of_pixels;
-
-            if ( (x_position_matrix > this->number_of_pixels) || (y_position_matrix > this->number_of_pixels) )
-            {
-                cerr << "Points coords out of range : (" << this->intervals[diag_no]->intervals[point_no].first << " , " << this->intervals[diag_no]->intervals[point_no].second << "). This is probably due to wrong values of min and max being setup. This point will be disregarded \n";
-                continue;
-            }
-
-            if ( dbg )
-            {
-                cerr << "The coords of the point are : (" << this->intervals[diag_no]->intervals[point_no].first << " , " << this->intervals[diag_no]->intervals[point_no].second << ") \n";
-                cerr << "The matrix coords are : (" << x_position_matrix << " , " << y_position_matrix << ") \n";
-                getchar();
-            }
-
-            if ( x_position_matrix == y_position_matrix )continue;
-
-            size_t lower_side = (size_t)std::max( (int)(y_position_matrix-this->diameter_of_gaussian) , (int)0 );
-            size_t upper_side = (size_t)std::min( (int)(y_position_matrix+this->diameter_of_gaussian) , (int)(this->number_of_pixels-1) );
-            size_t left_side = (size_t)std::max(  (int)(x_position_matrix - this->diameter_of_gaussian) , (int)0 );
-            size_t right_side = (size_t)std::min( (int)(x_position_matrix + this->diameter_of_gaussian) , (int)(this->number_of_pixels-1) );
-
-            if ( dbg )
-            {
-                cerr << "lower_side : " << lower_side << endl;
-                cerr << "upper_side : " << upper_side << endl;
-                cerr << "left_side : " << left_side << endl;
-                cerr <<  "right_side : " << right_side << endl;
-            }
+	//we will use this extra data structure to store the information about points in the diagram and their weights with respect to the weighting function.
+	std::vector< std::vector< std::vector< double > > > weights_of_points;
+	for ( aaa )
+	
 
-            for ( size_t y = lower_side ; y <= upper_side ; ++y )
-            {
-                for ( size_t x = left_side ; x <= right_side ; ++x )
-                {
-                    //cerr << "(" << x << "," << y << ")";
-                    //make it more ambitious:
-                    this->heat_map[x][y][diag_no]++;
-                }
-                //cerr << endl;
-            }
-            //getchar();
-        }
-    }
+	for ( size_t i = 0 ; i != this->number_of_pixels ; ++i )
+	{
+		for ( size_t j = 0 ; j != this->number_of_pixels ; ++j )
+		{
+			//compute the value of a heat map at a point [i,j].
+			aaa, todo
+		}
+	}
 }//construct
 
-template <typename T>
-Persistence_heat_maps<T>::Persistence_heat_maps( std::vector< Persistence_intervals<T>* > intervals_ , int number_of_pixels_ , int diameter_of_gaussian_ , T min_ , T max_ ):number_of_pixels(number_of_pixels_),diameter_of_gaussian(diameter_of_gaussian_)
+todo
+Persistence_heat_maps::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)( std::pair< double , double >& point_in_diagram ) = constant_function, bool erase_below_diagonal = false , int number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  )
 {
     this->construct( intervals_ , min_ , max_ );
 }
 
-template <typename T>
-Persistence_heat_maps<T>::Persistence_heat_maps( char* name_of_file_with_names_of_files_with_intervals , int number_of_pixels_ , int diameter_of_gaussian_ , T min_ , T max_):number_of_pixels(number_of_pixels_),diameter_of_gaussian(diameter_of_gaussian_)
-{
-    std::vector< Persistence_intervals<T>* > intervals;
-    std::vector<std::string> names = readFileNames( name_of_file_with_names_of_files_with_intervals );
-    for ( size_t file_no = 0 ; file_no != names.size() ; ++file_no )
-    {
-        cout << "Reading file : " << names[file_no] << endl;
-        Persistence_intervals<T>* interval = new Persistence_intervals<T>( (char*)names[file_no].c_str() );
-        intervals.push_back( interval );
-    }
-    this->construct( intervals , min_ , max_);
+todo
+Persistence_heat_maps::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)( std::pair< double , double >& point_in_diagram ) = constant_function, bool erase_below_diagonal = false , int number_of_pixels = 1000 , double min_ = -1 , double max_ = -1  )
+{    
+    std::vector< std::pair< double , double > > interval = read_standard_file( name_of_file_with_names_of_files_with_interval );   
+    this->construct( interval , min_ , max_);
 }
 
 
@@ -248,8 +214,8 @@ double vector_median(std::vector<double> vec)
     }
 }
 
-template <typename T>
-std::vector< std::vector<T> > Persistence_heat_maps<T>::check_and_initialize_maps( const std::vector<Persistence_heat_maps*>& maps )
+
+std::vector< std::vector<double> > Persistence_heat_maps::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.heat_maps.size() ; ++i )
@@ -265,21 +231,21 @@ std::vector< std::vector<T> > Persistence_heat_maps<T>::check_and_initialize_map
 			throw "Sizes of Persistence_heat_maps are not compatible. The program will terminate now \n";			
 		}
 	} 	
-	std::vector< std::vector<T> > heat_maps( maps.heat_maps.size() );
+	std::vector< std::vector<double> > heat_maps( maps.heat_maps.size() );
 	for ( size_t i = 0 ; i != maps.heat_maps.size() ; ++i )
     {
-		 std::vector<T> v( maps.heat_maps[0].size() , 0 );
+		 std::vector<double> v( maps.heat_maps[0].size() , 0 );
 		 heat_maps[i] = v;
 	}
 	return heat_maps;
 }
 
-template <typename T>
-void Persistence_heat_maps<T>::load_median( const std::vector<Persistence_heat_maps*>& maps )
+
+void Persistence_heat_maps::load_median( const std::vector<Persistence_heat_maps*>& maps )
 {
-	std::vector< std::vector<T> > heat_maps = this->check_and_initialize_maps( maps );
+	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
 	
-	std::vector<T> to_compute_median( maps.size() );
+	std::vector<double> to_compute_median( maps.size() );
     for ( size_t i = 0 ; i != heat_map.size() ; ++i )
     {
         for ( size_t j = 0 ; j != heat_map[i].size() ; ++j )
@@ -295,22 +261,22 @@ void Persistence_heat_maps<T>::load_median( const std::vector<Persistence_heat_m
 }
 
 
-template <typename T>
-void Persistence_heat_maps<T>::load_mean( const std::vector<Persistence_heat_maps*>& maps )
+
+void Persistence_heat_maps::load_mean( const std::vector<Persistence_heat_maps*>& maps )
 {
-	std::vector< std::vector<T> > heat_maps = this->check_and_initialize_maps( maps );
+	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
 	
-	std::vector<T> to_compute_median( maps.size() );
+	std::vector<double> to_compute_median( maps.size() );
     for ( size_t i = 0 ; i != heat_map.size() ; ++i )
     {
         for ( size_t j = 0 ; j != heat_map[i].size() ; ++j )
         {
-			T mean = 0;
+			double mean = 0;
 			for ( size_t aa = 0 ; aa != maps.size() ; ++aa )
 			{
 				mean += maps[i][j];
 			}
-            heat_maps[i][j] =  mean/(T)maps.size();
+            heat_maps[i][j] =  mean/(double)maps.size();
         }       
     }
 	this->heat_map = heat_maps;
@@ -319,10 +285,10 @@ void Persistence_heat_maps<T>::load_mean( const std::vector<Persistence_heat_map
 
 
 
-template <typename T>
-void Persistence_heat_maps<T>::load_percentage_of_active( const std::vector<Persistence_heat_maps*> maps , size_t cutoff )
+
+void Persistence_heat_maps::load_percentage_of_active( const std::vector<Persistence_heat_maps*> maps , size_t cutoff )
 {
-	std::vector< std::vector<T> > heat_maps = this->check_and_initialize_maps( maps );
+	std::vector< std::vector<double> > heat_maps = this->check_and_initialize_maps( maps );
 
     for ( size_t i = 0 ; i != this->heat_map.size() ; ++i )
     {
@@ -343,13 +309,12 @@ void Persistence_heat_maps<T>::load_percentage_of_active( const std::vector<Pers
             }
         }       
     }
-
     this->heat_map = heat_maps;
 }
 
 
-template <typename T>
-void Persistence_heat_maps<T>::plot( const char* filename )
+
+void Persistence_heat_maps::plot( const char* filename )
 {
 	 ofstream out;
     out.open( filename );
@@ -367,8 +332,8 @@ void Persistence_heat_maps<T>::plot( const char* filename )
 }
 
 
-template <typename T>
-void Persistence_heat_maps<T>::write_to_file( const char* filename )
+
+void Persistence_heat_maps::write_to_file( const char* filename )
 {
 	ofstream out;
 	out.open( filename );
@@ -383,8 +348,8 @@ void Persistence_heat_maps<T>::write_to_file( const char* filename )
 	out.close();
 }
 
-template <typename T>
-Persistence_heat_maps<T>::Persistence_heat_maps( const char* filename )
+
+Persistence_heat_maps::Persistence_heat_maps( const char* filename )
 {
 	bool dbg = true;
 	
@@ -407,10 +372,10 @@ Persistence_heat_maps<T>::Persistence_heat_maps( const char* filename )
         std::stringstream lineSS;
         lineSS << line;
         
-        std::vector< std::vector<T> > line_of_heat_map;    
+        std::vector< std::vector<double> > line_of_heat_map;    
         while ( lineSS.good() )
         {
-			T point;
+			double point;
 			lineSS >> point;
 			line_of_heat_map.push_back( point );
 			if ( dbg )
diff --git a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h
index ecb69a2f..e0db0070 100644
--- a/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h
+++ b/src/Gudhi_stat/include/gudhi/concretizations/Persistence_intervals.h
@@ -608,8 +608,8 @@ std::vector< double > Persistence_intervals::k_n_n( size_t k , size_t where_to_c
 

 std::pair<double,double> Persistence_intervals::min_max()

 {

-    double min_ = INT_MAX;

-    double max_ = -INT_MAX;

+    double min_ = std::numeric_limits<int>::max();

+    double max_ = -std::numeric_limits<int>::max();

     for ( size_t i = 0 ; i != this->intervals.size() ; ++i )

     {

         if ( this->intervals[i].first < min_ )min_ = this->intervals[i].first;

diff --git a/src/Gudhi_stat/test/CMakeLists.txt b/src/Gudhi_stat/test/CMakeLists.txt
index 7975e68a..6d728fd8 100644
--- a/src/Gudhi_stat/test/CMakeLists.txt
+++ b/src/Gudhi_stat/test/CMakeLists.txt
@@ -58,3 +58,15 @@ add_test(NAME persistence_lanscapes_test
   COMMAND ${CMAKE_CURRENT_BINARY_DIR}/persistence_lanscapes_test
     # XML format for Jenkins xUnit plugin 
     --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/persistence_lanscapes_test.xml --log_level=test_suite --report_level=no)
+
+
+
+
+add_executable ( persistence_lanscapes_on_grid_test persistence_lanscapes_on_grid_test.cpp )
+target_link_libraries(persistence_lanscapes_on_grid_test ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+# Unitary tests
+add_test(NAME persistence_lanscapes_on_grid_test
+  COMMAND ${CMAKE_CURRENT_BINARY_DIR}/persistence_lanscapes_on_grid_test
+    # XML format for Jenkins xUnit plugin 
+    --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/persistence_lanscapes_on_grid_test.xml --log_level=test_suite --report_level=no)