Merge remote-tracking branch 'origin/master' into tomato2

3 files changed, 217 insertions, 34 deletions

diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index 007abcb6..28fbe3af 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -27,19 +27,20 @@ __license__ = "MIT"
 
 cdef extern from "Cubical_complex_interface.h" namespace "Gudhi":
     cdef cppclass Bitmap_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface<>":
-        Bitmap_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells)
-        Bitmap_cubical_complex_base_interface(string perseus_file)
-        int num_simplices()
-        int dimension()
+        Bitmap_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells) nogil
+        Bitmap_cubical_complex_base_interface(string perseus_file) nogil
+        int num_simplices() nogil
+        int dimension() nogil
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<>>":
-        Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max)
-        void compute_persistence(int homology_coeff_field, double min_persistence)
-        vector[pair[int, pair[double, double]]] get_persistence()
-        vector[int] betti_numbers()
-        vector[int] persistent_betti_numbers(double from_value, double to_value)
-        vector[pair[double,double]] intervals_in_dimension(int dimension)
+        Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) nogil
+        void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+        vector[pair[int, pair[double, double]]] get_persistence() nogil
+        vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil
+        vector[int] betti_numbers() nogil
+        vector[int] persistent_betti_numbers(double from_value, double to_value) nogil
+        vector[pair[double,double]] intervals_in_dimension(int dimension) nogil
 
 # CubicalComplex python interface
 cdef class CubicalComplex:
@@ -79,7 +80,7 @@ cdef class CubicalComplex:
                   perseus_file=''):
         if ((dimensions is not None) and (top_dimensional_cells is not None)
             and (perseus_file == '')):
-            self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells)
+            self._construct_from_cells(dimensions, top_dimensional_cells)
         elif ((dimensions is None) and (top_dimensional_cells is not None)
             and (perseus_file == '')):
             top_dimensional_cells = np.array(top_dimensional_cells,
@@ -87,11 +88,11 @@ cdef class CubicalComplex:
                                              order = 'F')
             dimensions = top_dimensional_cells.shape
             top_dimensional_cells = top_dimensional_cells.ravel(order='F')
-            self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells)
+            self._construct_from_cells(dimensions, top_dimensional_cells)
         elif ((dimensions is None) and (top_dimensional_cells is None)
             and (perseus_file != '')):
             if os.path.isfile(perseus_file):
-                self.thisptr = new Bitmap_cubical_complex_base_interface(perseus_file.encode('utf-8'))
+                self._construct_from_file(perseus_file.encode('utf-8'))
             else:
                 raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT),
                                         perseus_file)
@@ -100,6 +101,14 @@ cdef class CubicalComplex:
                 "top_dimensional_cells or from a Perseus-style file name.",
                 file=sys.stderr)
 
+    def _construct_from_cells(self, vector[unsigned] dimensions, vector[double] top_dimensional_cells):
+        with nogil:
+            self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells)
+
+    def _construct_from_file(self, string filename):
+        with nogil:
+            self.thisptr = new Bitmap_cubical_complex_base_interface(filename)
+
     def __dealloc__(self):
         if self.thisptr != NULL:
             del self.thisptr
@@ -150,8 +159,11 @@ cdef class CubicalComplex:
         if self.pcohptr != NULL:
             del self.pcohptr
         assert self.__is_defined()
