plot persistence density functionnality with its documentation

Fix some documentation issue in other graphiacl tools
Add scipy as a dependency


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

Former-commit-id: 6481f1d6aa176bcf587f8875298a91cd22ac0319

4 files changed, 80 insertions, 5 deletions

diff --git a/src/cython/cython/persistence_graphical_tools.py b/src/cython/cython/persistence_graphical_tools.py
index 216ab8d6..7f604368 100755
--- a/src/cython/cython/persistence_graphical_tools.py
+++ b/src/cython/cython/persistence_graphical_tools.py
@@ -1,7 +1,9 @@
 import matplotlib.pyplot as plt
 import matplotlib.patches as mpatches
 import numpy as np
+from scipy.stats import kde
 import os
+import math
 
 """This file is part of the Gudhi Library. The Gudhi library
    (Geometric Understanding in Higher Dimensions) is a generic C++
@@ -79,6 +81,8 @@ def plot_persistence_barcode(persistence=[], persistence_file='', alpha=0.6,
     :param inf_delta: Infinity is placed at ((max_death - min_birth) x inf_delta).
         A reasonable value is between 0.05 and 0.5 - default is 0.1.
     :type inf_delta: float.
+    :param legend: Display the dimension color legend (default is False).
+    :type legend: boolean.
     :returns: A matplotlib object containing horizontal bar plot of persistence
         (launch `show()` method on it to display it).
     """
