Merge pull request #295 from yuichi-ike/weighted_rips

Weighted rips

7 files changed, 222 insertions, 0 deletions

diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib
index 3ea2f59f..677ed4df 100644
--- a/biblio/bibliography.bib
+++ b/biblio/bibliography.bib
@@ -1253,3 +1253,29 @@ year = "2011"
   year={2014},
   publisher={Springer}
 }
+
+@inproceedings{dtmfiltrations,
+  author    = {Hirokazu Anai and
+               Fr{\'{e}}d{\'{e}}ric Chazal and
+               Marc Glisse and
+               Yuichi Ike and
+               Hiroya Inakoshi and
+               Rapha{\"{e}}l Tinarrage and
+               Yuhei Umeda},
+  editor    = {Gill Barequet and
+               Yusu Wang},
+  title     = {DTM-Based Filtrations},
+  booktitle = {35th International Symposium on Computational Geometry, SoCG 2019,
+               June 18-21, 2019, Portland, Oregon, {USA}},
+  series    = {LIPIcs},
+  volume    = {129},
+  pages     = {58:1--58:15},
+  publisher = {Schloss Dagstuhl - Leibniz-Zentrum f{\"{u}}r Informatik},
+  year      = {2019},
+  url       = {https://doi.org/10.4230/LIPIcs.SoCG.2019.58},
+  doi       = {10.4230/LIPIcs.SoCG.2019.58},
+  timestamp = {Tue, 11 Feb 2020 15:52:14 +0100},
+  biburl    = {https://dblp.org/rec/conf/compgeom/AnaiCGIITU19.bib},
+  bibsource = {dblp computer science bibliography, https://dblp.org}
+}
+
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 976a8b52..ab08cd6d 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -57,6 +57,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}")
@@ -232,6 +233,7 @@ if(PYTHONINTERP_FOUND)
     file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/")
     file(COPY "gudhi/wasserstein" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -488,6 +490,11 @@ if(PYTHONINTERP_FOUND)
       add_gudhi_py_test(test_dtm)
     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/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/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")]]))
+