-        self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True)
-        self.pcohptr.compute_persistence(homology_coeff_field, min_persistence)
+        cdef int field = homology_coeff_field
+        cdef double minp = min_persistence
+        with nogil:
+            self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, 1)
+            self.pcohptr.compute_persistence(field, minp)
 
     def persistence(self, homology_coeff_field=11, min_persistence=0):
         """This function computes and returns the persistence of the complex.
@@ -170,6 +182,59 @@ cdef class CubicalComplex:
         self.compute_persistence(homology_coeff_field, min_persistence)
         return self.pcohptr.get_persistence()
 
+    def cofaces_of_persistence_pairs(self):
+        """A persistence interval is described by a pair of cells, one that creates the
+        feature and one that kills it. The filtration values of those 2 cells give coordinates
+        for a point in a persistence diagram, or a bar in a barcode. Structurally, in the
+        cubical complexes provided here, the filtration value of any cell is the minimum of the
+        filtration values of the maximal cells that contain it. Connecting persistence diagram
+        coordinates to the corresponding value in the input (i.e. the filtration values of
+        the top-dimensional cells) is useful for differentiation purposes.
+
+        This function returns a list of pairs of top-dimensional cells corresponding to
+        the persistence birth and death cells of the filtration. The cells are represented by
+        their indices in the input list of top-dimensional cells (and not their indices in the
+        internal datastructure that includes non-maximal cells). Note that when two adjacent
+        top-dimensional cells have the same filtration value, we arbitrarily return one of the two
+        when calling the function on one of their common faces.
+
+        :returns: The top-dimensional cells/cofaces of the positive and negative cells,
+            together with the corresponding homological dimension, in two lists of numpy arrays of integers.
+            The first list contains the regular persistence pairs, grouped by dimension.
+            It contains numpy arrays of shape [number_of_persistence_points, 2].
+            The indices of the arrays in the list correspond to the homological dimensions, and the
+            integers of each row in each array correspond to: (index of positive top-dimensional cell,
+            index of negative top-dimensional cell).
+            The second list contains the essential features, grouped by dimension.
+            It contains numpy arrays of shape [number_of_persistence_points, 1].
+            The indices of the arrays in the list correspond to the homological dimensions, and the
+            integers of each row in each array correspond to: (index of positive top-dimensional cell).
+        """
+
+        assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()"
+
+        cdef vector[vector[int]] persistence_result
+        output = [[],[]]
+        with nogil:
+            persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
+        pr = np.array(persistence_result)
+
+        ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
+        ess = pr[ess_ind]
+        max_h = max(ess[:,0])+1 if len(ess) > 0 else 0
+        for h in range(max_h):
+            hidxs = np.argwhere(ess[:,0] == h)[:,0]
+            output[1].append(ess[hidxs][:,1])
+
+        reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
+        reg = pr[reg_ind]
+        max_h = max(reg[:,0])+1 if len(reg) > 0 else 0
+        for h in range(max_h):
+            hidxs = np.argwhere(reg[:,0] == h)[:,0]
+            output[0].append(reg[hidxs][:,1:])
+
+        return output
+
     def betti_numbers(self):
         """This function returns the Betti numbers of the complex.
 
diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx
index 246a3a02..d353d2af 100644
--- a/src/python/gudhi/periodic_cubical_complex.pyx
+++ b/src/python/gudhi/periodic_cubical_complex.pyx
@@ -24,19 +24,20 @@ __license__ = "MIT"
 
 cdef extern from "Cubical_complex_interface.h" namespace "Gudhi":
     cdef cppclass Periodic_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface<Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>>":
-        Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions)
-        Periodic_cubical_complex_base_interface(string perseus_file)
-        int num_simplices()
-        int dimension()
+        Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions) nogil
+        Periodic_cubical_complex_base_interface(string perseus_file) nogil
+        int num_simplices() nogil
+        int dimension() nogil
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Periodic_cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>>>":
-        Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max)
-        void compute_persistence(int homology_coeff_field, double min_persistence)
-        vector[pair[int, pair[double, double]]] get_persistence()
-        vector[int] betti_numbers()
-        vector[int] persistent_betti_numbers(double from_value, double to_value)
-        vector[pair[double,double]] intervals_in_dimension(int dimension)
+        Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) nogil
+        void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+        vector[pair[int, pair[double, double]]] get_persistence() nogil
+        vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil
+        vector[int] betti_numbers() nogil
+        vector[int] persistent_betti_numbers(double from_value, double to_value) nogil
+        vector[pair[double,double]] intervals_in_dimension(int dimension) nogil
 
 # PeriodicCubicalComplex python interface
 cdef class PeriodicCubicalComplex:
