From 1d611ca664d78d3b0b53264c52b04a8396383140 Mon Sep 17 00:00:00 2001 From: pdlotko Date: Mon, 25 Sep 2017 20:14:45 +0000 Subject: Adding all that is needed for Rips complex construction from correlation matrix. Moduo the fact that I cannot see the doc (althoug make doxygen do not create any errors), this is ready for review. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@2712 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: e8e7acceaf2be6f17baa3d77c85e207a74a4fb33 --- .../doc/Intro_persistent_cohomology.h | 13 +- src/Persistent_cohomology/example/CMakeLists.txt | 6 + .../rips_correlation_matrix_persistence.cpp | 160 +++++++++++++++++++++ src/Rips_complex/doc/Intro_rips_complex.h | 23 +++ src/Rips_complex/example/CMakeLists.txt | 4 + ...eleton_rips_from_correlation_matrix_for_doc.txt | 17 +++ 6 files changed, 222 insertions(+), 1 deletion(-) create mode 100644 src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp create mode 100644 src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt (limited to 'src') diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index e17e5926..1fe048bc 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -160,7 +160,18 @@ persistence diagram with a family of field coefficients. \li Persistent_cohomology/rips_distance_matrix_persistence.cpp computes the Rips complex of a distance matrix and -outputs its persistence diagram. +outputs its persistence diagram. The file should contain lower diagonal distance matrix with semicolons as separators. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. + +\li +Persistent_cohomology/rips_correlation_matrix_persistence.cpp +computes the Rips complex of a correlation matrix and +outputs its persistence diagram. Note that no check is performed if +the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. The +input is to be given either as a lower triangular matrix. +Please refer to data/distance_matrix/lower_triangular_correlation_matrix.csv for an example of a file. \li Persistent_cohomology/alpha_complex_3d_persistence.cpp computes the persistent homology with diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index f47de4c3..8a21d038 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -8,6 +8,9 @@ add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simp add_executable(rips_distance_matrix_persistence rips_distance_matrix_persistence.cpp) target_link_libraries(rips_distance_matrix_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) +add_executable(rips_correlation_matrix_persistence rips_correlation_matrix_persistence.cpp) +target_link_libraries(rips_correlation_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable(rips_persistence rips_persistence.cpp) target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) @@ -24,6 +27,7 @@ if (TBB_FOUND) target_link_libraries(plain_homology ${TBB_LIBRARIES}) target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES}) target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES}) + target_link_libraries(rips_correlation_matrix_persistence ${TBB_LIBRARIES}) target_link_libraries(rips_persistence ${TBB_LIBRARIES}) target_link_libraries(rips_persistence_step_by_step ${TBB_LIBRARIES}) target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_LIBRARIES}) @@ -35,6 +39,8 @@ add_test(NAME Persistent_cohomology_example_from_simple_simplex_tree COMMAND $ "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0") +add_test(rips_distance_matrix ${CMAKE_CURRENT_BINARY_DIR}/rips_distance_matrix_persistence + ${CMAKE_SOURCE_DIR}/data/distance_matrix/full_correlation_matrix.csv.csv -r 1.0 -d 3 -p 3 -m 0) add_test(NAME Persistent_cohomology_example_from_rips_on_tore_3D COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") add_test(NAME Persistent_cohomology_example_from_rips_step_by_step_on_tore_3D COMMAND $ diff --git a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp new file mode 100644 index 00000000..6f2891fe --- /dev/null +++ b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp @@ -0,0 +1,160 @@ +/* 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, Vincent Rouvreau + * + * Copyright (C) 2016 INRIA + * + * 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 . + */ + +#include +#include +#include +#include + +#include + +#include +#include +#include // infinity + + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +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>; + +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); + +int main(int argc, char * argv[]) { + std::string csv_matrix_file; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); + + Correlation_matrix correlations = Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); + + //Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: + for ( size_t i = 0 ; i != correlations.size() ; ++i ) + { + for ( size_t j = 0 ; j != correlations[i].size() ; ++j ) + { + correlations[i][j] = 1-correlations[i][j]; + if ( correlations[i][j] < 0 ) + { + std::cerr << "The input matrix is not a correlation matrix. \n"; + throw "The input matrix is not a correlation matrix. \n"; + } + } + } + + Rips_complex rips_complex_from_file(correlations, threshold); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + rips_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + return 0; +} + +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) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + 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 ';'."); + + po::options_description visible("Allowed options", 100); + 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") + ("max-edge-length,r", + po::value(&threshold)->default_value(std::numeric_limits::infinity()), + "Maximal length 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); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + 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 << "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"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 124dfec9..401b689b 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -146,6 +146,29 @@ namespace rips_complex { * * \include Rips_complex/full_skeleton_rips_for_doc.txt * + * + * \section ripscorrelationematrix Correlation matrix + * + * Analogously to the case of distance matrix, Rips complexes can be also constructed based on correlation matrix. + * Given a correlation matrix M, comportment-wise 1-M is a distance matrix. + * This example builds the one skeleton graph from the given corelation matrix and threshold value. + * Then it creates a `Simplex_tree` with it. + * + * Then, it is asked to display information about the simplicial complex. + * + * \include Rips_complex/example_one_skeleton_rips_from_correlation_matrix.cpp + * + * When launching: + * + * \code $> ./example_one_skeleton_from_correlation_matrix + * \endcode + * + * the program output is: + * + * \include Rips_complex/one_skeleton_rips_from_correlation_matrix_for_doc.txt + * + * All the other constructions discussed for Rips complex for distance matrix can be also performed for Rips complexes construction from correlation matrices. + * * \copyright GNU General Public License v3. * \verbatim Contact: gudhi-users@lists.gforge.inria.fr \endverbatim */ diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt index 2940f164..f58ab455 100644 --- a/src/Rips_complex/example/CMakeLists.txt +++ b/src/Rips_complex/example/CMakeLists.txt @@ -9,12 +9,16 @@ add_executable ( Rips_complex_example_one_skeleton_from_points example_one_skele # Distance matrix add_executable ( Rips_complex_example_one_skeleton_from_distance_matrix example_one_skeleton_rips_from_distance_matrix.cpp ) +add_executable ( example_one_skeleton_rips_from_correlation_matrix example_one_skeleton_rips_from_correlation_matrix.cpp ) + + add_executable ( Rips_complex_example_from_csv_distance_matrix example_rips_complex_from_csv_distance_matrix_file.cpp ) if (TBB_FOUND) target_link_libraries(Rips_complex_example_from_off ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_points ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${TBB_LIBRARIES}) + target_link_libraries(example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${TBB_LIBRARIES}) endif() diff --git a/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt b/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt new file mode 100644 index 00000000..640d7083 --- /dev/null +++ b/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt @@ -0,0 +1,17 @@ +Rips complex is of dimension 1 - 15 simplices - 5 vertices. +Iterator on Rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 4 0 ) -> [0.11] + ( 2 1 ) -> [0.26] + ( 3 2 ) -> [0.28] + ( 4 3 ) -> [0.3] + ( 4 1 ) -> [0.39] + ( 2 0 ) -> [0.77] + ( 1 0 ) -> [0.94] + ( 4 2 ) -> [0.97] + ( 3 0 ) -> [0.99] + ( 3 1 ) -> [0.99] -- cgit v1.2.3 From f48ae2be7dcd552f44723b3a6b02cc9aae361aca Mon Sep 17 00:00:00 2001 From: pdlotko Date: Mon, 25 Sep 2017 20:28:02 +0000 Subject: a missing file... git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@2713 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: ecdf0da96de01e37758e157d84aee13c8825dc0e --- ...e_one_skeleton_rips_from_correlation_matrix.cpp | 77 ++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp (limited to 'src') 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 new file mode 100644 index 00000000..0acdfe83 --- /dev/null +++ b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp @@ -0,0 +1,77 @@ +#include +#include +#include + +#include +#include +#include +#include // for std::numeric_limits + +int main() { + // Type definitions + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex; + using Distance_matrix = std::vector>; + + // User defined correlation matrix is: + // |1 0.06 0.23 0.01 0.89| + // |0.06 1 0.74 0.01 0.61| + // |0.23 0.74 1 0.72 0.03| + // |0.01 0.01 0.72 1 0.7 | + // |0.89 0.61 0.03 0.7 1 | + + + Distance_matrix correlations; + correlations.push_back({}); + correlations.push_back({0.06}); + correlations.push_back({0.23, 0.74}); + correlations.push_back({0.01, 0.01, 0.72}); + correlations.push_back({0.89, 0.61, 0.03, 0.7}); + + // ---------------------------------------------------------------------------- + // Convert correlation matrix to a distance matrix: + // ---------------------------------------------------------------------------- + for ( size_t i = 0 ; i != correlations.size() ; ++i ) + { + for ( size_t j = 0 ; j != correlations[i].size() ; ++j ) + { + correlations[i][j] = 1-correlations[i][j]; + if ( correlations[i][j] < 0 ) + { + std::cerr << "The input matrix is not a correlation matrix. \n"; + throw "The input matrix is not a correlation matrix. \n"; + } + } + } + + //----------------------------------------------------------------------------- + // Now the correlation matrix is really the distance matrix and can be processed further. + //----------------------------------------------------------------------------- + Distance_matrix distances = correlations; + + 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 + // ---------------------------------------------------------------------------- + std::cout << "Rips complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + + return 0; +} -- cgit v1.2.3 From ef61b085afd77976a2c7fc5dfa13bc4b293b4f95 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 28 Sep 2017 13:43:58 +0000 Subject: Remove python rips_complex construction from files as it can lead to errors with correlation matrix Add examples for doxygen Cythonization of rips correlation matrix git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@2727 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 8aae33839fa27f9d26897e625904671b2c05e0e7 --- src/Persistent_cohomology/example/CMakeLists.txt | 2 +- .../rips_correlation_matrix_persistence.cpp | 92 +++++++++------------- src/Rips_complex/example/CMakeLists.txt | 11 ++- ...e_one_skeleton_rips_from_correlation_matrix.cpp | 39 ++++----- src/common/doc/main_page.h | 4 + src/cython/cython/off_reader.pyx | 1 + src/cython/cython/rips_complex.pyx | 35 +------- .../doc/persistence_graphical_tools_user.rst | 8 +- src/cython/doc/pyplots/diagram_persistence.py | 5 +- src/cython/doc/rips_complex_user.rst | 73 ++++++++++++++++- .../alpha_rips_persistence_bottleneck_distance.py | 5 +- ...ersistence_from_distance_matrix_file_example.py | 3 +- ...ex_diagram_persistence_from_off_file_example.py | 3 +- src/cython/include/Rips_complex_interface.h | 17 ---- 14 files changed, 155 insertions(+), 143 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 8a21d038..926cef6b 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -40,7 +40,7 @@ add_test(NAME Persistent_cohomology_example_from_simple_simplex_tree COMMAND $ "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0") add_test(rips_distance_matrix ${CMAKE_CURRENT_BINARY_DIR}/rips_distance_matrix_persistence - ${CMAKE_SOURCE_DIR}/data/distance_matrix/full_correlation_matrix.csv.csv -r 1.0 -d 3 -p 3 -m 0) + ${CMAKE_SOURCE_DIR}/data/correlation_matrix/full_correlation_matrix.csv.csv -r 1.0 -d 3 -p 3 -m 0) add_test(NAME Persistent_cohomology_example_from_rips_on_tore_3D COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") add_test(NAME Persistent_cohomology_example_from_rips_step_by_step_on_tore_3D COMMAND $ diff --git a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp index 6f2891fe..41cf915a 100644 --- a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp @@ -1,5 +1,5 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ +/* 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, Vincent Rouvreau @@ -31,24 +31,18 @@ #include #include // infinity - // Types definition using Simplex_tree = Gudhi::Simplex_tree; using Filtration_value = Simplex_tree::Filtration_value; 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 Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; using Correlation_matrix = std::vector>; -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); +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); -int main(int argc, char * argv[]) { +int main(int argc, char* argv[]) { std::string csv_matrix_file; std::string filediag; Filtration_value threshold; @@ -58,22 +52,20 @@ int main(int argc, char * argv[]) { program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); - Correlation_matrix correlations = Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); - - //Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: - for ( size_t i = 0 ; i != correlations.size() ; ++i ) - { - for ( size_t j = 0 ; j != correlations[i].size() ; ++j ) - { - correlations[i][j] = 1-correlations[i][j]; - if ( correlations[i][j] < 0 ) - { - std::cerr << "The input matrix is not a correlation matrix. \n"; - throw "The input matrix is not a correlation matrix. \n"; - } - } - } - + Correlation_matrix correlations = + Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); + + // Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + correlations[i][j] = 1 - correlations[i][j]; + if (correlations[i][j] < 0) { + std::cerr << "The input matrix is not a correlation matrix. \n"; + throw "The input matrix is not a correlation matrix. \n"; + } + } + } + Rips_complex rips_complex_from_file(correlations, threshold); // Construct the Rips complex in a Simplex Tree @@ -104,33 +96,28 @@ int main(int argc, char * argv[]) { return 0; } -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) { +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) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); - 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 distance 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") - ("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") - ("max-edge-length,r", - po::value(&threshold)->default_value(std::numeric_limits::infinity()), - "Maximal length 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"); + 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")( + "max-edge-length,r", + po::value(&threshold)->default_value(std::numeric_limits::infinity()), + "Maximal length 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); @@ -139,8 +126,7 @@ void program_options(int argc, char * argv[] all.add(visible).add(hidden); po::variables_map vm; - po::store(po::command_line_parser(argc, argv). - options(all).positional(pos).run(), vm); + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); po::notify(vm); if (vm.count("help") || !vm.count("input-file")) { diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt index f58ab455..fcb1eaee 100644 --- a/src/Rips_complex/example/CMakeLists.txt +++ b/src/Rips_complex/example/CMakeLists.txt @@ -9,23 +9,25 @@ add_executable ( Rips_complex_example_one_skeleton_from_points example_one_skele # Distance matrix add_executable ( Rips_complex_example_one_skeleton_from_distance_matrix example_one_skeleton_rips_from_distance_matrix.cpp ) -add_executable ( example_one_skeleton_rips_from_correlation_matrix example_one_skeleton_rips_from_correlation_matrix.cpp ) - - add_executable ( Rips_complex_example_from_csv_distance_matrix example_rips_complex_from_csv_distance_matrix_file.cpp ) +# Correlation matrix +add_executable ( Rips_complex_example_one_skeleton_rips_from_correlation_matrix example_one_skeleton_rips_from_correlation_matrix.cpp ) + if (TBB_FOUND) target_link_libraries(Rips_complex_example_from_off ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_points ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${TBB_LIBRARIES}) - target_link_libraries(example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES}) endif() add_test(NAME Rips_complex_example_one_skeleton_from_points COMMAND $) add_test(NAME Rips_complex_example_one_skeleton_from_distance_matrix COMMAND $) +add_test(NAME Rips_complex_example_one_skeleton_rips_from_correlation_matrix + COMMAND $) add_test(NAME Rips_complex_example_from_off_doc_12_1 COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "12.0" "1" "${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt") @@ -61,3 +63,4 @@ install(TARGETS Rips_complex_example_from_off DESTINATION bin) install(TARGETS Rips_complex_example_one_skeleton_from_points DESTINATION bin) install(TARGETS Rips_complex_example_one_skeleton_from_distance_matrix DESTINATION bin) install(TARGETS Rips_complex_example_from_csv_distance_matrix DESTINATION bin) +install(TARGETS Rips_complex_example_one_skeleton_rips_from_correlation_matrix DESTINATION bin) 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 0acdfe83..ae347a00 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 @@ -21,7 +21,6 @@ int main() { // |0.01 0.01 0.72 1 0.7 | // |0.89 0.61 0.03 0.7 1 | - Distance_matrix correlations; correlations.push_back({}); correlations.push_back({0.06}); @@ -32,24 +31,21 @@ int main() { // ---------------------------------------------------------------------------- // Convert correlation matrix to a distance matrix: // ---------------------------------------------------------------------------- - for ( size_t i = 0 ; i != correlations.size() ; ++i ) - { - for ( size_t j = 0 ; j != correlations[i].size() ; ++j ) - { - correlations[i][j] = 1-correlations[i][j]; - if ( correlations[i][j] < 0 ) - { - std::cerr << "The input matrix is not a correlation matrix. \n"; - throw "The input matrix is not a correlation matrix. \n"; - } - } - } - + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + correlations[i][j] = 1 - correlations[i][j]; + if (correlations[i][j] < 0) { + std::cerr << "The input matrix is not a correlation matrix. \n"; + throw "The input matrix is not a correlation matrix. \n"; + } + } + } + //----------------------------------------------------------------------------- - // Now the correlation matrix is really the distance matrix and can be processed further. + // Now the correlation matrix is really the distance matrix and can be processed further. //----------------------------------------------------------------------------- Distance_matrix distances = correlations; - + double threshold = 1.0; Rips_complex rips_complex_from_points(distances, threshold); @@ -58,18 +54,17 @@ int main() { // ---------------------------------------------------------------------------- // Display information about the one skeleton Rips complex // ---------------------------------------------------------------------------- - std::cout << "Rips complex is of dimension " << stree.dimension() << - " - " << stree.num_simplices() << " simplices - " << - stree.num_vertices() << " vertices." << std::endl; + std::cout << "Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; - std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << - std::endl; + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << std::endl; for (auto f_simplex : stree.filtration_simplex_range()) { std::cout << " ( "; for (auto vertex : stree.simplex_vertex_range(f_simplex)) { std::cout << vertex << " "; } - std::cout << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + std::cout << ") -> " + << "[" << stree.filtration(f_simplex) << "] "; std::cout << std::endl; } diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 1a7994a5..91535ee6 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -456,11 +456,15 @@ make doxygen * @example Persistent_cohomology/persistence_from_simple_simplex_tree.cpp * @example Persistent_cohomology/plain_homology.cpp * @example Persistent_cohomology/rips_multifield_persistence.cpp + * @example Persistent_cohomology/rips_correlation_matrix_persistence.cpp * @example Persistent_cohomology/rips_distance_matrix_persistence.cpp * @example Persistent_cohomology/rips_persistence.cpp * @example Persistent_cohomology/custom_persistence_sort.cpp * @example Persistent_cohomology/rips_persistence_step_by_step.cpp + * @example Rips_complex/example_one_skeleton_rips_from_correlation_matrix.cpp + * @example Rips_complex/example_one_skeleton_rips_from_distance_matrix.cpp * @example Rips_complex/example_one_skeleton_rips_from_points.cpp + * @example Rips_complex/example_rips_complex_from_csv_distance_matrix_file.cpp * @example Rips_complex/example_rips_complex_from_off_file.cpp * @example Simplex_tree/mini_simplex_tree.cpp * @example Simplex_tree/simple_simplex_tree.cpp diff --git a/src/cython/cython/off_reader.pyx b/src/cython/cython/off_reader.pyx index b6e107ef..266dae2c 100644 --- a/src/cython/cython/off_reader.pyx +++ b/src/cython/cython/off_reader.pyx @@ -46,4 +46,5 @@ def read_off(off_file=''): return read_points_from_OFF_file(str.encode(off_file)) else: print("file " + off_file + " not found.") + return [] diff --git a/src/cython/cython/rips_complex.pyx b/src/cython/cython/rips_complex.pyx index ad9b0a4d..73b154b8 100644 --- a/src/cython/cython/rips_complex.pyx +++ b/src/cython/cython/rips_complex.pyx @@ -34,8 +34,6 @@ __license__ = "GPL v3" cdef extern from "Rips_complex_interface.h" namespace "Gudhi": cdef cppclass Rips_complex_interface "Gudhi::rips_complex::Rips_complex_interface": Rips_complex_interface(vector[vector[double]] values, double threshold, bool euclidean) - # bool from_file is a workaround for cython to find the correct signature - Rips_complex_interface(string file_name, double threshold, bool euclidean, bool from_file) void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, int dim_max) # RipsComplex python interface @@ -49,7 +47,7 @@ cdef class RipsComplex: cdef Rips_complex_interface * thisptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, off_file='', distance_matrix=None, csv_file='', max_edge_length=float('inf')): + def __init__(self, points=None, distance_matrix=None, max_edge_length=float('inf')): """RipsComplex constructor. :param max_edge_length: Rips value. @@ -60,41 +58,14 @@ cdef class RipsComplex: Or - :param off_file: An OFF file style name. - :type off_file: string - - Or - :param distance_matrix: A distance matrix (full square or lower triangular). :type points: list of list of double - - Or - - :param csv_file: A csv file style name containing a full square or a - lower triangular distance matrix. - :type csv_file: string """ # The real cython constructor - def __cinit__(self, points=None, off_file='', distance_matrix=None, csv_file='', max_edge_length=float('inf')): - if off_file is not '': - if os.path.isfile(off_file): - self.thisptr = new Rips_complex_interface(str.encode(off_file), - max_edge_length, - True, - True) - else: - print("file " + off_file + " not found.") - elif csv_file is not '': - if os.path.isfile(csv_file): - self.thisptr = new Rips_complex_interface(str.encode(csv_file), - max_edge_length, - False, - True) - else: - print("file " + csv_file + " not found.") - elif distance_matrix is not None: + def __cinit__(self, points=None, distance_matrix=None, max_edge_length=float('inf')): + if distance_matrix is not None: self.thisptr = new Rips_complex_interface(distance_matrix, max_edge_length, False) else: if points is None: diff --git a/src/cython/doc/persistence_graphical_tools_user.rst b/src/cython/doc/persistence_graphical_tools_user.rst index 9033331f..a5523d23 100644 --- a/src/cython/doc/persistence_graphical_tools_user.rst +++ b/src/cython/doc/persistence_graphical_tools_user.rst @@ -58,8 +58,8 @@ This function can display the persistence result as a diagram: import gudhi - rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_1307.off', max_edge_length=0.2) + point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + '/data/points/tore3D_1307.off') + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=0.2) simplex_tree = rips_complex.create_simplex_tree(max_dimension=3) diag = simplex_tree.persistence() plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13) @@ -69,8 +69,8 @@ This function can display the persistence result as a diagram: import gudhi - rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_1307.off', max_edge_length=0.2) + point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + '/data/points/tore3D_1307.off') + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=0.2) simplex_tree = rips_complex.create_simplex_tree(max_dimension=3) diag = simplex_tree.persistence() plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13) diff --git a/src/cython/doc/pyplots/diagram_persistence.py b/src/cython/doc/pyplots/diagram_persistence.py index c2fbf801..ac20bf47 100755 --- a/src/cython/doc/pyplots/diagram_persistence.py +++ b/src/cython/doc/pyplots/diagram_persistence.py @@ -1,7 +1,8 @@ import gudhi -rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_1307.off', max_edge_length=0.2) +point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + \ + '/data/points/tore3D_1307.off') +rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=0.2) simplex_tree = rips_complex.create_simplex_tree(max_dimension=3) diag = simplex_tree.persistence() plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13) diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst index 96ba9944..f0e7bf2d 100644 --- a/src/cython/doc/rips_complex_user.rst +++ b/src/cython/doc/rips_complex_user.rst @@ -101,8 +101,8 @@ Finally, it is asked to display information about the Rips complex. .. testcode:: import gudhi - rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/alphacomplexdoc.off', max_edge_length=12.0) + point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + '/data/points/alphacomplexdoc.off') + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=12.0) simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ repr(simplex_tree.num_simplices()) + ' simplices - ' + \ @@ -206,8 +206,9 @@ Finally, it is asked to display information about the Rips complex. .. testcode:: import gudhi - rips_complex = gudhi.RipsComplex(csv_file=gudhi.__root_source_dir__ + \ - '/data/distance_matrix/full_square_distance_matrix.csv', max_edge_length=12.0) + distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=gudhi.__root_source_dir__ + \ + '/data/distance_matrix/full_square_distance_matrix.csv') + rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=12.0) simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ repr(simplex_tree.num_simplices()) + ' simplices - ' + \ @@ -240,3 +241,67 @@ the program output is: [0, 3] -> 9.43 [4, 6] -> 9.49 [3, 6] -> 11.00 + +Correlation matrix +--------------- + +Example from a correlation matrix +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Analogously to the case of distance matrix, Rips complexes can be also constructed based on correlation matrix. +Given a correlation matrix M, comportment-wise 1-M is a distance matrix. +This example builds the one skeleton graph from the given corelation matrix and threshold value. +Then it creates a :doc:`Simplex_tree ` with it. + +Finally, it is asked to display information about the simplicial complex. + +.. testcode:: + + import gudhi + import numpy as np + + # User defined correlation matrix is: + # |1 0.06 0.23 0.01 0.89| + # |0.06 1 0.74 0.01 0.61| + # |0.23 0.74 1 0.72 0.03| + # |0.01 0.01 0.72 1 0.7 | + # |0.89 0.61 0.03 0.7 1 | + correlation_matrix=np.array([[1., 0.06, 0.23, 0.01, 0.89], + [0.06, 1., 0.74, 0.01, 0.61], + [0.23, 0.74, 1., 0.72, 0.03], + [0.01, 0.01, 0.72, 1., 0.7], + [0.89, 0.61, 0.03, 0.7, 1.]], float) + + distance_matrix = np.ones((correlation_matrix.shape),float) - correlation_matrix + rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=1.0) + + simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) + result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ + repr(simplex_tree.num_simplices()) + ' simplices - ' + \ + repr(simplex_tree.num_vertices()) + ' vertices.' + print(result_str) + fmt = '%s -> %.2f' + for filtered_value in simplex_tree.get_filtration(): + print(fmt % tuple(filtered_value)) + +When launching (Rips maximal distance between 2 points is 12.0, is expanded +until dimension 1 - one skeleton graph in other words), the output is: + +.. testoutput:: + + Rips complex is of dimension 1 - 15 simplices - 5 vertices. + [0] -> 0.00 + [1] -> 0.00 + [2] -> 0.00 + [3] -> 0.00 + [4] -> 0.00 + [0, 4] -> 0.11 + [1, 2] -> 0.26 + [2, 3] -> 0.28 + [3, 4] -> 0.30 + [1, 4] -> 0.39 + [0, 2] -> 0.77 + [0, 1] -> 0.94 + [2, 4] -> 0.97 + [0, 3] -> 0.99 + [1, 3] -> 0.99 diff --git a/src/cython/example/alpha_rips_persistence_bottleneck_distance.py b/src/cython/example/alpha_rips_persistence_bottleneck_distance.py index ab5fc1e9..386f8457 100755 --- a/src/cython/example/alpha_rips_persistence_bottleneck_distance.py +++ b/src/cython/example/alpha_rips_persistence_bottleneck_distance.py @@ -45,13 +45,14 @@ args = parser.parse_args() with open(args.file, 'r') as f: first_line = f.readline() if (first_line == 'OFF\n') or (first_line == 'nOFF\n'): + point_cloud = gudhi.read_off(off_file=args.file) print("#####################################################################") print("RipsComplex creation from points read in a OFF file") message = "RipsComplex with max_edge_length=" + repr(args.threshold) print(message) - rips_complex = gudhi.RipsComplex(off_file=args.file, + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=args.threshold) rips_stree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) @@ -67,7 +68,7 @@ with open(args.file, 'r') as f: message = "AlphaComplex with max_edge_length=" + repr(args.threshold) print(message) - alpha_complex = gudhi.AlphaComplex(off_file=args.file) + alpha_complex = gudhi.AlphaComplex(points=point_cloud) alpha_stree = alpha_complex.create_simplex_tree(max_alpha_square=(args.threshold * args.threshold)) message = "Number of simplices=" + repr(alpha_stree.num_simplices()) diff --git a/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py index 3baebd17..fa82a2f3 100755 --- a/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py +++ b/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py @@ -50,7 +50,8 @@ print("RipsComplex creation from distance matrix read in a csv file") message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) print(message) -rips_complex = gudhi.RipsComplex(csv_file=args.file, max_edge_length=args.max_edge_length) +distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=args.file) +rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=args.max_edge_length) simplex_tree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) message = "Number of simplices=" + repr(simplex_tree.num_simplices()) diff --git a/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py index 5951eedf..544b68c9 100755 --- a/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py +++ b/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py @@ -53,7 +53,8 @@ with open(args.file, 'r') as f: message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) print(message) - rips_complex = gudhi.RipsComplex(off_file=args.file, max_edge_length=args.max_edge_length) + point_cloud = gudhi.read_off(off_file=args.file) + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=args.max_edge_length) simplex_tree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) message = "Number of simplices=" + repr(simplex_tree.num_simplices()) diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 02985727..f26befbc 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -25,9 +25,7 @@ #include #include -#include #include -#include #include "Simplex_tree_interface.h" @@ -56,21 +54,6 @@ class Rips_complex_interface { } } - Rips_complex_interface(const std::string& file_name, double threshold, bool euclidean, bool from_file = true) { - if (euclidean) { - // Rips construction where file_name is an OFF file - Gudhi::Points_off_reader off_reader(file_name); - rips_complex_ = new Rips_complex::Filtration_value>(off_reader.get_point_cloud(), - threshold, - Gudhi::Euclidean_distance()); - } else { - // Rips construction where values is a distance matrix - Distance_matrix distances = - Gudhi::read_lower_triangular_matrix_from_csv_file::Filtration_value>(file_name); - rips_complex_ = new Rips_complex::Filtration_value>(distances, threshold); - } - } - ~Rips_complex_interface() { delete rips_complex_; } -- cgit v1.2.3 From af51e15aadd1c97da1767a2404fd1d60a105c191 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 28 Sep 2017 13:59:04 +0000 Subject: Fix some doc issue git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@2728 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 9a9745d18a37b53618aa3054114f26a8bdbc0009 --- .../doc/Intro_persistent_cohomology.h | 16 +++++++++------- src/cython/doc/rips_complex_user.rst | 2 +- 2 files changed, 10 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index 1fe048bc..ece5e6c3 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -160,18 +160,20 @@ persistence diagram with a family of field coefficients. \li Persistent_cohomology/rips_distance_matrix_persistence.cpp computes the Rips complex of a distance matrix and -outputs its persistence diagram. The file should contain lower diagonal distance matrix with semicolons as separators. +outputs its persistence diagram. + +The file should contain square or lower triangular distance matrix with semicolons as separators. The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. \li Persistent_cohomology/rips_correlation_matrix_persistence.cpp -computes the Rips complex of a correlation matrix and -outputs its persistence diagram. Note that no check is performed if -the matrix given as the input is a correlation matrix. -It is the user responsibility to ensure that this is the case. The -input is to be given either as a lower triangular matrix. -Please refer to data/distance_matrix/lower_triangular_correlation_matrix.csv for an example of a file. +computes the Rips complex of a correlation matrix and outputs its persistence diagram. + +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. The input is to be given either as a square or a lower +triangular matrix. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. \li Persistent_cohomology/alpha_complex_3d_persistence.cpp computes the persistent homology with diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst index f0e7bf2d..b80ff7fe 100644 --- a/src/cython/doc/rips_complex_user.rst +++ b/src/cython/doc/rips_complex_user.rst @@ -197,7 +197,7 @@ Example from csv file ^^^^^^^^^^^^^^^^^^^^^ This example builds the :doc:`Rips_complex ` from the given -points in an OFF file, and max_edge_length value. +distance matrix in a csv file, and max_edge_length value. Then it creates a :doc:`Simplex_tree ` with it. Finally, it is asked to display information about the Rips complex. -- cgit v1.2.3 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 (limited to 'src') 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 From 4a5332dfb88ec27157c82f77df87f588f7016736 Mon Sep 17 00:00:00 2001 From: pdlotko Date: Thu, 9 Nov 2017 14:48:39 +0000 Subject: Answer to Vincent's comments. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@2859 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 1a049ce98ccf893fc6af6c2631e508d3a9dc6c09 --- .../rips_correlation_matrix_persistence.cpp | 63 +++++++++------------- 1 file changed, 26 insertions(+), 37 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp index 94d5b8d4..676ef793 100644 --- a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp @@ -42,17 +42,17 @@ 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); + Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence); int main(int argc, char* argv[]) { std::string csv_matrix_file; std::string filediag; - Filtration_value threshold; + Filtration_value correlation_min; int dim_max; int p; Filtration_value min_persistence; - program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); + program_options(argc, argv, csv_matrix_file, filediag, correlation_min, dim_max, p, min_persistence); Correlation_matrix correlations = Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); @@ -68,13 +68,19 @@ int main(int argc, char* argv[]) { } } - // If the threshold, being minimal correlation is in the range [0,1], - // change it to 1-threshold - if ( ( threshold>=0 ) && ( threshold<=1 ) ) + Filtration_value threshold; + //If the correlation_min, being minimal corelation is in the range [0,1], + //change it to 1-correlation_min + if ( ( correlation_min>=0 ) && ( correlation_min<=1 ) ) { - threshold = 1-threshold; + threshold = 1-correlation_min; } - + else + { + std::cout << "Wrong value of the treshold corelation (should be between 0 and 1). The program will now terminate.\n"; + return 1; + } + Rips_complex rips_complex_from_file(correlations, threshold); // Construct the Rips complex in a Simplex Tree @@ -93,9 +99,9 @@ int main(int argc, char* argv[]) { pcoh.init_coefficients(p); //compute persistence pcoh.compute_persistent_cohomology(min_persistence); - - - // invert the persistence diagram + + + //invert the persistence diagram auto pairs = pcoh.get_persistent_pairs(); std::vector< intervals_common > processed_persistence_intervals; processed_persistence_intervals.reserve( pairs.size() ); @@ -108,12 +114,12 @@ int main(int argc, char* argv[]) { 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 + //and write them to a file if (filediag.empty()) { write_persistence_intervals_to_stream(processed_persistence_intervals); } else { @@ -125,37 +131,20 @@ int main(int argc, char* argv[]) { } 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) { + Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); - // hidden.add_options()( - // "input-file", po::value(&csv_matrix_file), - // "Name of file containing a correlation 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 correlation matrix. Can be square or lower triangular matrix. Separator is ';'."); - + "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")( - // "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")( - // "max-edge-length,r", - // po::value(&threshold)->default_value(std::numeric_limits::infinity()), - // "Maximal length 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"); 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-correlation,c", - po::value(&threshold)->default_value(std::numeric_limits::infinity()), - "Minimal correlation of an edge for the Rips complex construction.") + ("min-edge-corelation,c", + po::value(&correlation_min)->default_value(0), + "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), @@ -176,7 +165,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 correlation 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"; -- cgit v1.2.3 From baa98e772b9d97d8b7c22fc8292f54b5047f0569 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 2 Feb 2018 16:49:52 +0000 Subject: Fix readme and doc issue after merge git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3213 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 4fd050a6f686044aa6f6dcc290025e321e0dba1b --- .../doc/Intro_persistent_cohomology.h | 13 ++ .../rips_correlation_matrix_persistence.cpp | 179 --------------------- src/Rips_complex/utilities/CMakeLists.txt | 8 + .../rips_correlation_matrix_persistence.cpp | 179 +++++++++++++++++++++ src/Rips_complex/utilities/ripscomplex.md | 23 ++- 5 files changed, 222 insertions(+), 180 deletions(-) delete mode 100644 src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp create mode 100644 src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp (limited to 'src') diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index 4dbe82c7..3113a22c 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -162,6 +162,19 @@ persistence diagram with a family of field coefficients. Rips_complex/rips_distance_matrix_persistence.cpp computes the Rips complex of a distance matrix and outputs its persistence diagram. +The file should contain square or lower triangular distance matrix with semicolons as separators. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. + +\li +Rips_complex/rips_correlation_matrix_persistence.cpp +computes the Rips complex of a correlation matrix and outputs its persistence diagram. + +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. The input is to be given either as a square or a lower +triangular matrix. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. + \li Alpha_complex/alpha_complex_3d_persistence.cpp computes the persistent homology with \f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. diff --git a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp b/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp deleted file mode 100644 index 676ef793..00000000 --- a/src/Persistent_cohomology/example/rips_correlation_matrix_persistence.cpp +++ /dev/null @@ -1,179 +0,0 @@ -/* 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, Vincent Rouvreau - * - * Copyright (C) 2016 INRIA - * - * 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 . - */ - -#include -#include -#include -#include -#include - -#include - -#include -#include -#include // infinity - -// Types definition -using Simplex_tree = Gudhi::Simplex_tree; -using Filtration_value = Simplex_tree::Filtration_value; -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& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence); - -int main(int argc, char* argv[]) { - std::string csv_matrix_file; - std::string filediag; - Filtration_value correlation_min; - int dim_max; - int p; - Filtration_value min_persistence; - - program_options(argc, argv, csv_matrix_file, filediag, correlation_min, dim_max, p, min_persistence); - - Correlation_matrix correlations = - Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); - - // Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: - for (size_t i = 0; i != correlations.size(); ++i) { - for (size_t j = 0; j != correlations[i].size(); ++j) { - correlations[i][j] = 1 - correlations[i][j]; - if (correlations[i][j] < 0) { - std::cerr << "The input matrix is not a correlation matrix. \n"; - throw "The input matrix is not a correlation matrix. \n"; - } - } - } - - Filtration_value threshold; - //If the correlation_min, being minimal corelation is in the range [0,1], - //change it to 1-correlation_min - if ( ( correlation_min>=0 ) && ( correlation_min<=1 ) ) - { - threshold = 1-correlation_min; - } - else - { - std::cout << "Wrong value of the treshold corelation (should be between 0 and 1). The program will now terminate.\n"; - return 1; - } - - Rips_complex rips_complex_from_file(correlations, threshold); - - // Construct the Rips complex in a Simplex Tree - Simplex_tree simplex_tree; - - rips_complex_from_file.create_complex(simplex_tree, dim_max); - std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; - std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - simplex_tree.initialize_filtration(); - - // Compute the persistence diagram of the complex - Persistent_cohomology pcoh(simplex_tree); - // initializes the coefficient field for homology - pcoh.init_coefficients(p); - //compute persistence - pcoh.compute_persistent_cohomology(min_persistence); - - - //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()) { - write_persistence_intervals_to_stream(processed_persistence_intervals); - } else { - std::ofstream out(filediag); - write_persistence_intervals_to_stream(processed_persistence_intervals,out); - out.close(); - } - return 0; -} - -void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, - Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence) { - namespace po = boost::program_options; - po::options_description hidden("Hidden options"); - 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") - ("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(&correlation_min)->default_value(0), - "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); - - po::options_description all; - all.add(visible).add(hidden); - - po::variables_map vm; - po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); - po::notify(vm); - - 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 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"; - std::cout << "b and d are respectively the birth and death of the feature and \n"; - std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; - - std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; - std::cout << visible << std::endl; - std::abort(); - } -} diff --git a/src/Rips_complex/utilities/CMakeLists.txt b/src/Rips_complex/utilities/CMakeLists.txt index baa571fa..b99fc808 100644 --- a/src/Rips_complex/utilities/CMakeLists.txt +++ b/src/Rips_complex/utilities/CMakeLists.txt @@ -7,9 +7,13 @@ target_link_libraries(rips_distance_matrix_persistence ${Boost_PROGRAM_OPTIONS_L add_executable(rips_persistence rips_persistence.cpp) target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) +add_executable(rips_correlation_matrix_persistence rips_correlation_matrix_persistence.cpp) +target_link_libraries(rips_correlation_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + if (TBB_FOUND) target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES}) target_link_libraries(rips_persistence ${TBB_LIBRARIES}) + target_link_libraries(rips_correlation_matrix_persistence ${TBB_LIBRARIES}) endif() add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $ @@ -17,5 +21,9 @@ add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") +add_test(NAME Rips_complex_utility_from_rips_correlation_matrix COMMAND $ + "${CMAKE_SOURCE_DIR}/data/correlation_matrix/lower_triangular_correlation_matrix.csv" "-c" "0.3" "-d" "3" "-p" "3" "-m" "0") + install(TARGETS rips_distance_matrix_persistence DESTINATION bin) install(TARGETS rips_persistence DESTINATION bin) +install(TARGETS rips_correlation_matrix_persistence DESTINATION bin) diff --git a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp new file mode 100644 index 00000000..676ef793 --- /dev/null +++ b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp @@ -0,0 +1,179 @@ +/* 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, Vincent Rouvreau + * + * Copyright (C) 2016 INRIA + * + * 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 . + */ + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +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& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string csv_matrix_file; + std::string filediag; + Filtration_value correlation_min; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, correlation_min, dim_max, p, min_persistence); + + Correlation_matrix correlations = + Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); + + // Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + correlations[i][j] = 1 - correlations[i][j]; + if (correlations[i][j] < 0) { + std::cerr << "The input matrix is not a correlation matrix. \n"; + throw "The input matrix is not a correlation matrix. \n"; + } + } + } + + Filtration_value threshold; + //If the correlation_min, being minimal corelation is in the range [0,1], + //change it to 1-correlation_min + if ( ( correlation_min>=0 ) && ( correlation_min<=1 ) ) + { + threshold = 1-correlation_min; + } + else + { + std::cout << "Wrong value of the treshold corelation (should be between 0 and 1). The program will now terminate.\n"; + return 1; + } + + Rips_complex rips_complex_from_file(correlations, threshold); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + rips_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + //compute persistence + pcoh.compute_persistent_cohomology(min_persistence); + + + //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()) { + write_persistence_intervals_to_stream(processed_persistence_intervals); + } else { + std::ofstream out(filediag); + write_persistence_intervals_to_stream(processed_persistence_intervals,out); + out.close(); + } + return 0; +} + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + 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") + ("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(&correlation_min)->default_value(0), + "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); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); + po::notify(vm); + + 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 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"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md index 4291fae7..3f064e67 100644 --- a/src/Rips_complex/utilities/ripscomplex.md +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -39,11 +39,32 @@ Same as `rips_persistence` but taking a distance matrix as input. **Usage** -`rips_persistence [options] ` +`rips_distance_matrix_persistence [options] ` where `` is the path to the file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. **Example** `rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` + + +## rips_correlation_matrix_persistence ## + +Same as `rips_distance_matrix_persistence` but taking a correlation matrix as input. + +**Usage** + +`rips_correlation_matrix_persistence [options] ` + +where +`` is the path to the file containing a correlation matrix. Can be square or lower triangular matrix. Separator is ';'. +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. + +**Example** + +`rips_correlation_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` -- cgit v1.2.3 From 6d759732282c4b79a4c93a0dbceb4667af2c8940 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Mon, 5 Feb 2018 08:45:01 +0000 Subject: clang format files git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3216 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 09a40b92cc8e0cc1a08f962de14daed975338d1c --- src/Rips_complex/doc/Intro_rips_complex.h | 3 +- ...e_one_skeleton_rips_from_correlation_matrix.cpp | 11 ++- .../rips_correlation_matrix_persistence.cpp | 83 ++++++++++------------ 3 files changed, 45 insertions(+), 52 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 0c1acba7..3f02206b 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -167,7 +167,8 @@ namespace rips_complex { * * \include Rips_complex/one_skeleton_rips_from_correlation_matrix_for_doc.txt * - * All the other constructions discussed for Rips complex for distance matrix can be also performed for Rips complexes construction from correlation matrices. + * All the other constructions discussed for Rips complex for distance matrix can be also performed for Rips complexes + * construction from correlation matrices. * */ /** @} */ // end defgroup rips_complex 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 d1ccbf31..0da1dc20 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,16 @@ int main() { } //----------------------------------------------------------------------------- - // Now the correlation matrix is a 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. + // 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; @@ -63,7 +62,7 @@ int main() { // 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). + // 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/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp index 676ef793..3fe0edb0 100644 --- a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp +++ b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp @@ -39,7 +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 >; +using intervals_common = Gudhi::Persistence_interval_common; void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence); @@ -69,16 +69,14 @@ int main(int argc, char* argv[]) { } Filtration_value threshold; - //If the correlation_min, being minimal corelation is in the range [0,1], - //change it to 1-correlation_min - if ( ( correlation_min>=0 ) && ( correlation_min<=1 ) ) - { - threshold = 1-correlation_min; - } - else - { - std::cout << "Wrong value of the treshold corelation (should be between 0 and 1). The program will now terminate.\n"; - return 1; + // If the correlation_min, being minimal corelation is in the range [0,1], + // change it to 1-correlation_min + if ((correlation_min >= 0) && (correlation_min <= 1)) { + threshold = 1 - correlation_min; + } else { + std::cout + << "Wrong value of the treshold corelation (should be between 0 and 1). The program will now terminate.\n"; + return 1; } Rips_complex rips_complex_from_file(correlations, threshold); @@ -97,34 +95,30 @@ int main(int argc, char* argv[]) { Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(p); - //compute persistence + // compute persistence pcoh.compute_persistent_cohomology(min_persistence); - - //invert the persistence diagram + // 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) - ); + std::vector 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() ); + // sort the processed intervals: + std::sort(processed_persistence_intervals.begin(), processed_persistence_intervals.end()); - //and write them to a file + // and write them to a file if (filediag.empty()) { write_persistence_intervals_to_stream(processed_persistence_intervals); } else { std::ofstream out(filediag); - write_persistence_intervals_to_stream(processed_persistence_intervals,out); + write_persistence_intervals_to_stream(processed_persistence_intervals, out); out.close(); } return 0; @@ -134,23 +128,22 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std:: Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); - 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 ';'."); + 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") - ("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(&correlation_min)->default_value(0), - "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"); + 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(&correlation_min)->default_value(0), + "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); -- cgit v1.2.3 From 3787c1b1473a590515a62bfb0b616ce17045f731 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Mon, 5 Feb 2018 10:02:52 +0000 Subject: Do like C++ example git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3217 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: cc963c77e30c649e7428afc609bac4ba1f515a37 --- ...rsistence_from_correlation_matrix_file_example.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py index 79f6df75..ad990fdc 100755 --- a/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py +++ b/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import gudhi +import sys import argparse """This file is part of the Gudhi Library. The Gudhi library @@ -37,23 +38,29 @@ parser = argparse.ArgumentParser(description='RipsComplex creation from ' '- Constructs a Rips complex with the ' 'correlation matrix from the given csv file.') parser.add_argument("-f", "--file", type=str, required=True) -parser.add_argument("-e", "--max_edge_length", type=float, default=0.5) +parser.add_argument("-c", "--min_edge_correlation", type=float, default=0.5) parser.add_argument("-d", "--max_dimension", type=int, default=1) parser.add_argument("-b", "--band_boot", type=float, default=0.) parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams') args = parser.parse_args() +if not (-1. < args.min_edge_correlation < 1.): + print("Wrong value of the treshold corelation (should be between -1 and 1).") + sys.exit(1) + print("#####################################################################") print("RipsComplex creation from correlation matrix read in a csv file") -message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) +message = "RipsComplex with min_edge_correlation=" + repr(args.min_edge_correlation) print(message) correlation_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=args.file) # Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: -distance_matrix = [[1-correlation_matrix[i][j] for j in range(len(correlation_matrix[0]))] for i in range(len(correlation_matrix))] -rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=args.max_edge_length) +distance_matrix = [[1.-correlation_matrix[i][j] for j in range(len(correlation_matrix[i]))] for i in range(len(correlation_matrix))] + +rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, + max_edge_length=1.-args.min_edge_correlation) simplex_tree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) message = "Number of simplices=" + repr(simplex_tree.num_simplices()) @@ -64,6 +71,9 @@ diag = simplex_tree.persistence() print("betti_numbers()=") print(simplex_tree.betti_numbers()) +# invert the persistence diagram +invert_diag = [(diag[pers][0],(1.-diag[pers][1][0], 1.-diag[pers][1][1])) for pers in range(len(diag))] + if args.no_diagram == False: - pplot = gudhi.plot_persistence_diagram(diag, band_boot=args.band_boot) + pplot = gudhi.plot_persistence_diagram(invert_diag, band_boot=args.band_boot) pplot.show() -- cgit v1.2.3 From fe97343e5b53fe9d93e6aab3312021c53f2f0923 Mon Sep 17 00:00:00 2001 From: pdlotko Date: Thu, 8 Feb 2018 09:57:19 +0000 Subject: Commit to remove one file. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3231 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5795ef0322c1ba328e8bee068009a27f18b7224c --- src/common/include/gudhi/file_writer.h | 19 +-- .../include/gudhi/writing_persistence_to_file.h | 158 +++++++++++++++++++++ 2 files changed, 168 insertions(+), 9 deletions(-) create mode 100644 src/common/include/gudhi/writing_persistence_to_file.h (limited to 'src') diff --git a/src/common/include/gudhi/file_writer.h b/src/common/include/gudhi/file_writer.h index 1b59ae46..53c83533 100644 --- a/src/common/include/gudhi/file_writer.h +++ b/src/common/include/gudhi/file_writer.h @@ -41,24 +41,24 @@ 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() ){} + 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()){} + 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){} + 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_) + (this->dimension_ == i2.dimension_) && (this->arith_element_ == i2.arith_element_) ); } @@ -116,9 +116,9 @@ public: friend std::ostream& operator<<(std::ostream& out, const Persistence_interval_common& it) { - if ( it.Arith_element_ != std::numeric_limits::max() ) + if ( it.arith_element_ != std::numeric_limits::max() ) { - out << it.Arith_element_ << " "; + out << it.arith_element_ << " "; } if ( it.dimension_ != std::numeric_limits::max() ) { @@ -132,7 +132,7 @@ private: Filtration_type birth_; Filtration_type death_; unsigned dimension_; - Coefficient_field Arith_element_; + Coefficient_field arith_element_; };//Persistence_interval_common @@ -141,12 +141,13 @@ private: **/ template void write_persistence_intervals_to_stream( -const std::vector< Persistence_interval_common >& intervals, +//const std::vector< Persistence_interval_common >& intervals, + std::ostream& out = std::cout ) { for ( auto interval : intervals ) { - out << interval << std::endl; + out << interval << "\n"; } }//write_persistence_intervals_to_stream diff --git a/src/common/include/gudhi/writing_persistence_to_file.h b/src/common/include/gudhi/writing_persistence_to_file.h new file mode 100644 index 00000000..53c83533 --- /dev/null +++ b/src/common/include/gudhi/writing_persistence_to_file.h @@ -0,0 +1,158 @@ +/* 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 << "\n"; + } +}//write_persistence_intervals_to_stream + + + +} + +#endif //FILE_WRITER_ -- cgit v1.2.3 From 82764fe1fa117f5a3d1cce0eb7dc7dddd7af974f Mon Sep 17 00:00:00 2001 From: pdlotko Date: Thu, 8 Feb 2018 09:57:39 +0000 Subject: The file has been removed. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3232 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 4e6e6b2c1939506e19bb76a7bb507f869160577b --- src/common/include/gudhi/file_writer.h | 158 --------------------------------- 1 file changed, 158 deletions(-) delete mode 100644 src/common/include/gudhi/file_writer.h (limited to 'src') diff --git a/src/common/include/gudhi/file_writer.h b/src/common/include/gudhi/file_writer.h deleted file mode 100644 index 53c83533..00000000 --- a/src/common/include/gudhi/file_writer.h +++ /dev/null @@ -1,158 +0,0 @@ -/* 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 << "\n"; - } -}//write_persistence_intervals_to_stream - - - -} - -#endif //FILE_WRITER_ -- cgit v1.2.3 From b5667b282e8cee9e1d2bfcf90bdfbbf2bc0030b2 Mon Sep 17 00:00:00 2001 From: pdlotko Date: Thu, 8 Feb 2018 10:00:57 +0000 Subject: random commit. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3233 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 9a3ccdd8a000bd3fd79bab61b35e50bd707cfb29 --- .../include/gudhi/Bitmap_cubical_complex.h | 2 +- ...e_one_skeleton_rips_from_correlation_matrix.cpp | 23 +++++++++++----------- .../include/gudhi/writing_persistence_to_file.h | 10 +++++----- 3 files changed, 18 insertions(+), 17 deletions(-) (limited to 'src') diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index 969daba6..770eb55f 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -383,7 +383,7 @@ class Bitmap_cubical_complex : public T { std::vector bdry = this->get_boundary_of_a_cell(sh); if (globalDbg) { std::cerr << "std::pair endpoints( Simplex_handle sh )\n"; - std::cerr << "bdry.size() : " << bdry.size() << std::endl; + std::cerr << "bdry.size() : " << bdry.size() << "\n"; } // this method returns two first elements from the boundary of sh. if (bdry.size() < 2) 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 0da1dc20..f66c4b04 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 @@ -31,26 +31,27 @@ int main() { // ---------------------------------------------------------------------------- // Convert correlation matrix to a distance matrix: // ---------------------------------------------------------------------------- + double threshold = 0; for (size_t i = 0; i != correlations.size(); ++i) { for (size_t j = 0; j != correlations[i].size(); ++j) { - correlations[i][j] = 1 - correlations[i][j]; - if (correlations[i][j] < 0) { - std::cerr << "The input matrix is not a correlation matrix. \n"; - throw "The input matrix is not a correlation matrix. \n"; - } + //Here we check if our data comes from corelation matrix. + if ( (correlations[i][j]<-1) || (correlations[i][j]>1) ) + { + std::cerr << "The input matrix is not a correlation matrix. The program will now terminate.\n"; + throw "The input matrix is not a correlation matrix. The program will now terminate.\n"; + } + correlations[i][j] = 1 - correlations[i][j]; } + //Here we make sure that we will get the treshold value equal to maximal + //distance in the matrix. + if ( correlations[i][j] > threshold )threshold = correlations[i][j]; } //----------------------------------------------------------------------------- // 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); diff --git a/src/common/include/gudhi/writing_persistence_to_file.h b/src/common/include/gudhi/writing_persistence_to_file.h index 53c83533..5767c06e 100644 --- a/src/common/include/gudhi/writing_persistence_to_file.h +++ b/src/common/include/gudhi/writing_persistence_to_file.h @@ -20,8 +20,8 @@ * along with this program. If not, see . */ -#ifndef FILE_WRITER_ -#define FILE_WRITER_ +#ifndef WRITING_PERSISTENCE_TO_FILE_H +#define WRITING_PERSISTENCE_TO_FILE_H #include #include @@ -141,8 +141,8 @@ private: **/ template void write_persistence_intervals_to_stream( -//const std::vector< Persistence_interval_common >& intervals, - +const std::vector< Persistence_interval_common >& intervals, +//TODO: change to ranges when it is clear how to do that. std::ostream& out = std::cout ) { for ( auto interval : intervals ) @@ -155,4 +155,4 @@ void write_persistence_intervals_to_stream( } -#endif //FILE_WRITER_ +#endif //WRITING_PERSISTENCE_TO_FILE_H -- cgit v1.2.3 From 82c0652b06949b0b002781688565d7ecf30f04fe Mon Sep 17 00:00:00 2001 From: pdlotko Date: Thu, 8 Feb 2018 10:37:21 +0000 Subject: Fixed comments from Marc and Vincent. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3234 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 8b5e08c38cd05512ba8b2824f13a62d35f39bac3 --- ...e_one_skeleton_rips_from_correlation_matrix.cpp | 10 +- .../rips_correlation_matrix_persistence.cpp | 40 ++++---- .../include/gudhi/writing_persistence_to_file.h | 102 ++++++++++++--------- 3 files changed, 86 insertions(+), 66 deletions(-) (limited to 'src') 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 f66c4b04..a34ce15c 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 @@ -40,11 +40,11 @@ int main() { std::cerr << "The input matrix is not a correlation matrix. The program will now terminate.\n"; throw "The input matrix is not a correlation matrix. The program will now terminate.\n"; } - correlations[i][j] = 1 - correlations[i][j]; - } - //Here we make sure that we will get the treshold value equal to maximal - //distance in the matrix. - if ( correlations[i][j] > threshold )threshold = correlations[i][j]; + correlations[i][j] = 1 - correlations[i][j]; + //Here we make sure that we will get the treshold value equal to maximal + //distance in the matrix. + if ( correlations[i][j] > threshold )threshold = correlations[i][j]; + } } //----------------------------------------------------------------------------- diff --git a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp index 3fe0edb0..95bce491 100644 --- a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp +++ b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp @@ -24,7 +24,7 @@ #include #include #include -#include +#include #include @@ -57,28 +57,23 @@ int main(int argc, char* argv[]) { Correlation_matrix correlations = Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); + Filtration_value threshold = 0; + // Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: for (size_t i = 0; i != correlations.size(); ++i) { for (size_t j = 0; j != correlations[i].size(); ++j) { correlations[i][j] = 1 - correlations[i][j]; - if (correlations[i][j] < 0) { - std::cerr << "The input matrix is not a correlation matrix. \n"; - throw "The input matrix is not a correlation matrix. \n"; + //Here we make sure that the values of corelations lie between -1 and 1. + //If not, we throw an exception. + if ((correlations[i][j] < -1) || (correlations[i][j] > 1)) + { + std::cerr << "The input matrix is not a correlation matrix. The program will now terminate. \n"; + throw "The input matrix is not a correlation matrix. The program will now terminate. \n"; } + if ( correlations[i][j] > threshold ) threshold = correlations[i][j]; } } - Filtration_value threshold; - // If the correlation_min, being minimal corelation is in the range [0,1], - // change it to 1-correlation_min - if ((correlation_min >= 0) && (correlation_min <= 1)) { - threshold = 1 - correlation_min; - } else { - std::cout - << "Wrong value of the treshold corelation (should be between 0 and 1). The program will now terminate.\n"; - return 1; - } - Rips_complex rips_complex_from_file(correlations, threshold); // Construct the Rips complex in a Simplex Tree @@ -95,10 +90,16 @@ int main(int argc, char* argv[]) { Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(p); - // compute persistence + //compute persistence pcoh.compute_persistent_cohomology(min_persistence); - - // invert the persistence diagram + + + //invert the persistence diagram. The reason for this procedure is the following: + //The input to the program is a corelation matrix M. When processing it, it is + //turned into 1-M and the obtained persistence intervals are in '1-M' units. + //Below we reverse every (birth,death) pair into (1-birth, 1-death) pair + //so that the input and the output to the program is expressed in the same + //units. auto pairs = pcoh.get_persistent_pairs(); std::vector processed_persistence_intervals; processed_persistence_intervals.reserve(pairs.size()); @@ -118,8 +119,7 @@ int main(int argc, char* argv[]) { write_persistence_intervals_to_stream(processed_persistence_intervals); } else { std::ofstream out(filediag); - write_persistence_intervals_to_stream(processed_persistence_intervals, out); - out.close(); + write_persistence_intervals_to_stream(processed_persistence_intervals, out); } return 0; } diff --git a/src/common/include/gudhi/writing_persistence_to_file.h b/src/common/include/gudhi/writing_persistence_to_file.h index 5767c06e..5457cf48 100644 --- a/src/common/include/gudhi/writing_persistence_to_file.h +++ b/src/common/include/gudhi/writing_persistence_to_file.h @@ -39,21 +39,35 @@ template class Persistence_interval_common { public: + /** + * Constructor taking as an input birth and death of the pair. + **/ Persistence_interval_common( Filtration_type birth , Filtration_type death ): birth_(birth),death_(death),dimension_(std::numeric_limits::max), arith_element_(std::numeric_limits::max() ){} + /** + * Constructor taking as an input birth, death and dimension of the pair. + **/ Persistence_interval_common( Filtration_type birth , Filtration_type death, unsigned dim ): birth_(birth),death_(death),dimension_(dim), arith_element_(std::numeric_limits::max()){} + /** + * Constructor taking as an input birth, death, dimension of the pair as well + * as the number p such that this interval is present over Z_p field. + **/ Persistence_interval_common( Filtration_type birth , Filtration_type death, unsigned dim , Coefficient_field field ): birth_(birth),death_(death),dimension_(dim), arith_element_(field){} - + /** + * Operator to compare two persistence pairs. During the comparision all the + * fields: birth, death, dimensiona and arith_element_ are taken into account + * and they all have to be equal for two pairs to be equal. + **/ inline bool operator == ( const Persistence_interval_common &i2) { return ( @@ -62,6 +76,9 @@ public: ); } + /** + * Check if two persistence paris are not equal. + **/ inline bool operator != ( const Persistence_interval_common &i2) { return (!((*this)==i2)); @@ -69,49 +86,52 @@ public: /** + * Operator to compare objects of a type Persistence_interval_common. + * One intervals is smaller than the other if it has lower persistence. * 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; - } - } - } - } - } + { + return fabs( this->death_-this->birth_ ) < fabs( i2.death_-i2.birth_ ); + //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) -- cgit v1.2.3 From a6b1719067462f93a7ee3b7f8f632b85fe27117d Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 3 Apr 2018 08:40:05 +0000 Subject: Add warnings for rips persistence from correlation matrix (points under the diagonal). Use range instead of vector in write_persistence_intervals_to_stream git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/rips_complex_from_correlation_matrix@3326 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: ea0fdc32c91f1927c5a9f2a3864360895ce8128a --- src/Rips_complex/doc/Intro_rips_complex.h | 3 + src/Rips_complex/utilities/ripscomplex.md | 5 + .../include/gudhi/writing_persistence_to_file.h | 212 ++++++++------------- src/cython/doc/rips_complex_user.rst | 5 + ...istence_from_correlation_matrix_file_example.py | 5 + 5 files changed, 93 insertions(+), 137 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 3f02206b..496c4218 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -170,6 +170,9 @@ namespace rips_complex { * All the other constructions discussed for Rips complex for distance matrix can be also performed for Rips complexes * construction from correlation matrices. * + * @warning As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical + * tool will not work properly, this is a known issue. + * */ /** @} */ // end defgroup rips_complex diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md index 3f064e67..857e6293 100644 --- a/src/Rips_complex/utilities/ripscomplex.md +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -68,3 +68,8 @@ Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv **Example** `rips_correlation_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` + +**Warning** + +As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical tool will not work +properly, this is a known issue. diff --git a/src/common/include/gudhi/writing_persistence_to_file.h b/src/common/include/gudhi/writing_persistence_to_file.h index 5457cf48..5020b5fb 100644 --- a/src/common/include/gudhi/writing_persistence_to_file.h +++ b/src/common/include/gudhi/writing_persistence_to_file.h @@ -29,150 +29,88 @@ namespace Gudhi { - /** -* This is a class to store persistence intervals. Its main purpose is to +* 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. +* of writing a collection of persistence intervals to file. **/ -template -class Persistence_interval_common -{ -public: - /** - * Constructor taking as an input birth and death of the pair. - **/ - Persistence_interval_common( Filtration_type birth , Filtration_type death ): - birth_(birth),death_(death),dimension_(std::numeric_limits::max), - arith_element_(std::numeric_limits::max() ){} - - /** - * Constructor taking as an input birth, death and dimension of the pair. - **/ - Persistence_interval_common( Filtration_type birth , Filtration_type death, - unsigned dim ): - birth_(birth),death_(death),dimension_(dim), - arith_element_(std::numeric_limits::max()){} - - /** - * Constructor taking as an input birth, death, dimension of the pair as well - * as the number p such that this interval is present over Z_p field. - **/ - Persistence_interval_common( Filtration_type birth , Filtration_type death, - unsigned dim , Coefficient_field field ): - birth_(birth),death_(death),dimension_(dim), - arith_element_(field){} - - /** - * Operator to compare two persistence pairs. During the comparision all the - * fields: birth, death, dimensiona and arith_element_ are taken into account - * and they all have to be equal for two pairs to be equal. - **/ - 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_) - ); - } - - /** - * Check if two persistence paris are not equal. - **/ - inline bool operator != ( const Persistence_interval_common &i2) - { - return (!((*this)==i2)); - } - - - /** - * Operator to compare objects of a type Persistence_interval_common. - * One intervals is smaller than the other if it has lower persistence. - * Note that this operator do not take Arith_element into account when doing comparisions. - **/ - inline bool operator < ( const Persistence_interval_common &i2) - { - return fabs( this->death_-this->birth_ ) < fabs( i2.death_-i2.birth_ ); - //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 +template +class Persistence_interval_common { + public: + /** + * Constructor taking as an input birth and death of the pair. + **/ + Persistence_interval_common(Filtration_type birth, Filtration_type death) + : birth_(birth), + death_(death), + dimension_(std::numeric_limits::max), + arith_element_(std::numeric_limits::max()) {} + /** + * Constructor taking as an input birth, death and dimension of the pair. + **/ + Persistence_interval_common(Filtration_type birth, Filtration_type death, unsigned dim) + : birth_(birth), death_(death), dimension_(dim), arith_element_(std::numeric_limits::max()) {} -/** - * This function write a vector to a stream -**/ -template -void write_persistence_intervals_to_stream( -const std::vector< Persistence_interval_common >& intervals, -//TODO: change to ranges when it is clear how to do that. - std::ostream& out = std::cout ) -{ - for ( auto interval : intervals ) - { - out << interval << "\n"; - } -}//write_persistence_intervals_to_stream + /** +* Constructor taking as an input birth, death, dimension of the pair as well +* as the number p such that this interval is present over Z_p field. +**/ + Persistence_interval_common(Filtration_type birth, Filtration_type death, unsigned dim, Coefficient_field field) + : birth_(birth), death_(death), dimension_(dim), arith_element_(field) {} + + /** + * Operator to compare two persistence pairs. During the comparision all the + * fields: birth, death, dimensiona and arith_element_ are taken into account + * and they all have to be equal for two pairs to be equal. + **/ + 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_)); + } + + /** + * Check if two persistence paris are not equal. + **/ + inline bool operator!=(const Persistence_interval_common& i2) { return (!((*this) == i2)); } + /** + * Operator to compare objects of a type Persistence_interval_common. + * One intervals is smaller than the other if it has lower persistence. + * Note that this operator do not take Arith_element into account when doing comparisions. + **/ + inline bool operator<(const Persistence_interval_common& i2) { + return fabs(this->death_ - this->birth_) < fabs(i2.death_ - i2.birth_); + } + 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_; +}; + +/** + * This function write a vector to a stream +**/ +template +void write_persistence_intervals_to_stream(const Persistence_interval_range& intervals, + std::ostream& out = std::cout) { + for (auto interval : intervals) { + out << interval << "\n"; + } +} } -#endif //WRITING_PERSISTENCE_TO_FILE_H +#endif // WRITING_PERSISTENCE_TO_FILE_H diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst index b80ff7fe..7738aef0 100644 --- a/src/cython/doc/rips_complex_user.rst +++ b/src/cython/doc/rips_complex_user.rst @@ -305,3 +305,8 @@ until dimension 1 - one skeleton graph in other words), the output is: [2, 4] -> 0.97 [0, 3] -> 0.99 [1, 3] -> 0.99 + +.. note:: + As persistence diagrams points will be under the diagonal, + bottleneck distance and persistence graphical tool will not work properly, + this is a known issue. diff --git a/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py index ad990fdc..aa82ef71 100755 --- a/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py +++ b/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py @@ -49,6 +49,11 @@ if not (-1. < args.min_edge_correlation < 1.): print("Wrong value of the treshold corelation (should be between -1 and 1).") sys.exit(1) +print("#####################################################################") +print("Caution: as persistence diagrams points will be under the diagonal,") +print("bottleneck distance and persistence graphical tool will not work") +print("properly, this is a known issue.") + print("#####################################################################") print("RipsComplex creation from correlation matrix read in a csv file") -- cgit v1.2.3