@@ -154,6 +158,8 @@ def plot_persistence_diagram(persistence=[], persistence_file='', alpha=0.6,
     :param inf_delta: Infinity is placed at ((max_death - min_birth) x inf_delta).
         A reasonable value is between 0.05 and 0.5 - default is 0.1.
     :type inf_delta: float.
+    :param legend: Display the dimension color legend (default is False).
+    :type legend: boolean.
     :returns: A matplotlib object containing diagram plot of persistence
         (launch `show()` method on it to display it).
     """
@@ -174,7 +180,6 @@ def plot_persistence_diagram(persistence=[], persistence_file='', alpha=0.6,
         persistence = sorted(persistence, key=lambda life_time: life_time[1][1]-life_time[1][0], reverse=True)[:max_plots]
 
     (min_birth, max_death) = __min_birth_max_death(persistence, band)
-    ind = 0
     delta = ((max_death - min_birth) * inf_delta)
     # Replace infinity values with max_death + delta for diagram to be more
     # readable
@@ -201,7 +206,6 @@ def plot_persistence_diagram(persistence=[], persistence_file='', alpha=0.6,
             # Infinite death case for diagram to be nicer
             plt.scatter(interval[1][0], infinity, alpha=alpha,
                         color = palette[interval[0]])
-        ind = ind + 1
 
     if legend:
         dimensions = list(set(item[0] for item in persistence))
@@ -213,3 +217,69 @@ def plot_persistence_diagram(persistence=[], persistence_file='', alpha=0.6,
     # Ends plot on infinity value and starts a little bit before min_birth
     plt.axis([axis_start, infinity, axis_start, infinity + delta])
     return plt
+
+def plot_persistence_density(persistence=[], persistence_file='', nbins=300,
+                             max_plots=1000, cmap=plt.cm.hot_r, legend=False):
+    """This function plots the persistence density from persistence values list
+    or from a :doc:`persistence file <fileformats>`. Be aware that this
+    function does not distinguish the dimension, it is up to you to select the
+    required one.
+
+    :param persistence: Persistence values list.
+    :type persistence: list of tuples(dimension, tuple(birth, death)).
+    :param persistence_file: A :doc:`persistence file <fileformats>` style name
+        (reset persistence if both are set).
+    :type persistence_file: string
+    :param nbins: Evaluate a gaussian kde on a regular grid of nbins x nbins
+        over data extents (default is 300)
+    :type nbins: int.
+    :param max_plots: number of maximal plots to be displayed
+        Set it to 0 to see all, Default value is 1000.
+        (persistence will be sorted by life time if max_plots is set)
+    :type max_plots: int.
+    :param cmap: A matplotlib colormap (default is matplotlib.pyplot.cm.hot_r).
+    :type cmap: cf. matplotlib colormap.
+    :param legend: Display the color bar values (default is False).
+    :type legend: boolean.
+    :returns: A matplotlib object containing diagram plot of persistence
+        (launch `show()` method on it to display it).
+    """
+    if persistence_file is not '':
+        if os.path.isfile(persistence_file):
+            # Reset persistence
+            persistence = []
+            diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file)
+            for key in diag.keys():
+                for persistence_interval in diag[key]:
+                    persistence.append((key, persistence_interval))
+        else:
+            print("file " + persistence_file + " not found.")
+            return None
+
+    if max_plots > 0 and max_plots < len(persistence):
+        # Sort by life time, then takes only the max_plots elements
+        persistence = sorted(persistence, key=lambda life_time: life_time[1][1]-life_time[1][0], reverse=True)[:max_plots]
+
+    # Set as numpy array birth and death (remove undefined values - inf and NaN)
+    birth = np.asarray([(interval[1][0]) for interval in persistence if (math.isfinite(interval[1][1]) and math.isfinite(interval[1][0]))])
+    death = np.asarray([(interval[1][1]) for interval in persistence if (math.isfinite(interval[1][1]) and math.isfinite(interval[1][0]))])
+
+    # line display of equation : birth = death
+    x = np.linspace(death.min(), birth.max(), 1000)
+    plt.plot(x, x, color='k', linewidth=1.0)
+
+    # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
+    k = kde.gaussian_kde([birth,death])
+    xi, yi = np.mgrid[birth.min():birth.max():nbins*1j, death.min():death.max():nbins*1j]
+    zi = k(np.vstack([xi.flatten(), yi.flatten()]))
+
+    # Make the plot
+    plt.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap)
+
+    if legend:
+        plt.colorbar()
+
+    plt.title('Persistence density')
+    plt.xlabel('Birth')
+    plt.ylabel('Death')
+    return plt
diff --git a/src/cython/doc/persistence_graphical_tools_ref.rst b/src/cython/doc/persistence_graphical_tools_ref.rst
index a2c6bcef..54aff4bc 100644
--- a/src/cython/doc/persistence_graphical_tools_ref.rst
+++ b/src/cython/doc/persistence_graphical_tools_ref.rst
@@ -9,3 +9,4 @@ Persistence graphical tools reference manual
 .. autofunction:: gudhi.__min_birth_max_death
 .. autofunction:: gudhi.plot_persistence_barcode
 .. autofunction:: gudhi.plot_persistence_diagram
+.. autofunction:: gudhi.plot_persistence_density
diff --git a/src/cython/doc/persistence_graphical_tools_user.rst b/src/cython/doc/persistence_graphical_tools_user.rst
index 292915eb..b21de06a 100644
--- a/src/cython/doc/persistence_graphical_tools_user.rst
+++ b/src/cython/doc/persistence_graphical_tools_user.rst
@@ -43,6 +43,9 @@ This function can display the persistence result as a diagram:
         legend=True)
     plt.show()
 
+Persistence density
+-------------------
+
 If you want more information on a specific dimension, for instance:
 
 .. plot::
@@ -56,7 +59,7 @@ If you want more information on a specific dimension, for instance:
         gudhi.read_persistence_intervals_grouped_by_dimension(persistence_file=\
             persistence_file)
     dim = 1
-    # Display all points with some transparency
-    plt = gudhi.plot_persistence_diagram([(dim,interval) for interval in diag[dim]],
-        max_plots=0, alpha=0.1)
+    # Display persistence density
+    plt = gudhi.plot_persistence_density([(dim,interval) for interval in diag[dim]],
+        max_plots=0, legend=True)
     plt.show()
diff --git a/src/cython/setup.py.in b/src/cython/setup.py.in
index ee381a1b..22b2954e 100644
--- a/src/cython/setup.py.in
+++ b/src/cython/setup.py.in
@@ -49,6 +49,7 @@ setup(
     install_requires = [
         "matplotlib",
         "numpy",
+        "scipy",
         "cython",
         ],
 )