@@ -80,9 +81,7 @@ cdef class PeriodicCubicalComplex:
                   periodic_dimensions=None, perseus_file=''):
         if ((dimensions is not None) and (top_dimensional_cells is not None)
             and (periodic_dimensions is not None) and (perseus_file == '')):
-            self.thisptr = new Periodic_cubical_complex_base_interface(dimensions,
-                                                                       top_dimensional_cells,
-                                                                       periodic_dimensions)
+            self._construct_from_cells(dimensions, top_dimensional_cells, periodic_dimensions)
         elif ((dimensions is None) and (top_dimensional_cells is not None)
             and (periodic_dimensions is not None) and (perseus_file == '')):
             top_dimensional_cells = np.array(top_dimensional_cells,
@@ -90,13 +89,11 @@ cdef class PeriodicCubicalComplex:
                                              order = 'F')
             dimensions = top_dimensional_cells.shape
             top_dimensional_cells = top_dimensional_cells.ravel(order='F')
-            self.thisptr = new Periodic_cubical_complex_base_interface(dimensions,
-                                                                       top_dimensional_cells,
-                                                                       periodic_dimensions)
+            self._construct_from_cells(dimensions, top_dimensional_cells, periodic_dimensions)
         elif ((dimensions is None) and (top_dimensional_cells is None)
             and (periodic_dimensions is None) and (perseus_file != '')):
             if os.path.isfile(perseus_file):
-                self.thisptr = new Periodic_cubical_complex_base_interface(perseus_file.encode('utf-8'))
+                self._construct_from_file(perseus_file.encode('utf-8'))
             else:
                 print("file " + perseus_file + " not found.", file=sys.stderr)
         else:
@@ -105,6 +102,14 @@ cdef class PeriodicCubicalComplex:
               "top_dimensional_cells and periodic_dimensions or from "
               "a Perseus-style file name.", file=sys.stderr)
 
+    def _construct_from_cells(self, vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions):
+        with nogil:
+            self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, top_dimensional_cells, periodic_dimensions)
+
+    def _construct_from_file(self, string filename):
+        with nogil:
+            self.thisptr = new Periodic_cubical_complex_base_interface(filename)
+
     def __dealloc__(self):
         if self.thisptr != NULL:
             del self.thisptr
@@ -155,8 +160,11 @@ cdef class PeriodicCubicalComplex:
         if self.pcohptr != NULL:
             del self.pcohptr
         assert self.__is_defined()
