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

12 files changed, 459 insertions, 45 deletions

diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index e6a78a6d..f8f80ded 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -13,6 +13,8 @@
 
 #include <gudhi/Bitmap_cubical_complex/counter.h>
 
+#include <boost/config.hpp>
+
 #include <iostream>
 #include <vector>
 #include <string>
@@ -110,6 +112,16 @@ class Bitmap_cubical_complex_base {
   virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
 
   /**
+   * This function finds a top-dimensional cell that is incident to the input cell and has
+   * the same filtration value. In case several cells are suitable, an arbitrary one is
+   * returned. Note that the input parameter can be a cell of any dimension (vertex, edge, etc).
+   * On the other hand, the output is always indicating the position of
+   * a top-dimensional cube in the data structure.
+   * \pre The filtration values are assigned as per `impose_lower_star_filtration()`.
+   **/
+  inline size_t get_top_dimensional_coface_of_a_cell(size_t splx);
+
+  /**
   * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of
   * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence
   * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$.
@@ -603,6 +615,19 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells
 }
 
 template <typename T>
+size_t Bitmap_cubical_complex_base<T>::get_top_dimensional_coface_of_a_cell(size_t splx) {
+  if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;}
+  else{
+    for (auto v : this->get_coboundary_of_a_cell(splx)){
+      if(this->get_cell_data(v) == this->get_cell_data(splx)){
+        return this->get_top_dimensional_coface_of_a_cell(v);
+      }
+    }
+  }
+  BOOST_UNREACHABLE_RETURN(-2);
+}
+
+template <typename T>
 Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions,
                                                             const std::vector<T>& top_dimensional_cells) {
   this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells);
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 7186dc2b..a8872d73 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -58,6 +58,7 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'representations', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'point_cloud', ")
+    set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'weighted_rips_complex', ")
 
     add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}")
     add_gudhi_debug_info("Cython version ${CYTHON_VERSION}")
@@ -235,6 +236,7 @@ if(PYTHONINTERP_FOUND)
     file(COPY "gudhi/wasserstein" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/clustering" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi" FILES_MATCHING PATTERN "*.py")
+    file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -496,6 +498,11 @@ if(PYTHONINTERP_FOUND)
       add_gudhi_py_test(test_tomato)
     endif()
 
+    # Weighted Rips
+    if(SCIPY_FOUND)
+      add_gudhi_py_test(test_weighted_rips)
+    endif()
+
     # Set missing or not modules
     set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES")
   else(CYTHON_FOUND)
diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst
index 22b5616c..5f3e46c1 100644
--- a/src/python/doc/rips_complex_ref.rst
+++ b/src/python/doc/rips_complex_ref.rst
@@ -12,3 +12,16 @@ Rips complex reference manual
    :show-inheritance:
 
    .. automethod:: gudhi.RipsComplex.__init__
+
+.. _weighted-rips-complex-reference-manual:
+
+======================================
+Weighted Rips complex reference manual
+======================================
+
+.. autoclass:: gudhi.weighted_rips_complex.WeightedRipsComplex
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+   .. automethod:: gudhi.weighted_rips_complex.WeightedRipsComplex.__init__
diff --git a/src/python/doc/rips_complex_sum.inc b/src/python/doc/rips_complex_sum.inc
index 6feb74cd..f7580714 100644
--- a/src/python/doc/rips_complex_sum.inc
+++ b/src/python/doc/rips_complex_sum.inc
@@ -11,6 +11,9 @@
    |                                                                |                                                                        |                                                                      |
    |                                                                | This complex can be built from a point cloud and a distance function,  |                                                                      |
    |                                                                | or from a distance matrix.                                             |                                                                      |
+   |                                                                |                                                                        |                                                                      |
+   |                                                                | Weighted Rips complex constructs a simplicial complex from a distance  |                                                                      |
+   |                                                                | matrix and weights on vertices.                                        |                                                                      |
    +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+
    | * :doc:`rips_complex_user`                                     | * :doc:`rips_complex_ref`                                                                                                                     |
    +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------+
diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst
index 8efb12e6..819568be 100644
--- a/src/python/doc/rips_complex_user.rst
+++ b/src/python/doc/rips_complex_user.rst
@@ -347,3 +347,54 @@ until dimension 1 - one skeleton graph in other words), the output is:
     points in the persistence diagram will be under the diagonal, and
     bottleneck distance and persistence graphical tool will not work properly,
     this is a known issue.
+
+Weighted Rips Complex
+---------------------
+
+`WeightedRipsComplex <rips_complex_ref.html#weighted-rips-complex-reference-manual>`_ builds a simplicial complex from a distance matrix and weights on vertices.
+
+
+Example from a distance matrix and weights
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The following example computes the weighted Rips filtration associated with a distance matrix and weights on vertices.
+
+.. testcode::
+
+    from gudhi.weighted_rips_complex import WeightedRipsComplex
+    dist = [[], [1]]
+    weights = [1, 100]
+    w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights)
+    st = w_rips.create_simplex_tree(max_dimension=2)
+    print(list(st.get_filtration()))
+
+The output is:
+
+.. testoutput::
+
+    [([0], 2.0), ([1], 200.0), ([0, 1], 200.0)]
+
+Example from a point cloud combined with DistanceToMeasure
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-DTM-filtrations.ipynb>`_. 
+
+.. testcode::
+
+    import numpy as np
+    from scipy.spatial.distance import cdist
+    from gudhi.point_cloud.dtm import DistanceToMeasure
+    from gudhi.weighted_rips_complex import WeightedRipsComplex
+    pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]])
+    dist = cdist(pts,pts)
+    dtm = DistanceToMeasure(2, q=2, metric="precomputed")
+    r = dtm.fit_transform(dist)
+    w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r)
+    st = w_rips.create_simplex_tree(max_dimension=2)
+    print(st.persistence())
+
+The output is:
+
+.. testoutput::
+
+    [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))]
diff --git a/src/python/doc/wasserstein_distance_sum.inc b/src/python/doc/wasserstein_distance_sum.inc
index f9308e5e..c41de017 100644
--- a/src/python/doc/wasserstein_distance_sum.inc
+++ b/src/python/doc/wasserstein_distance_sum.inc
@@ -1,14 +1,12 @@
 .. table::
    :widths: 30 40 30
 
-   +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+
-   | .. figure::                                                     | The q-Wasserstein distance measures the similarity between two       | :Author: Theo Lacombe                                            |
-   |      ../../doc/Bottleneck_distance/perturb_pd.png               | persistence diagrams using the sum of all edges lengths (instead of  |                                                                  |
-   |      :figclass: align-center                                    | the maximum). It allows to define sophisticated objects such as      | :Since: GUDHI 3.1.0                                              |
-   |                                                                 | barycenters of a family of persistence diagrams.                     |                                                                  |
-   |                                                                 |                                                                      | :License: MIT                                                    |
-   |                                                                 |                                                                      |                                                                  |
-   |                                                                 |                                                                      | :Requires: Python Optimal Transport (POT) :math:`\geq` 0.5.1     |
-   +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+
-   | * :doc:`wasserstein_distance_user`                              |                                                                                                                                         |
-   +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------+
+   +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------+
+   | .. figure::                                                     | The q-Wasserstein distance measures the similarity between two       | :Author: Theo Lacombe, Marc Glisse      |
+   |      ../../doc/Bottleneck_distance/perturb_pd.png               | persistence diagrams using the sum of all edges lengths (instead of  |                                         |
+   |      :figclass: align-center                                    | the maximum). It allows to define sophisticated objects such as      | :Since: GUDHI 3.1.0                     |
+   |                                                                 | barycenters of a family of persistence diagrams.                     |                                         |
+   |                                                                 |                                                                      | :License: MIT, BSD-3-Clause             |
+   +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------+
+   | * :doc:`wasserstein_distance_user`                              |                                                                                                                |
+   +-----------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------+
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
+        
diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h
index 0de9bd5c..e5a3dfba 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -13,9 +13,11 @@
 
 #include <gudhi/Persistent_cohomology.h>
 
+#include <cstdlib>
 #include <vector>
 #include <utility>  // for std::pair
 #include <algorithm>  // for sort
+#include <unordered_map>
 
 namespace Gudhi {
 
@@ -67,6 +69,48 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     return persistence;
   }
 
+  // This function computes the top-dimensional cofaces associated to the positive and negative
+  // simplices of a cubical complex. The output format is a vector of vectors of three integers,
+  // which are [homological dimension, index of top-dimensional coface of positive simplex,
+  // index of top-dimensional coface of negative simplex]. If the topological feature is essential,
+  // then the index of top-dimensional coface of negative simplex is arbitrarily set to -1.
+  std::vector<std::vector<int>> cofaces_of_cubical_persistence_pairs() {
+
+    // Warning: this function is meant to be used with CubicalComplex only!!
+
+    auto&& pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
+      persistent_cohomology::Field_Zp>::get_persistent_pairs();
+
+    // Gather all top-dimensional cells and store their simplex handles
+    std::vector<std::size_t> max_splx;
+    for (auto splx : stptr_->top_dimensional_cells_range())
+      max_splx.push_back(splx);
+    // Sort these simplex handles and compute the ordering function
+    // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data
+    std::unordered_map<std::size_t, int> order;
+    //std::sort(max_splx.begin(), max_splx.end());
+    for (unsigned int i = 0; i < max_splx.size(); i++)  order.emplace(max_splx[i], i);
+
+    std::vector<std::vector<int>> persistence_pairs;
+    for (auto pair : pairs) {
+      int h = stptr_->dimension(get<0>(pair));
+      // Recursively get the top-dimensional cell / coface associated to the persistence generator
+      std::size_t face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair));
+      // Retrieve the index of the corresponding top-dimensional cell in the input data
+      int splx0 = order[face0];
+
+      int splx1 = -1;
+      if (get<1>(pair) != stptr_->null_simplex()){
+        // Recursively get the top-dimensional cell / coface associated to the persistence generator
+        std::size_t face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair));
+        // Retrieve the index of the corresponding top-dimensional cell in the input data
+        splx1 = order[face1];
+      }
+      persistence_pairs.push_back({ h, splx0, splx1 });
+    }
+    return persistence_pairs;
+  }
+
   std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() {
     std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs;
     auto const& pairs = Base::get_persistent_pairs();
diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py
index fce4875c..d0e4e9e8 100755
--- a/src/python/test/test_cubical_complex.py
+++ b/src/python/test/test_cubical_complex.py
@@ -147,3 +147,30 @@ def test_connected_sublevel_sets():
                                  periodic_dimensions = periodic_dimensions)
     assert cub.persistence() == [(0, (2.0, float("inf")))]
     assert cub.betti_numbers() == [1, 0, 0]
+
+def test_cubical_generators():
+    cub = CubicalComplex(top_dimensional_cells = [[0, 0, 0], [0, 1, 0], [0, 0, 0]])
+    cub.persistence()
+    g = cub.cofaces_of_persistence_pairs()
+    assert len(g[0]) == 2
+    assert len(g[1]) == 1
+    assert np.array_equal(g[0][0], np.empty(shape=[0,2]))
+    assert np.array_equal(g[0][1], np.array([[7, 4]]))
+    assert np.array_equal(g[1][0], np.array([8]))
+
+def test_cubical_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_death():
+    cubCpx = CubicalComplex(dimensions=[1,2], top_dimensional_cells=[0.0, 1.0])
+    Diag = cubCpx.persistence(homology_coeff_field=2, min_persistence=0)
+    pairs = cubCpx.cofaces_of_persistence_pairs()
+    assert pairs[0] == []
+    assert np.array_equal(pairs[1][0], np.array([0]))
+
+def test_periodic_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_death():
+    perCubCpx = PeriodicCubicalComplex(dimensions=[1,2], top_dimensional_cells=[0.0, 1.0],
+                                       periodic_dimensions=[True, True])
+    Diag = perCubCpx.persistence(homology_coeff_field=2, min_persistence=0)
+    pairs = perCubCpx.cofaces_of_persistence_pairs()
+    assert pairs[0] == []
+    assert np.array_equal(pairs[1][0], np.array([0]))
+    assert np.array_equal(pairs[1][1], np.array([0, 1]))
+    assert np.array_equal(pairs[1][2], np.array([1]))
diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips.py
new file mode 100644
index 00000000..7ef48333
--- /dev/null
+++ b/src/python/test/test_weighted_rips.py
@@ -0,0 +1,63 @@
+""" 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):       Yuichi Ike and Masatoshi Takenouchi
+
+    Copyright (C) 2020 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+from gudhi.weighted_rips_complex import WeightedRipsComplex
+from gudhi.point_cloud.dtm import DistanceToMeasure
+import numpy as np
+from math import sqrt
+from scipy.spatial.distance import cdist
+import pytest
+
+def test_non_dtm_rips_complex():
+    dist = [[], [1]]
+    weights = [1, 100]
+    w_rips = WeightedRipsComplex(distance_matrix=dist, weights=weights)
+    st = w_rips.create_simplex_tree(max_dimension=2)
+    assert st.filtration([0,1]) == pytest.approx(200.0)
+
+def test_compatibility_with_rips():
+    distance_matrix = [[0], [1, 0], [1, sqrt(2), 0], [sqrt(2), 1, 1, 0]]
+    w_rips = WeightedRipsComplex(distance_matrix=distance_matrix,max_filtration=42)
+    st = w_rips.create_simplex_tree(max_dimension=1)
+    assert list(st.get_filtration()) == [
+        ([0], 0.0),
+        ([1], 0.0),
+        ([2], 0.0),
+        ([3], 0.0),
+        ([0, 1], 1.0),
+        ([0, 2], 1.0),
+        ([1, 3], 1.0),
+        ([2, 3], 1.0),
+        ([1, 2], sqrt(2)),
+        ([0, 3], sqrt(2)),
+    ]
+
+def test_compatibility_with_filtered_rips():
+    distance_matrix = [[0], [1, 0], [1, sqrt(2), 0], [sqrt(2), 1, 1, 0]]
+    w_rips = WeightedRipsComplex(distance_matrix=distance_matrix,max_filtration=1.0)
+    st = w_rips.create_simplex_tree(max_dimension=1)
+    
+    assert st.__is_defined() == True
+    assert st.__is_persistence_defined() == False
+
+    assert st.num_simplices() == 8
+    assert st.num_vertices() == 4
+
+def test_dtm_rips_complex():
+    pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]])
+    dist = cdist(pts,pts)
+    dtm = DistanceToMeasure(2, q=2, metric="precomputed")
+    r = dtm.fit_transform(dist)
+    w_rips = WeightedRipsComplex(distance_matrix=dist, weights=r)
+    st = w_rips.create_simplex_tree(max_dimension=2)
+    st.persistence()
+    persistence_intervals0 = st.persistence_intervals_in_dimension(0)
+    assert persistence_intervals0 == pytest.approx(np.array([[3.16227766, 5.39834564],[3.16227766, 5.39834564], [3.16227766, float("inf")]]))
+