Merge branch sparserips-python

Fix doxygen bad link on sparsify point set function
Fix some lines too long


git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@3961 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 94946e2dd7853bd65d6790546ad2ef8a35d298c1

10 files changed, 185 insertions, 62 deletions

diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib
index 16f43d6f..a8c040ed 100644
--- a/biblio/bibliography.bib
+++ b/biblio/bibliography.bib
@@ -1073,9 +1073,10 @@ language={English}
 @article{buchet16efficient,
   title = {Efficient and Robust Persistent Homology for Measures},
   author = {Micka\"{e}l Buchet and Fr\'{e}d\'{e}ric Chazal and Steve Y. Oudot and Donald Sheehy},
-  booktitle = {Computational Geometry: Theory and Applications},
+  journal = {Computational Geometry: Theory and Applications},
   volume = {58},
   pages = {70--96},
+  doi = "https://doi.org/10.1016/j.comgeo.2016.07.001",
   year = {2016}
 }
 
@@ -1083,6 +1084,7 @@ language={English}
   author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy},
   booktitle = {Proceedings of the Canadian Conference on Computational Geometry},
   title = {A Geometric Perspective on Sparse Filtrations},
+  url = {http://research.cs.queensu.ca/cccg2015/CCCG15-papers/01.pdf},
   year = {2015}
 }
 
@@ -1101,5 +1103,6 @@ language={English}
   volume = {49},
   number = {4},
   pages = {778--796},
+  doi = "10.1007/s00454-013-9513-1",
   year = {2013}
 }
diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h
index 712d3b6e..a2537036 100644
--- a/src/Rips_complex/doc/Intro_rips_complex.h
+++ b/src/Rips_complex/doc/Intro_rips_complex.h
@@ -39,7 +39,7 @@ namespace rips_complex {
  * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a>
  * is an abstract simplicial complex
  * defined on a finite metric space, where each simplex corresponds to a subset
- * of point whose diameter is smaller that some given threshold.
+ * of points whose diameter is smaller that some given threshold.
  * Varying the threshold, we can also see the Rips complex as a filtration of
  * the \f$(n-1)-\f$dimensional simplex, where the filtration value of each
  * simplex is the diameter of the corresponding subset of points.
@@ -53,7 +53,10 @@ namespace rips_complex {
  * The number of simplices in the full Rips complex is exponential in the
  * number of vertices, it is thus usually restricted, by excluding all the
  * simplices with filtration value larger than some threshold, and keeping only
- * the dim_max-skeleton.
+ * the dim_max-skeleton. It may also be a good idea to start by making the
+ * point set sparser, for instance with
+ * `Gudhi::subsampling::sparsify_point_set()`, since small clusters of points
+ * have a disproportionate cost without affecting the persistence diagram much.
  *
  * In order to build this complex, the algorithm first builds the graph.
  * The filtration value of each edge is computed from a user-given distance
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 3ab23c12..3339606c 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -1035,6 +1035,7 @@ class Simplex_tree {
    * The Simplex_tree must contain no simplex of dimension bigger than
    * 1 when calling the method. */
   void expansion(int max_dim) {
+    if (max_dim <= 1) return;
     dimension_ = max_dim;
     for (Dictionary_it root_it = root_.members_.begin();
          root_it != root_.members_.end(); ++root_it) {
diff --git a/src/cython/cython/rips_complex.pyx b/src/cython/cython/rips_complex.pyx
index 30ca4443..7c83241c 100644
--- a/src/cython/cython/rips_complex.pyx
+++ b/src/cython/cython/rips_complex.pyx
@@ -33,7 +33,11 @@ __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)
+        Rips_complex_interface()
+        void init_points(vector[vector[double]] values, double threshold)
+        void init_matrix(vector[vector[double]] values, double threshold)
+        void init_points_sparse(vector[vector[double]] values, double threshold, double sparse)
+        void init_matrix_sparse(vector[vector[double]] values, double threshold, double sparse)
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, int dim_max)
 
 # RipsComplex python interface
@@ -44,10 +48,11 @@ cdef class RipsComplex:
     function, or a distance matrix.
     """
 
-    cdef Rips_complex_interface * thisptr
+    cdef Rips_complex_interface thisref
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=None, distance_matrix=None, max_edge_length=float('inf')):
+    def __init__(self, points=None, distance_matrix=None,
+                 max_edge_length=float('inf'), sparse=None):
         """RipsComplex constructor.
 
         :param max_edge_length: Rips value.
@@ -59,29 +64,38 @@ cdef class RipsComplex:
         Or
 
         :param distance_matrix: A distance matrix (full square or lower
-        triangular).
+            triangular).
         :type points: list of list of double
+
+        And in both cases
+
+        :param sparse: If this is not None, it switches to building a sparse
+            Rips and represents the approximation parameter epsilon.
+        :type sparse: float
         """
 
     # The real cython constructor
-    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)
+    def __cinit__(self, points=None, distance_matrix=None,
+                  max_edge_length=float('inf'), sparse=None):
+        if sparse is not None:
+          if distance_matrix is not None:
+              self.thisref.init_matrix_sparse(distance_matrix,
+                                              max_edge_length,
+                                              sparse)
+          else:
+              if points is None:
+                  # Empty Rips construction
+                  points=[]
+              self.thisref.init_points_sparse(points, max_edge_length, sparse)
         else:
-            if points is None:
-                # Empty Rips construction
-                points=[]
-            self.thisptr = new Rips_complex_interface(points, max_edge_length, True)
-
-
-    def __dealloc__(self):
-        if self.thisptr != NULL:
-            del self.thisptr
+          if distance_matrix is not None:
+              self.thisref.init_matrix(distance_matrix, max_edge_length)
+          else:
+              if points is None:
+                  # Empty Rips construction
+                  points=[]
+              self.thisref.init_points(points, max_edge_length)
 
-    def __is_defined(self):
-        """Returns true if RipsComplex pointer is not NULL.
-         """
-        return self.thisptr != NULL
 
     def create_simplex_tree(self, max_dimension=1):
         """
@@ -92,5 +106,5 @@ cdef class RipsComplex:
         :rtype: SimplexTree
         """
         simplex_tree = SimplexTree()
-        self.thisptr.create_simplex_tree(simplex_tree.thisptr, max_dimension)
+        self.thisref.create_simplex_tree(simplex_tree.thisptr, max_dimension)
         return simplex_tree
diff --git a/src/cython/doc/examples.rst b/src/cython/doc/examples.rst
index 1f02f8a2..edbc2f72 100644
--- a/src/cython/doc/examples.rst
+++ b/src/cython/doc/examples.rst
@@ -22,6 +22,7 @@ Examples
     * :download:`rips_complex_diagram_persistence_from_off_file_example.py <../example/rips_complex_diagram_persistence_from_off_file_example.py>`
     * :download:`rips_complex_diagram_persistence_from_distance_matrix_file_example.py <../example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py>`
     * :download:`rips_persistence_diagram.py <../example/rips_persistence_diagram.py>`
+    * :download:`sparse_rips_persistence_diagram.py <../example/sparse_rips_persistence_diagram.py>`
     * :download:`random_cubical_complex_persistence_example.py <../example/random_cubical_complex_persistence_example.py>`
     * :download:`coordinate_graph_induced_complex.py <../example/coordinate_graph_induced_complex.py>`
     * :download:`functional_graph_induced_complex.py <../example/functional_graph_induced_complex.py>`
diff --git a/src/cython/doc/rips_complex_sum.inc b/src/cython/doc/rips_complex_sum.inc
index 5616bfa9..ea26769a 100644
--- a/src/cython/doc/rips_complex_sum.inc
+++ b/src/cython/doc/rips_complex_sum.inc
@@ -1,6 +1,6 @@
-=================================================================  ===================================  ===================================
-:Author: Clément Maria, Pawel Dlotko, Vincent Rouvreau             :Introduced in: GUDHI 2.0.0          :Copyright: GPL v3
-=================================================================  ===================================  ===================================
+=====================================================================  ===========================  ===================================
+:Author: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse    :Introduced in: GUDHI 2.0.0          :Copyright: GPL v3
+=====================================================================  ===========================  ===================================
 
 +----------------------------------------------------------------+------------------------------------------------------------------------+
 | .. figure::                                                    | Rips complex is a simplicial complex constructed from a one skeleton   |
diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst
index a8c06cf9..e814b4c3 100644
--- a/src/cython/doc/rips_complex_user.rst
+++ b/src/cython/doc/rips_complex_user.rst
@@ -7,27 +7,27 @@ Rips complex user manual
 Definition
 ----------
 
-=======================================================  =====================================  =====================================
-:Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau  :Introduced in: GUDHI 2.0.0            :Copyright: GPL v3
-=======================================================  =====================================  =====================================
+====================================================================  ================================  ======================
+:Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse  :Introduced in: GUDHI 2.0.0       :Copyright: GPL v3
+====================================================================  ================================  ======================
 
 +-------------------------------------------+----------------------------------------------------------------------+
 | :doc:`rips_complex_user`                  | :doc:`rips_complex_ref`                                              |
 +-------------------------------------------+----------------------------------------------------------------------+
 
-`Rips complex <https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex>`_ is a one skeleton graph that allows to
-construct a simplicial complex from it. The input can be a point cloud with a given distance function, or a distance
-matrix.
+The `Rips complex <https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex>`_ is a simplicial complex that
+generalizes proximity (:math:`\varepsilon`-ball) graphs to higher dimensions. The vertices correspond to the input
+points, and a simplex is present if and only if its diameter is smaller than some parameter α.  Considering all
+parameters α defines a filtered simplicial complex, where the filtration value of a simplex is its diameter.
+The filtration can be restricted to values α smaller than some threshold, to reduce its size.
 
-The filtration value of each edge is computed from a user-given distance function, or directly from the distance
-matrix.
+The input discrete metric space can be provided as a point cloud plus a distance function, or as a distance matrix.
 
-All edges that have a filtration value strictly greater than a given threshold value are not inserted into the complex.
+When creating a simplicial complex from the graph, :doc:`RipsComplex <rips_complex_ref>` first builds the graph and
+inserts it into the data structure. It then expands the simplicial complex (adds the simplices corresponding to cliques)
+when required. The expansion can be stopped at dimension `max_dimension`, by default 1.
 
-When creating a simplicial complex from this one skeleton graph, Rips inserts the one skeleton graph into the data
-structure, and then expands the simplicial complex when required.
-
-Vertex name correspond to the index of the point in the given range (aka. the point cloud).
+A vertex name corresponds to the index of the point in the given range (aka. the point cloud).
 
 .. figure::
     ../../doc/Rips_complex/rips_complex_representation.png
@@ -38,8 +38,27 @@ Vertex name correspond to the index of the point in the given range (aka. the po
 On this example, as edges (4,5), (4,6) and (5,6) are in the complex, simplex (4,5,6) is added with the filtration value
 set with :math:`max(filtration(4,5), filtration(4,6), filtration(5,6))`. And so on for simplex (0,1,2,3).
 
-If the Rips_complex interfaces are not detailed enough for your need, please refer to rips_persistence_step_by_step.cpp
-example, where the graph construction over the Simplex_tree is more detailed.
+If the `RipsComplex` interfaces are not detailed enough for your need, please refer to rips_persistence_step_by_step.cpp
+C++ example, where the graph construction over the Simplex_tree is more detailed.
+
+A Rips complex can easily become huge, even if we limit the length of the edges
+and the dimension of the simplices. One easy trick, before building a Rips
+complex on a point cloud, is to call `sparsify_point_set` which removes points
+that are too close to each other. This does not change its persistence diagram
+by more than the length used to define "too close".
+
+A more general technique is to use a sparse approximation of the Rips
+introduced by Don Sheehy :cite:`sheehy13linear`. We are using the version
+described in :cite:`buchet16efficient` (except that we multiply all filtration
+values by 2, to match the usual Rips complex), which proves a
+:math:`\frac{1+\varepsilon}{1-\varepsilon}`-interleaving, although in practice the
+error is usually smaller.  A more intuitive presentation of the idea is
+available in :cite:`cavanna15geometric`, and in a video
+:cite:`cavanna15visualizing`. Passing an extra argument `sparse=0.3` at the
+construction of a `RipsComplex` object asks it to build a sparse Rips with
+parameter :math:`\varepsilon=0.3`, while the default `sparse=None` builds the
+regular Rips complex.
+
 
 Point cloud
 -----------
@@ -47,7 +66,7 @@ Point cloud
 Example from a point cloud
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the one skeleton graph from the given points, and max_edge_length value.
+This example builds the neighborhood graph from the given points, up to max_edge_length.
 Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it.
 
 Finally, it is asked to display information about the simplicial complex.
@@ -56,7 +75,7 @@ Finally, it is asked to display information about the simplicial complex.
 
     import gudhi
     rips_complex = gudhi.RipsComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]],
-        max_edge_length=12.0)
+                                     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()) + ' - ' + \
@@ -92,10 +111,20 @@ until dimension 1 - one skeleton graph in other words), the output is:
     [4, 6] -> 9.49
     [3, 6] -> 11.00
 
+Notice that if we use
+
+.. code-block:: python
+
+    rips_complex = gudhi.RipsComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]],
+                                     max_edge_length=12.0, sparse=2)
+
+asking for a very sparse version (theory only gives some guarantee on the meaning of the output if `sparse<1`),
+2 to 5 edges disappear, depending on the random vertex used to start the sparsification.
+
 Example from OFF file
 ^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the :doc:`Rips_complex <rips_complex_ref>` from the given
+This example builds the :doc:`RipsComplex <rips_complex_ref>` from the given
 points in an OFF file, and max_edge_length value.
 Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it.
 
@@ -200,7 +229,7 @@ until dimension 1 - one skeleton graph in other words), the output is:
 Example from csv file
 ^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the :doc:`Rips_complex <rips_complex_ref>` from the given
+This example builds the :doc:`RipsComplex <rips_complex_ref>` from the given
 distance matrix in a csv file, and max_edge_length value.
 Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it.
 
diff --git a/src/cython/example/sparse_rips_persistence_diagram.py b/src/cython/example/sparse_rips_persistence_diagram.py
new file mode 100755
index 00000000..d58c244c
--- /dev/null
+++ b/src/cython/example/sparse_rips_persistence_diagram.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+
+import gudhi
+
+"""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):       Marc Glisse
+
+   Copyright (C) 2018 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__ = "Marc Glisse"
+__copyright__ = "Copyright (C) 2018 Inria"
+__license__ = "GPL v3"
+
+print("#####################################################################")
+print("Sparse RipsComplex creation from points")
+rips = gudhi.RipsComplex(points=[[0, 0], [0, 0.1], [1, 0], [0, 1], [1, 1]],
+                         max_edge_length=42, sparse=.5)
+
+simplex_tree = rips.create_simplex_tree(max_dimension=2)
+
+
+diag = simplex_tree.persistence(homology_coeff_field=2, min_persistence=0)
+print("diag=", diag)
+
+pplot = gudhi.plot_persistence_diagram(diag)
+pplot.show()
diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h
index 8b6c9c35..1a6e2477 100644
--- a/src/cython/include/Rips_complex_interface.h
+++ b/src/cython/include/Rips_complex_interface.h
@@ -25,8 +25,11 @@
 
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Rips_complex.h>
+#include <gudhi/Sparse_rips_complex.h>
 #include <gudhi/distance_functions.h>
 
+#include <boost/optional.hpp>
+
 #include "Simplex_tree_interface.h"
 
 #include <iostream>
@@ -43,28 +46,40 @@ class Rips_complex_interface {
   using Distance_matrix = std::vector<std::vector<Simplex_tree_interface<>::Filtration_value>>;
 
  public:
-  Rips_complex_interface(const std::vector<std::vector<double>>& values, double threshold, bool euclidean) {
-    if (euclidean) {
-      // Rips construction where values is a vector of points
-      rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold,
-                                                                                   Gudhi::Euclidean_distance());
-    } else {
-      // Rips construction where values is a distance matrix
-      rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold);
-    }
+  void init_points(const std::vector<std::vector<double>>& points, double threshold) {
+    rips_complex_.emplace(points, threshold, Gudhi::Euclidean_distance());
+  }
+  void init_matrix(const std::vector<std::vector<double>>& matrix, double threshold) {
+    rips_complex_.emplace(matrix, threshold);
   }
 
-  ~Rips_complex_interface() {
-    delete rips_complex_;
+  void init_points_sparse(const std::vector<std::vector<double>>& points, double threshold, double epsilon) {
+    sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon);
+    threshold_ = threshold;
+  }
+  void init_matrix_sparse(const std::vector<std::vector<double>>& matrix, double threshold, double epsilon) {
+    sparse_rips_complex_.emplace(matrix, epsilon);
+    threshold_ = threshold;
   }
 
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, int dim_max) {
-    rips_complex_->create_complex(*simplex_tree, dim_max);
+    if (rips_complex_)
+      rips_complex_->create_complex(*simplex_tree, dim_max);
+    else {
+      sparse_rips_complex_->create_complex(*simplex_tree, dim_max);
+      // This pruning should be done much earlier! It isn't that useful for sparse Rips,
+      // but it would be inconsistent not to do it.
+      simplex_tree->prune_above_filtration(threshold_);
+    }
     simplex_tree->initialize_filtration();
   }
 
  private:
