diff options
Diffstat (limited to 'cython/cython')
| -rwxr-xr-x | cython/cython/persistence_graphical_tools.py | 76 | ||||
| -rw-r--r-- | cython/cython/reader_utils.pyx | 95 | ||||
| -rw-r--r-- | cython/cython/simplex_tree.pyx | 81 |
3 files changed, 204 insertions, 48 deletions
diff --git a/cython/cython/persistence_graphical_tools.py b/cython/cython/persistence_graphical_tools.py index a984633e..fb837e29 100755 --- a/cython/cython/persistence_graphical_tools.py +++ b/cython/cython/persistence_graphical_tools.py @@ -1,11 +1,12 @@ import matplotlib.pyplot as plt import numpy as np +import os """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): Vincent Rouvreau + Author(s): Vincent Rouvreau, Bertrand Michel Copyright (C) 2016 INRIA @@ -23,15 +24,17 @@ import numpy as np along with this program. If not, see <http://www.gnu.org/licenses/>. """ -__author__ = "Vincent Rouvreau" +__author__ = "Vincent Rouvreau, Bertrand Michel" __copyright__ = "Copyright (C) 2016 INRIA" __license__ = "GPL v3" -def __min_birth_max_death(persistence): +def __min_birth_max_death(persistence, band_boot=0.): """This function returns (min_birth, max_death) from the persistence. :param persistence: The persistence to plot. :type persistence: list of tuples(dimension, tuple(birth, death)). + :param band_boot: bootstrap band + :type band_boot: float. :returns: (float, float) -- (min_birth, max_death). """ # Look for minimum birth date and maximum death date for plot optimisation @@ -45,6 +48,8 @@ def __min_birth_max_death(persistence): max_death = float(interval[1][0]) if float(interval[1][0]) < min_birth: min_birth = float(interval[1][0]) + if band_boot > 0.: + max_death += band_boot return (min_birth, max_death) """ @@ -59,7 +64,7 @@ def show_palette_values(alpha=0.6): :param alpha: alpha value in [0.0, 1.0] for horizontal bars (default is 0.6). :type alpha: float. - :returns: plot -- An horizontal bar plot of dimensions color. + :returns: plot the dimension palette values. """ colors = [] for color in palette: @@ -70,18 +75,38 @@ def show_palette_values(alpha=0.6): plt.barh(y_pos, y_pos + 1, align='center', alpha=alpha, color=colors) plt.ylabel('Dimension') plt.title('Dimension palette values') + return plt - plt.show() - -def plot_persistence_barcode(persistence, alpha=0.6): +def plot_persistence_barcode(persistence=[], persistence_file='', alpha=0.6, max_barcodes=0): """This function plots the persistence bar code. :param persistence: The persistence to plot. :type persistence: list of tuples(dimension, tuple(birth, death)). + :param persistence_file: A persistence file style name (reset persistence if both are set). + :type persistence_file: string :param alpha: alpha value in [0.0, 1.0] for horizontal bars (default is 0.6). :type alpha: float. + :param max_barcodes: number of maximal barcodes to be displayed + (persistence will be sorted by life time if max_barcodes is set) + :type max_barcodes: int. :returns: plot -- An horizontal bar plot of persistence. """ + if persistence_file is not '': + if os.path.isfile(persistence_file): + # Reset persistence + persistence = [] + diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file) + for key in diag.keys(): + for persistence_interval in diag[key]: + persistence.append((key, persistence_interval)) + else: + print("file " + persistence_file + " not found.") + return None + + if max_barcodes > 0 and max_barcodes < len(persistence): + # Sort by life time, then takes only the max_plots elements + persistence = sorted(persistence, key=lambda life_time: life_time[1][1]-life_time[1][0], reverse=True)[:max_barcodes] + (min_birth, max_death) = __min_birth_max_death(persistence) ind = 0 delta = ((max_death - min_birth) / 10.0) @@ -106,18 +131,40 @@ def plot_persistence_barcode(persistence, alpha=0.6): plt.title('Persistence barcode') # Ends plot on infinity value and starts a little bit before min_birth plt.axis([axis_start, infinity, 0, ind]) - plt.show() + return plt -def plot_persistence_diagram(persistence, alpha=0.6): - """This function plots the persistence diagram. +def plot_persistence_diagram(persistence=[], persistence_file='', alpha=0.6, band_boot=0., max_plots=0): + """This function plots the persistence diagram with an optional confidence band. :param persistence: The persistence to plot. :type persistence: list of tuples(dimension, tuple(birth, death)). + :param persistence_file: A persistence file style name (reset persistence if both are set). + :type persistence_file: string :param alpha: alpha value in [0.0, 1.0] for points and horizontal infinity line (default is 0.6). :type alpha: float. - :returns: plot -- An diagram plot of persistence. + :param band_boot: bootstrap band (not displayed if :math:`\leq` 0.) + :type band_boot: float. + :param max_plots: number of maximal plots to be displayed + :type max_plots: int. + :returns: plot -- A diagram plot of persistence. """ - (min_birth, max_death) = __min_birth_max_death(persistence) + if persistence_file is not '': + if os.path.isfile(persistence_file): + # Reset persistence + persistence = [] + diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file) + for key in diag.keys(): + for persistence_interval in diag[key]: + persistence.append((key, persistence_interval)) + else: + print("file " + persistence_file + " not found.") + return None + + if max_plots > 0 and max_plots < len(persistence): + # Sort by life time, then takes only the max_plots elements + persistence = sorted(persistence, key=lambda life_time: life_time[1][1]-life_time[1][0], reverse=True)[:max_plots] + + (min_birth, max_death) = __min_birth_max_death(persistence, band_boot) ind = 0 delta = ((max_death - min_birth) / 10.0) # Replace infinity values with max_death + delta for diagram to be more @@ -131,6 +178,9 @@ def plot_persistence_diagram(persistence, alpha=0.6): plt.plot(x, x, color='k', linewidth=1.0) plt.plot(x, [infinity] * len(x), linewidth=1.0, color='k', alpha=alpha) plt.text(axis_start, infinity, r'$\infty$', color='k', alpha=alpha) + # bootstrap band + if band_boot > 0.: + plt.fill_between(x, x, x+band_boot, alpha=alpha, facecolor='red') # Draw points in loop for interval in reversed(persistence): @@ -149,4 +199,4 @@ def plot_persistence_diagram(persistence, alpha=0.6): plt.ylabel('Death') # Ends plot on infinity value and starts a little bit before min_birth plt.axis([axis_start, infinity, axis_start, infinity + delta]) - plt.show() + return plt diff --git a/cython/cython/reader_utils.pyx b/cython/cython/reader_utils.pyx new file mode 100644 index 00000000..3a17c5a0 --- /dev/null +++ b/cython/cython/reader_utils.pyx @@ -0,0 +1,95 @@ +from cython cimport numeric +from libcpp.vector cimport vector +from libcpp.string cimport string +from libcpp.map cimport map +from libcpp.pair cimport pair +import os + +"""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): Vincent Rouvreau + + Copyright (C) 2017 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 <http://www.gnu.org/licenses/>. +""" + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2017 INRIA" +__license__ = "GPL v3" + +cdef extern from "Reader_utils_interface.h" namespace "Gudhi": + vector[vector[double]] read_matrix_from_csv_file(string off_file, char separator) + map[int, vector[pair[double, double]]] read_pers_intervals_grouped_by_dimension(string filename) + vector[pair[double, double]] read_pers_intervals_in_dimension(string filename, int only_this_dim) + +def read_lower_triangular_matrix_from_csv_file(csv_file='', separator=';'): + """Read lower triangular matrix from a CSV style file. + + :param csv_file: A CSV file style name. + :type csv_file: string + :param separator: The value separator in the CSV file. Default value is ';' + :type separator: char + + :returns: The lower triangular matrix. + :rtype: vector[vector[double]] + """ + if csv_file is not '': + if os.path.isfile(csv_file): + return read_matrix_from_csv_file(str.encode(csv_file), ord(separator[0])) + print("file " + csv_file + " not set or not found.") + return [] + +def read_persistence_intervals_grouped_by_dimension(persistence_file=''): + """Reads a file containing persistence intervals. + Each line might contain 2, 3 or 4 values: [[field] dimension] birth death + The return value is an `map[dim, vector[pair[birth, death]]]` + where `dim` is an `int`, `birth` a `double`, and `death` a `double`. + Note: the function does not check that birth <= death. + + :param persistence_file: A persistence file style name. + :type persistence_file: string + + :returns: The persistence pairs grouped by dimension. + :rtype: map[int, vector[pair[double, double]]] + """ + if persistence_file is not '': + if os.path.isfile(persistence_file): + return read_pers_intervals_grouped_by_dimension(str.encode(persistence_file)) + print("file " + persistence_file + " not set or not found.") + return [] + +def read_persistence_intervals_in_dimension(persistence_file='', only_this_dim=-1): + """Reads a file containing persistence intervals. + Each line might contain 2, 3 or 4 values: [[field] dimension] birth death + If `only_this_dim` = -1, dimension is ignored and all lines are returned. + If `only_this_dim` is >= 0, only the lines where dimension = `only_this_dim` + (or where dimension is not specified) are returned. + The return value is an `vector[pair[birth, death]]` + where `birth` a `double`, and `death` a `double`. + Note: the function does not check that birth <= death. + + :param persistence_file: A persistence file style name. + :type persistence_file: string + + :returns: The persistence pairs grouped by dimension. + :rtype: map[int, vector[pair[double, double]]] + """ + if persistence_file is not '': + if os.path.isfile(persistence_file): + return read_pers_intervals_in_dimension(str.encode(persistence_file), only_this_dim) + print("file " + persistence_file + " not set or not found.") + return [] diff --git a/cython/cython/simplex_tree.pyx b/cython/cython/simplex_tree.pyx index 9d40a8b5..45487158 100644 --- a/cython/cython/simplex_tree.pyx +++ b/cython/cython/simplex_tree.pyx @@ -2,6 +2,7 @@ from cython cimport numeric from libcpp.vector cimport vector from libcpp.utility cimport pair from libcpp cimport bool +from libcpp.string cimport string """This file is part of the Gudhi Library. The Gudhi library (Geometric Understanding in Higher Dimensions) is a generic C++ @@ -35,9 +36,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>": Simplex_tree() - double filtration() double simplex_filtration(vector[int] simplex) - void set_filtration(double filtration) void initialize_filtration() int num_vertices() int num_simplices() @@ -61,6 +60,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": vector[int] betti_numbers() vector[int] persistent_betti_numbers(double from_value, double to_value) vector[pair[double,double]] intervals_in_dimension(int dimension) + void write_output_diagram(string diagram_file_name) # SimplexTree python interface cdef class SimplexTree: @@ -113,14 +113,6 @@ cdef class SimplexTree: """ return self.thisptr.simplex_filtration(simplex) - def set_filtration(self, filtration): - """This function sets the main simplicial complex filtration value. - - :param filtration: The filtration value. - :type filtration: float. - """ - self.thisptr.set_filtration(<double> filtration) - def initialize_filtration(self): """This function initializes and sorts the simplicial complex filtration vector. @@ -183,10 +175,10 @@ cdef class SimplexTree: :returns: true if the simplex was found, false otherwise. :rtype: bool """ - cdef vector[int] complex + cdef vector[int] csimplex for i in simplex: - complex.push_back(i) - return self.thisptr.find_simplex(complex) + csimplex.push_back(i) + return self.thisptr.find_simplex(csimplex) def insert(self, simplex, filtration=0.0): """This function inserts the given N-simplex and its subfaces with the @@ -200,10 +192,10 @@ cdef class SimplexTree: :returns: true if the simplex was found, false otherwise. :rtype: bool """ - cdef vector[int] complex + cdef vector[int] csimplex for i in simplex: - complex.push_back(i) - return self.thisptr.insert_simplex_and_subfaces(complex, + csimplex.push_back(i) + return self.thisptr.insert_simplex_and_subfaces(csimplex, <double>filtration) def get_filtration(self): @@ -232,35 +224,35 @@ cdef class SimplexTree: :returns: The (simplices of the) skeleton of a maximum dimension. :rtype: list of tuples(simplex, filtration) """ - cdef vector[pair[vector[int], double]] skeletons \ + cdef vector[pair[vector[int], double]] skeleton \ = self.thisptr.get_skeleton(<int>dimension) ct = [] - for filtered_complex in skeletons: + for filtered_simplex in skeleton: v = [] - for vertex in filtered_complex.first: + for vertex in filtered_simplex.first: v.append(vertex) - ct.append((v, filtered_complex.second)) + ct.append((v, filtered_simplex.second)) return ct def get_star(self, simplex): - """This function returns the stars of a given N-simplex. + """This function returns the star of a given N-simplex. :param simplex: The N-simplex, represented by a list of vertex. :type simplex: list of int. :returns: The (simplices of the) star of a simplex. :rtype: list of tuples(simplex, filtration) """ - cdef vector[int] complex + cdef vector[int] csimplex for i in simplex: - complex.push_back(i) - cdef vector[pair[vector[int], double]] stars \ - = self.thisptr.get_star(complex) + csimplex.push_back(i) + cdef vector[pair[vector[int], double]] star \ + = self.thisptr.get_star(csimplex) ct = [] - for filtered_complex in stars: + for filtered_simplex in star: v = [] - for vertex in filtered_complex.first: + for vertex in filtered_simplex.first: v.append(vertex) - ct.append((v, filtered_complex.second)) + ct.append((v, filtered_simplex.second)) return ct def get_cofaces(self, simplex, codimension): @@ -275,17 +267,17 @@ cdef class SimplexTree: :returns: The (simplices of the) cofaces of a simplex :rtype: list of tuples(simplex, filtration) """ - cdef vector[int] complex + cdef vector[int] csimplex for i in simplex: - complex.push_back(i) + csimplex.push_back(i) cdef vector[pair[vector[int], double]] cofaces \ - = self.thisptr.get_cofaces(complex, <int>codimension) + = self.thisptr.get_cofaces(csimplex, <int>codimension) ct = [] - for filtered_complex in cofaces: + for filtered_simplex in cofaces: v = [] - for vertex in filtered_complex.first: + for vertex in filtered_simplex.first: v.append(vertex) - ct.append((v, filtered_complex.second)) + ct.append((v, filtered_simplex.second)) return ct def remove_maximal_simplex(self, simplex): @@ -385,7 +377,7 @@ cdef class SimplexTree: complex in a specific dimension. :param dimension: The specific dimension. - :type from_value: int. + :type dimension: int. :returns: The persistence intervals. :rtype: list of pair of float @@ -399,3 +391,22 @@ cdef class SimplexTree: print("intervals_in_dim function requires persistence function" " to be launched first.") return intervals_result + + def write_persistence_diagram(self, persistence_file=''): + """This function writes the persistence intervals of the simplicial + complex in a user given file name. + + :param persistence_file: The specific dimension. + :type persistence_file: string. + + :note: intervals_in_dim function requires persistence function to be + launched first. + """ + if self.pcohptr != NULL: + if persistence_file != '': + self.pcohptr.write_output_diagram(str.encode(persistence_file)) + else: + print("persistence_file must be specified") + else: + print("intervals_in_dim function requires persistence function" + " to be launched first.") |