a few more changes

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

Former-commit-id: d12b2413620003883e82b3022f8d933aa05e2856

3 files changed, 220 insertions, 8 deletions

diff --git a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp
index 41cf915a..6f12dada 100644
--- a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp
+++ b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp
@@ -24,6 +24,7 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Persistent_cohomology.h>
 #include <gudhi/reader_utils.h>
+#include <gudhi/file_writer.h>
 
 #include <boost/program_options.hpp>
 
@@ -38,6 +39,7 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
 using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
 using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
 using Correlation_matrix = std::vector<std::vector<Filtration_value>>;
+using intervals_common = Gudhi::Persistence_interval_common< double , int >;
 
 void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag,
                      Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence);
@@ -66,6 +68,13 @@ int main(int argc, char* argv[]) {
     }
   }
 
+  //If the treshold, being minimal corelation is in the range [0,1], 
+  //change it to 1-threshold
+  if ( ( threshold>=0 ) && ( threshold<=1 ) )
+  {
+	  threshold = 1-threshold;
+  }
+  
   Rips_complex rips_complex_from_file(correlations, threshold);
 
   // Construct the Rips complex in a Simplex Tree
@@ -82,15 +91,34 @@ int main(int argc, char* argv[]) {
   Persistent_cohomology pcoh(simplex_tree);
   // initializes the coefficient field for homology
   pcoh.init_coefficients(p);
-
+  //compute persistence
   pcoh.compute_persistent_cohomology(min_persistence);
-
-  // Output the diagram in filediag
+    
+  
+  //invert the persistence diagram 
+  auto pairs = pcoh.get_persistent_pairs();
+  std::vector< intervals_common > processed_persistence_intervals;
+  processed_persistence_intervals.reserve( pairs.size() );
+  for (auto pair :pairs )
+  {
+	  double birth = 1-simplex_tree.filtration( get<0>(pair) );
+	  double death = 1-simplex_tree.filtration( get<1>(pair) );
+	  unsigned dimension = (unsigned)simplex_tree.dimension( get<0>(pair) );
+	  int field = get<2>(pair);
+	  processed_persistence_intervals.push_back(
+	  intervals_common(birth, death,dimension,field)
+	  );
+  }	  
+  
+  //sort the processed intervals:
+  std::sort( processed_persistence_intervals.begin() , processed_persistence_intervals.end() );
+
+  //and write them to a file 
   if (filediag.empty()) {
-    pcoh.output_diagram();
+    write_persistence_intervals_to_stream(processed_persistence_intervals);
   } else {
     std::ofstream out(filediag);
-    pcoh.output_diagram(out);
+    write_persistence_intervals_to_stream(processed_persistence_intervals,out);
     out.close();
   }
   return 0;
@@ -103,6 +131,9 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::
   hidden.add_options()(
       "input-file", po::value<std::string>(&csv_matrix_file),
       "Name of file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'.");
+  hidden.add_options()
+      ("input-file", po::value<std::string>(&csv_matrix_file),
+       "Name of file containing a corelation matrix. Can be square or lower triangular matrix. Separator is ';'.");
 
   po::options_description visible("Allowed options", 100);
   visible.add_options()("help,h", "produce help message")(
@@ -118,6 +149,19 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::
       "min-persistence,m", po::value<Filtration_value>(&min_persistence),
       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
       "intervals");
+  visible.add_options()
+      ("help,h", "produce help message")
+      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+       "Name of file in which the persistence diagram is written. Default print in std::cout")
+      ("min-edge-corelation,c",
+       po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+       "Minimal corelation of an edge for the Rips complex construction.")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+       "Maximal dimension of the Rips complex we want to compute.")
+      ("field-charac,p", po::value<int>(&p)->default_value(11),
+       "Characteristic p of the coefficient field Z/pZ for computing homology.")
+      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
+       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
 
   po::positional_options_description pos;
   pos.add("input-file", 1);