-        self.pcohptr = new Periodic_cubical_complex_persistence_interface(self.thisptr, True)
-        self.pcohptr.compute_persistence(homology_coeff_field, min_persistence)
+        cdef int field = homology_coeff_field
+        cdef double minp = min_persistence
+        with nogil:
+            self.pcohptr = new Periodic_cubical_complex_persistence_interface(self.thisptr, 1)
+            self.pcohptr.compute_persistence(field, minp)
 
     def persistence(self, homology_coeff_field=11, min_persistence=0):
         """This function computes and returns the persistence of the complex.
@@ -175,6 +183,57 @@ cdef class PeriodicCubicalComplex:
         self.compute_persistence(homology_coeff_field, min_persistence)
         return self.pcohptr.get_persistence()
 
+    def cofaces_of_persistence_pairs(self):
+        """A persistence interval is described by a pair of cells, one that creates the
+        feature and one that kills it. The filtration values of those 2 cells give coordinates
+        for a point in a persistence diagram, or a bar in a barcode. Structurally, in the
+        cubical complexes provided here, the filtration value of any cell is the minimum of the
+        filtration values of the maximal cells that contain it. Connecting persistence diagram
+        coordinates to the corresponding value in the input (i.e. the filtration values of
+        the top-dimensional cells) is useful for differentiation purposes.
+
+        This function returns a list of pairs of top-dimensional cells corresponding to
+        the persistence birth and death cells of the filtration. The cells are represented by
+        their indices in the input list of top-dimensional cells (and not their indices in the
+        internal datastructure that includes non-maximal cells). Note that when two adjacent
+        top-dimensional cells have the same filtration value, we arbitrarily return one of the two
+        when calling the function on one of their common faces.
+
+        :returns: The top-dimensional cells/cofaces of the positive and negative cells,
+            together with the corresponding homological dimension, in two lists of numpy arrays of integers.
+            The first list contains the regular persistence pairs, grouped by dimension.
+            It contains numpy arrays of shape [number_of_persistence_points, 2].
+            The indices of the arrays in the list correspond to the homological dimensions, and the
+            integers of each row in each array correspond to: (index of positive top-dimensional cell,
+            index of negative top-dimensional cell).
+            The second list contains the essential features, grouped by dimension.
+            It contains numpy arrays of shape [number_of_persistence_points, 1].
+            The indices of the arrays in the list correspond to the homological dimensions, and the
+            integers of each row in each array correspond to: (index of positive top-dimensional cell).
+        """
+        assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()"
+        cdef vector[vector[int]] persistence_result
+
+        output = [[],[]]
+        with nogil:
+            persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
+        pr = np.array(persistence_result)
+
+        ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
+        ess = pr[ess_ind]
+        max_h = max(ess[:,0])+1 if len(ess) > 0 else 0
+        for h in range(max_h):
+            hidxs = np.argwhere(ess[:,0] == h)[:,0]
+            output[1].append(ess[hidxs][:,1])
+
+        reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
+        reg = pr[reg_ind]
+        max_h = max(reg[:,0])+1 if len(reg) > 0 else 0
+        for h in range(max_h):
+            hidxs = np.argwhere(reg[:,0] == h)[:,0]
+            output[0].append(reg[hidxs][:,1:])
+        return output
+
     def betti_numbers(self):
         """This function returns the Betti numbers of the complex.
 
diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py
new file mode 100644
index 00000000..0541572b
--- /dev/null
+++ b/src/python/gudhi/weighted_rips_complex.py
@@ -0,0 +1,59 @@
+# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+# Author(s):       Raphaël Tinarrage, Yuichi Ike, Masatoshi Takenouchi
+#
+# Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd.
+#
+# Modification(s):
+#   - YYYY/MM Author: Description of the modification
+
+from gudhi import SimplexTree
+
+class WeightedRipsComplex:
+    """
+    Class to generate a weighted Rips complex from a distance matrix and weights on vertices, 
+    in the way described in :cite:`dtmfiltrations`.
+    Remark that all the filtration values are doubled compared to the definition in the paper 
+    for the consistency with RipsComplex.
+    """
+    def __init__(self, 
+                distance_matrix, 
+                weights=None,
+                max_filtration=float('inf')):
+        """
+        Args:
+            distance_matrix (Sequence[Sequence[float]]): distance matrix (full square or lower triangular).
+            weights (Sequence[float]): (one half of) weight for each vertex.
+            max_filtration (float): specifies the maximal filtration value to be considered.      
+        """
+        self.distance_matrix = distance_matrix
+        if weights is not None:
+            self.weights = weights
+        else:
+            self.weights = [0] * len(distance_matrix)
+        self.max_filtration = max_filtration
+            
+    def create_simplex_tree(self, max_dimension):
+        """
+        Args:
+            max_dimension (int): graph expansion until this given dimension.
+        """
+        dist = self.distance_matrix
+        F = self.weights
+        num_pts = len(dist)
+        
+        st = SimplexTree()
+        
+        for i in range(num_pts):
+            if 2*F[i] <= self.max_filtration:
+                st.insert([i], 2*F[i])
+        for i in range(num_pts):
+            for j in range(i):
+                value = max(2*F[i], 2*F[j], dist[i][j] + F[i] + F[j])
+                # max is needed when F is not 1-Lipschitz
+                if value <= self.max_filtration:
+                    st.insert([i,j], filtration=value)
+                    
+        st.expansion(max_dimension) 
+        return st
+