-  Rips_complex<Simplex_tree_interface<>::Filtration_value>* rips_complex_;
+  // std::variant would work, but we don't require C++17 yet, and boost::variant is not super convenient.
+  // Anyway, storing a graph would make more sense. Or changing the interface completely so there is no such storage.
+  boost::optional<Rips_complex<Simplex_tree_interface<>::Filtration_value>> rips_complex_;
+  boost::optional<Sparse_rips_complex<Simplex_tree_interface<>::Filtration_value>> sparse_rips_complex_;
+  double threshold_;
 };
 
 }  // namespace rips_complex
diff --git a/src/cython/test/test_rips_complex.py b/src/cython/test/test_rips_complex.py
index c37b5400..05dfcaf7 100755
--- a/src/cython/test/test_rips_complex.py
+++ b/src/cython/test/test_rips_complex.py
@@ -30,7 +30,6 @@ __license__ = "GPL v3"
 
 def test_empty_rips():
     rips_complex = RipsComplex()
-    assert rips_complex.__is_defined() == True
 
 def test_rips_from_points():
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
@@ -68,12 +67,26 @@ def test_filtered_rips_from_points():
     assert simplex_tree.num_simplices() == 8
     assert simplex_tree.num_vertices() == 4
 
+def test_sparse_filtered_rips_from_points():
+    point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
+    filtered_rips = RipsComplex(points=point_list, max_edge_length=1.0,
+                                sparse=.001)
+
+    simplex_tree = filtered_rips.create_simplex_tree(max_dimension=1)
+
+    assert simplex_tree.__is_defined() == True
+    assert simplex_tree.__is_persistence_defined() == False
+
+    assert simplex_tree.num_simplices() == 8
+    assert simplex_tree.num_vertices() == 4
+
 def test_rips_from_distance_matrix():
     distance_matrix = [[0],
                        [1, 0],
                        [1, sqrt(2), 0],
                        [sqrt(2), 1, 1, 0]]
-    rips_complex = RipsComplex(distance_matrix=distance_matrix, max_edge_length=42)
+    rips_complex = RipsComplex(distance_matrix=distance_matrix,
+                               max_edge_length=42)
 
     simplex_tree = rips_complex.create_simplex_tree(max_dimension=1)
 
@@ -100,7 +113,8 @@ def test_filtered_rips_from_distance_matrix():
                        [1, 0],
                        [1, sqrt(2), 0],
                        [sqrt(2), 1, 1, 0]]
-    filtered_rips = RipsComplex(distance_matrix=distance_matrix, max_edge_length=1.0)
+    filtered_rips = RipsComplex(distance_matrix=distance_matrix,
+                                max_edge_length=1.0)
 
     simplex_tree = filtered_rips.create_simplex_tree(max_dimension=1)