@@ -132,7 +176,7 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::
   if (vm.count("help") || !vm.count("input-file")) {
     std::cout << std::endl;
     std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of a Rips complex defined on a set of distance matrix.\n \n";
+    std::cout << "of a Rips complex defined on a corelation matrix.\n \n";
     std::cout << "The output diagram contains one bar per line, written with the convention: \n";
     std::cout << "   p   dim b d \n";
     std::cout << "where dim is the dimension of the homological feature,\n";
diff --git a/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp
index ae347a00..d1ccbf31 100644
--- a/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp
+++ b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp
@@ -42,17 +42,28 @@ int main() {
   }
 
   //-----------------------------------------------------------------------------
-  // Now the correlation matrix is really the distance matrix and can be processed further.
+  // Now the correlation matrix is a distance matrix and can be processed further. 
   //-----------------------------------------------------------------------------
   Distance_matrix distances = correlations;
 
+  //------------------------------------------------------------------------------
+  //Note that this treshold mean that the points in the distance 1, i.e. corelation
+  //0 will be connected. 
+  //------------------------------------------------------------------------------
   double threshold = 1.0;
+  
+  
   Rips_complex rips_complex_from_points(distances, threshold);
 
   Simplex_tree stree;
   rips_complex_from_points.create_complex(stree, 1);
   // ----------------------------------------------------------------------------
-  // Display information about the one skeleton Rips complex
+  // Display information about the one skeleton Rips complex. Note that
+  // the filtration displayed here comes from the distance matrix computed
+  // above, which is 1 - initial correlation matrix. Only this way, we obtain
+  // a complex with filtration. If a correlation matrix is used instead, we would
+  // have a reverse filtration (i.e. filtration of boundary of each simplex S
+  // is greater or equal to the filtration of S). 
   // ----------------------------------------------------------------------------
   std::cout << "Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - "
             << stree.num_vertices() << " vertices." << std::endl;
diff --git a/src/common/include/gudhi/file_writer.h b/src/common/include/gudhi/file_writer.h
new file mode 100644
index 00000000..1b59ae46
--- /dev/null
+++ b/src/common/include/gudhi/file_writer.h
@@ -0,0 +1,157 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2017  Swansea University, UK
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef FILE_WRITER_
+#define FILE_WRITER_
+
+#include <iostream>
+#include <string>
+#include <limits>
+
+namespace Gudhi {
+
+
+/**
+* This is a class to store persistence intervals. Its main purpose is to 
+* exchange data in between different packages and provide unified way
+* of writing a collection of persistence intervals to file. 
+**/
+template <typename Filtration_type , typename Coefficient_field>
+class Persistence_interval_common
+{
+public:	
+	Persistence_interval_common( Filtration_type birth , Filtration_type death ):
+	birth_(birth),death_(death),dimension_(std::numeric_limits<unsigned>::max),
+	Arith_element_(std::numeric_limits<Coefficient_field>::max() ){} 
+	
+	Persistence_interval_common( Filtration_type birth , Filtration_type death,
+	unsigned dim ):
+	birth_(birth),death_(death),dimension_(dim),
+	Arith_element_(std::numeric_limits<Coefficient_field>::max()){}
+	
+	Persistence_interval_common( Filtration_type birth , Filtration_type death,
+	unsigned dim , Coefficient_field field ):
+	birth_(birth),death_(death),dimension_(dim),
+	Arith_element_(field){}
+	
+	
+	inline bool operator == ( const Persistence_interval_common &i2)
+	{
+		return ( 
+		        (this->birth_ == i2.birth_) && (this->death_ == i2.death_) && 
+		        (this->dimension_ == i2.dimension_) && (this->Arith_element_ == i2.Arith_element_)
+		       );
+	}
+	
+	inline bool operator != ( const Persistence_interval_common &i2)
+	{
+		return (!((*this)==i2));
+	}
+	
+	
+	/**
+	 * Note that this operator do not take Arith_element into account when doing comparisions.
+	**/ 
+	inline bool operator < ( const Persistence_interval_common &i2)
+	{				
+		if ( this->birth_ < i2.birth_ )
+		{
+			return true;
+		}
+		else
+		{
+			if ( this->birth_ > i2.birth_ )
+			{
+				return false;
+			}
+			else
+			{
+				//in this case this->birth_ == i2.birth_
+				if ( this->death_ > i2.death_ )
+				{
+					return true;
+				}
+				else
+				{
+					if ( this->death_ < i2.death_ )
+					{
+						return false;
+					}
+					else
+					{
+						//in this case this->death_ == i2.death_
+						if ( this->dimension_ < i2.dimension_ )
+						{
+							return true;
+						}
+						else
+						{
+							//in this case this->dimension >= i2.dimension
+							return false;
+						}
+					}
+				}
+			}
+		}
+	}
+	
+	friend std::ostream& operator<<(std::ostream& out, const Persistence_interval_common& it)
+	{
+		if ( it.Arith_element_ != std::numeric_limits<Coefficient_field>::max() ) 
+		{
+			out << it.Arith_element_ << "  ";	
+		}		
+		if ( it.dimension_ != std::numeric_limits<unsigned>::max() )
+		{
+			out << it.dimension_ << " ";
+		}				
+		out << it.birth_ << " " << it.death_ << " ";				
+		return out;
+	}
+	
+private:
+	Filtration_type birth_;
+	Filtration_type death_;
+	unsigned dimension_;
+	Coefficient_field Arith_element_;
+};//Persistence_interval_common
+
+
+/**
+ * This function write a vector<Persistence_interval_common> to a stream
+**/ 
+template <typename Filtration_type , typename Coefficient_field>
+void write_persistence_intervals_to_stream( 
+const std::vector< Persistence_interval_common<Filtration_type,Coefficient_field> >& intervals,
+                                          std::ostream& out = std::cout )
+{
+	for ( auto interval : intervals )
+	{
+		out << interval << std::endl;
+	}
+}//write_persistence_intervals_to_stream
+
+
+
+}
+
+#endif //FILE_WRITER_