From 0729a55a67503c068e4843c63930d6b29e76f7ac Mon Sep 17 00:00:00 2001 From: pdlotko Date: Fri, 3 Nov 2017 05:10:17 +0000 Subject: 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 --- .../rips_correlation_matrix_persistence.cpp | 56 +++++++- ...e_one_skeleton_rips_from_correlation_matrix.cpp | 15 +- src/common/include/gudhi/file_writer.h | 157 +++++++++++++++++++++ 3 files changed, 220 insertions(+), 8 deletions(-) create mode 100644 src/common/include/gudhi/file_writer.h 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 #include #include +#include #include @@ -38,6 +39,7 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; using Correlation_matrix = std::vector>; +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(&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(&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(&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(&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(&threshold)->default_value(std::numeric_limits::infinity()), + "Minimal corelation of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::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 . + */ + +#ifndef FILE_WRITER_ +#define FILE_WRITER_ + +#include +#include +#include + +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 +class Persistence_interval_common +{ +public: + Persistence_interval_common( Filtration_type birth , Filtration_type death ): + birth_(birth),death_(death),dimension_(std::numeric_limits::max), + Arith_element_(std::numeric_limits::max() ){} + + Persistence_interval_common( Filtration_type birth , Filtration_type death, + unsigned dim ): + birth_(birth),death_(death),dimension_(dim), + Arith_element_(std::numeric_limits::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::max() ) + { + out << it.Arith_element_ << " "; + } + if ( it.dimension_ != std::numeric_limits::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 to a stream +**/ +template +void write_persistence_intervals_to_stream( +const std::vector< Persistence_interval_common >& intervals, + std::ostream& out = std::cout ) +{ + for ( auto interval : intervals ) + { + out << interval << std::endl; + } +}//write_persistence_intervals_to_stream + + + +} + +#endif //FILE_WRITER_ -- cgit v1.2.3