Add of tangential complex doc

Separate simplex tree from alpha complex



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

Former-commit-id: 1cf4a35e0a099501eb1cb6b9809041dd1a1e2617

12 files changed, 266 insertions, 306 deletions

diff --git a/src/cython/cython/alpha_complex.pyx b/src/cython/cython/alpha_complex.pyx
index 85771300..9bf7b257 100644
--- a/src/cython/cython/alpha_complex.pyx
+++ b/src/cython/cython/alpha_complex.pyx
@@ -3,6 +3,7 @@ from libcpp.vector cimport vector
 from libcpp.utility cimport pair
 from libcpp.string cimport string
 from libcpp cimport bool
+from cython.operator cimport dereference as deref
 import os
 
 """This file is part of the Gudhi Library. The Gudhi library
@@ -33,9 +34,9 @@ __license__ = "GPL v3"
 
 cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
     cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface":
-        Alpha_complex_interface(vector[vector[double]] points, double max_alpha_square)
+        Alpha_complex_interface(vector[vector[double]] points)
         # bool from_file is a workaround for cython to find the correct signature
-        Alpha_complex_interface(string off_file, double max_alpha_square, bool from_file)
+        Alpha_complex_interface(string off_file, bool from_file)
         double filtration()
         double simplex_filtration(vector[int] simplex)
         void set_filtration(double filtration)
@@ -57,6 +58,7 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
         vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence)
         vector[int] get_betti_numbers()
         vector[int] get_persistent_betti_numbers(double from_value, double to_value)
+        void create_simplex_tree(Simplex_tree_interface_full_featured simplex_tree, double max_alpha_square)
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
@@ -81,7 +83,7 @@ cdef class AlphaComplex:
     cdef Alpha_complex_interface * thisptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=None, off_file='', max_alpha_square=float('inf')):
+    def __init__(self, points=None, off_file=''):
         """AlphaComplex constructor.
 
         :param points: A list of points in d-Dimension.
@@ -91,22 +93,17 @@ cdef class AlphaComplex:
 
         :param off_file: An OFF file style name.
         :type off_file: string
-
-        :param max_alpha_square: Maximum Alpha square value. Default is :math:`\infty`
-        :type max_alpha_square: double
         """
 
     # The real cython constructor
-    def __cinit__(self, points=[], off_file='', max_alpha_square=float('inf')):
+    def __cinit__(self, points=[], off_file=''):
         if off_file is not '':
             if os.path.isfile(off_file):
-                self.thisptr = new Alpha_complex_interface(off_file,
-                    max_alpha_square, True)
+                self.thisptr = new Alpha_complex_interface(off_file, True)
             else:
                 print("file " + off_file + " not found.")
         else:
-            self.thisptr = new Alpha_complex_interface(points,
-                max_alpha_square)
+            self.thisptr = new Alpha_complex_interface(points)
                 
 
     def __dealloc__(self):
@@ -118,195 +115,6 @@ cdef class AlphaComplex:
          """
         return self.thisptr != NULL
 
-    def get_filtration(self):
-        """This function returns the main simplicial complex filtration value.
-
-        :returns:  float -- the simplicial complex filtration value.
-        """
-        return self.thisptr.filtration()
-
-    def filtration(self, simplex):
-        """This function returns the simplicial complex filtration value for a
-        given N-simplex.
-
-        :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int
-        :returns:  float -- the simplicial complex filtration value.
-        """
-        return self.thisptr.simplex_filtration(simplex)
-
-    def set_filtration(self, filtration):
-        """This function sets the main simplicial complex filtration value.
-
-        :param filtration: The filtration value.
-        :type filtration: float.
-        """
-        self.thisptr.set_filtration(<double> filtration)
-
-    def initialize_filtration(self):
-        """This function initializes and sorts the simplicial complex
-        filtration vector.
-
-        .. note::
-
-            This function must be launched before persistence, betti_numbers,
-            persistent_betti_numbers or get_filtered_tree after inserting or
-            removing simplices.
-        """
-        self.thisptr.initialize_filtration()
-
-    def num_vertices(self):
-        """This function returns the number of vertices of the simplicial
-        complex.
-
-        :returns:  int -- the simplicial complex number of vertices.
-        """
-        return self.thisptr.num_vertices()
-
-    def num_simplices(self):
-        """This function returns the number of simplices of the simplicial
-        complex.
-
-        :returns:  int -- the simplicial complex number of simplices.
-        """
-        return self.thisptr.num_simplices()
-
-    def dimension(self):
-        """This function returns the dimension of the simplicial complex.
-
-        :returns:  int -- the simplicial complex dimension.
-        """
-        return self.thisptr.dimension()
-
-    def set_dimension(self, dimension):
-        """This function sets the dimension of the simplicial complex.
-
-        :param dimension: The new dimension value.
-        :type dimension: int.
-        """
-        self.thisptr.set_dimension(<int>dimension)
-
-    def find(self, simplex):
-        """This function returns if the N-simplex was found in the simplicial
-        complex or not.
-
-        :param simplex: The N-simplex to find, represented by a list of vertex.
-        :type simplex: list of int.
-        :returns:  bool -- true if the simplex was found, false otherwise.
-        """
-        cdef vector[int] complex
-        for i in simplex:
-            complex.push_back(i)
-        return self.thisptr.find_simplex(complex)
-
-    def insert(self, simplex, filtration=0.0):
-        """This function inserts the given N-simplex with the given filtration
-        value (default value is '0.0').
-
-        :param simplex: The N-simplex to insert, represented by a list of
-            vertex.
-        :type simplex: list of int.
-        :param filtration: The filtration value of the simplex.
-        :type filtration: float.
-        :returns:  bool -- true if the simplex was found, false otherwise.
-        """
-        cdef vector[int] complex
-        for i in simplex:
-            complex.push_back(i)
-        return self.thisptr.insert_simplex_and_subfaces(complex,
-                                                        <double>filtration)
-
-    def get_filtered_tree(self):
-        """This function returns the tree sorted by increasing filtration
-        values.
-
-        :returns:  list of tuples(simplex, filtration) -- the tree sorted by
-            increasing filtration values.
-        """
-        cdef vector[pair[vector[int], double]] coface_tree \
-            = self.thisptr.get_filtered_tree()
-        ct = []
-        for filtered_complex in coface_tree:
-            v = []
-            for vertex in filtered_complex.first:
-                v.append(vertex)
-            ct.append((v, filtered_complex.second))
-        return ct
-
-    def get_skeleton_tree(self, dimension):
-        """This function returns the tree skeleton of a maximum given
-        dimension.
-
-        :param dimension: The skeleton dimension value.
-        :type dimension: int.
-        :returns:  list of tuples(simplex, filtration) -- the skeleton tree
-            of a maximum dimension.
-        """
-        cdef vector[pair[vector[int], double]] coface_tree \
-            = self.thisptr.get_skeleton_tree(<int>dimension)
-        ct = []
-        for filtered_complex in coface_tree:
-            v = []
-            for vertex in filtered_complex.first:
-                v.append(vertex)
-            ct.append((v, filtered_complex.second))
-        return ct
-
-    def get_star_tree(self, simplex):
-        """This function returns the star tree of a given N-simplex.
-
-        :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
-        :returns:  list of tuples(simplex, filtration) -- the star tree of a
-            simplex.
-        """
-        cdef vector[int] complex
-        for i in simplex:
-            complex.push_back(i)
-        cdef vector[pair[vector[int], double]] coface_tree \
-            = self.thisptr.get_star_tree(complex)
-        ct = []
-        for filtered_complex in coface_tree:
-            v = []
-            for vertex in filtered_complex.first:
-                v.append(vertex)
-            ct.append((v, filtered_complex.second))
-        return ct
-
-    def get_coface_tree(self, simplex, codimension):
-        """This function returns the coface tree of a given N-simplex with a
-        given codimension.
-
-        :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
-        :param codimension: The codimension. If codimension = 0, all cofaces
-            are returned (equivalent of get_star_tree function)
-        :type codimension: int.
-        :returns: list of tuples(simplex, filtration) -- the coface tree of a
-            simplex.
-        """
-        cdef vector[int] complex
-        for i in simplex:
-            complex.push_back(i)
-        cdef vector[pair[vector[int], double]] coface_tree \
-            = self.thisptr.get_coface_tree(complex, <int>codimension)
-        ct = []
-        for filtered_complex in coface_tree:
-            v = []
-            for vertex in filtered_complex.first:
-                v.append(vertex)
-            ct.append((v, filtered_complex.second))
-        return ct
-
-    def remove_maximal_simplex(self, simplex):
-        """This function removes a given maximal N-simplex from the simplicial
-        complex.
-
-        :param simplex: The N-simplex, represented by a list of vertex.
-        :type simplex: list of int.
-        """
-        self.thisptr.remove_maximal_simplex(simplex)
-
     def get_point(self, vertex):
         """This function returns the point corresponding to a given vertex.
 
@@ -317,47 +125,14 @@ cdef class AlphaComplex:
         cdef vector[double] point = self.thisptr.get_point(vertex)
         return point
 
-    def persistence(self, homology_coeff_field=11, min_persistence=0.0):
-        """This function returns the persistence of the simplicial complex.
-
-        :param homology_coeff_field: The homology coefficient field. Must be a
-            prime number
-        :type homology_coeff_field: int.
-        :param min_persistence: The minimum persistence value to take into
-            account (strictly greater than min_persistence). Default value is
-            0.0.
-            Sets min_persistence to -1.0 to see all values.
-        :type min_persistence: float.
-        :note: list of pairs(dimension, pair(birth, death)) -- the
-            persistence of the simplicial complex.
-        """
-        return self.thisptr.get_persistence(homology_coeff_field, min_persistence)
-
-    def betti_numbers(self):
-        """This function returns the Betti numbers of the simplicial complex.
-
-        :returns: list of int -- The Betti numbers ([B0, B1, ..., Bn]).
-
-        :note: betti_numbers function requires persistence function to be
-            launched first.
-        """
-        return self.thisptr.get_betti_numbers()
-
-    def persistent_betti_numbers(self, from_value, to_value):
-        """This function returns the persistent Betti numbers of the
-        simplicial complex.
-
-        :param from_value: The persistence birth limit to be added in the
-            numbers (persistent birth <= from_value).
-        :type from_value: float.
-        :param to_value: The persistence death limit to be added in the
-            numbers (persistent death > to_value).
-        :type to_value: float.
-
-        :returns: list of int -- The persistent Betti numbers ([B0, B1, ...,
-            Bn]).
+    def create_simplex_tree(self, SimplexTree simplex_tree, max_alpha_square=float('inf')):
+        """This function creates the given simplex tree from the Delaunay
+        Triangulation.
 
-        :note: persistent_betti_numbers function requires persistence
-            function to be launched first.
+        :param simplex_tree: The simplex tree to create (must be empty)
+        :type simplex tree: gudhi.SimplexTree.
+        :param max_alpha_square: The maximum alpha square threshold the
+            simplices shall not exceed. Default is set to infinity.
+        :type max_alpha_square: float.
         """
-        return self.thisptr.get_persistent_betti_numbers(from_value, to_value)
+        self.thisptr.create_simplex_tree(deref(simplex_tree.thisptr), max_alpha_square)
\ No newline at end of file
diff --git a/src/cython/cython/simplex_tree.pyx b/src/cython/cython/simplex_tree.pyx
index bf91e812..cf7cbe26 100644
--- a/src/cython/cython/simplex_tree.pyx
+++ b/src/cython/cython/simplex_tree.pyx
@@ -75,7 +75,7 @@ cdef class SimplexTree:
     cdef Simplex_tree_persistence_interface * pcohptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=[], max_alpha_square=float('inf')):
+    def __init__(self):
         """SimplexTree constructor.
         """
 
diff --git a/src/cython/doc/alpha_complex_user.rst b/src/cython/doc/alpha_complex_user.rst
index b6b6e550..5ad3a9e7 100644
--- a/src/cython/doc/alpha_complex_user.rst
+++ b/src/cython/doc/alpha_complex_user.rst
@@ -22,15 +22,17 @@ This example builds the Delaunay triangulation from the given points, and initia
 
 .. testcode::
 
-   import gudhi
-   alpha_complex = gudhi.AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]],
-                                      max_alpha_square=60.0)
-   result_str = 'Alpha complex is of dimension ' + repr(alpha_complex.dimension()) + ' - ' + \
-       repr(alpha_complex.num_simplices()) + ' simplices - ' + \
-       repr(alpha_complex.num_vertices()) + ' vertices.'
-   print(result_str)
-   for fitered_value in alpha_complex.get_filtered_tree():
-       print(fitered_value)
+    import gudhi
+    alpha_complex = gudhi.AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]])
+
+    simplex_tree = gudhi.SimplexTree()
+    alpha_complex.create_simplex_tree(simplex_tree, max_alpha_square=60.0)
+    result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \
+        repr(simplex_tree.num_simplices()) + ' simplices - ' + \
+        repr(simplex_tree.num_vertices()) + ' vertices.'
+    print(result_str)
+    for fitered_value in simplex_tree.get_filtered_tree():
+        print(fitered_value)
 
 The output is:
 
@@ -152,15 +154,16 @@ Then, it is asked to display information about the alpha complex:
 
 .. testcode::
 
-   import gudhi
-   alpha_complex = gudhi.AlphaComplex(off_file='alphacomplexdoc.off',
-                                      max_alpha_square=59.0)
-   result_str = 'Alpha complex is of dimension ' + repr(alpha_complex.dimension()) + ' - ' + \
-       repr(alpha_complex.num_simplices()) + ' simplices - ' + \
-       repr(alpha_complex.num_vertices()) + ' vertices.'
-   print(result_str)
-   for fitered_value in alpha_complex.get_filtered_tree():
-       print(fitered_value)
+    import gudhi
+    alpha_complex = gudhi.AlphaComplex(off_file='alphacomplexdoc.off')
+    simplex_tree = gudhi.SimplexTree()
+    alpha_complex.create_simplex_tree(simplex_tree, max_alpha_square=59.0)
+    result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \
+        repr(simplex_tree.num_simplices()) + ' simplices - ' + \
+        repr(simplex_tree.num_vertices()) + ' vertices.'
+    print(result_str)
+    for fitered_value in simplex_tree.get_filtered_tree():
+        print(fitered_value)
 
 the program output is:
 
diff --git a/src/cython/doc/persistence_graphical_tools_user.rst b/src/cython/doc/persistence_graphical_tools_user.rst
index 0cb61cf1..b23ad5e6 100644
--- a/src/cython/doc/persistence_graphical_tools_user.rst
+++ b/src/cython/doc/persistence_graphical_tools_user.rst
@@ -41,6 +41,7 @@ This function can display the persistence result as a diagram:
     import gudhi
     
     alpha_complex = gudhi.AlphaComplex(off_file='tore3D_300.off')
-    alpha_complex.initialize_filtration()
-    diag = alpha_complex.persistence()
+    simplex_tree = gudhi.SimplexTree()
+    alpha_complex.create_simplex_tree(simplex_tree)
+    diag = simplex_tree.persistence()
     gudhi.diagram_persistence(diag)
diff --git a/src/cython/doc/tangential_complex_ref.rst b/src/cython/doc/tangential_complex_ref.rst
new file mode 100644
index 00000000..35589475
--- /dev/null
+++ b/src/cython/doc/tangential_complex_ref.rst
@@ -0,0 +1,10 @@
+===================================
+Tangential complex reference manual
+===================================
+
+.. autoclass:: gudhi.TangentialComplex
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+   .. automethod:: gudhi.TangentialComplex.__init__
diff --git a/src/cython/doc/tangential_complex_sum.rst b/src/cython/doc/tangential_complex_sum.rst
new file mode 100644
index 00000000..a2c0be0e
--- /dev/null
+++ b/src/cython/doc/tangential_complex_sum.rst
@@ -0,0 +1,16 @@
+=====================================  =====================================  =====================================
+:Author: Clément Jamin                 :Introduced in: GUDHI 1.4.0            :Copyright: GPL v3
+=====================================  =====================================  =====================================
+:Requires: CGAL &ge; 4.8.0             Eigen3
+=====================================  =====================================  =====================================
+
++-------------------------------------------+----------------------------------------------------------------------+
+| .. image::                                | A Tangential Delaunay complex is a simplicial complex designed to    |
+|      img/tc_examples.png                  | reconstruct a :math:`k`-dimensional manifold embedded in :math:`d`-  |
+|                                           | dimensional Euclidean space. The input is a point sample coming from |
+|                                           | an unknown manifold. The running time depends only linearly on the   |
+|                                           | extrinsic dimension :math:`d` and exponentially on the intrinsic     |
+|                                           | dimension :math:`k`.                                                 |
++-------------------------------------------+----------------------------------------------------------------------+
+| :doc:`tangential_complex_user`            | :doc:`tangential_complex_ref`                                             |
++-------------------------------------------+----------------------------------------------------------------------+
diff --git a/src/cython/doc/tangential_complex_user.rst b/src/cython/doc/tangential_complex_user.rst
new file mode 100644
index 00000000..588de08c
--- /dev/null
+++ b/src/cython/doc/tangential_complex_user.rst
@@ -0,0 +1,144 @@
+==============================
+Tangential complex user manual
+==============================
+.. include:: tangential_complex_sum.rst
+
+Definition
+----------
+
+A Tangential Delaunay complex is a simplicial complex designed to reconstruct a
+:math:`k`-dimensional smooth manifold embedded in :math:`d`-dimensional
+Euclidean space. The input is a point sample coming from an unknown manifold,
+which means that the points lie close to a structure of "small" intrinsic
+dimension. The running time depends only linearly on the extrinsic dimension
+:math:`d` and exponentially on the intrinsic dimension :math:`k`.
+
+An extensive description of the Tangential complex can be found in
+:cite:`tangentialcomplex2014`).
+
+What is a Tangential Complex?
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Let us start with the description of the Tangential complex of a simple
+example, with :math:`k = 1` and :math:`d = 2`. The input data is 4 points
+:math:`P` located on a curve embedded in 2D.
+
+.. figure:: img/tc_example_01.png
+    :alt: The input
+    :figclass: align-center
+    The input
+
+For each point :math:`p`, estimate its tangent subspace :math:`T_p` (e.g.
+using PCA). 
+
+.. figure:: img/tc_example_02.png
+    :alt: The estimated normals
+    :figclass: align-center
+    The estimated normals
+
+Let us add the Voronoi diagram of the points in orange. For each point
+:math:`p`, construct its star in the Delaunay triangulation of :math:`P`
+restricted to :math:`T_p`.
+
+.. figure:: img/tc_example_03.png
+    :alt: The Voronoi diagram
+    :figclass: align-center
+    The Voronoi diagram
+
+The Tangential Delaunay complex is the union of those stars.
+
+In practice, neither the ambient Voronoi diagram nor the ambient Delaunay
+triangulation is computed. Instead, local :math:`k`-dimensional regular
+triangulations are computed with a limited number of points as we only need the
+star of each point. More details can be found in :cite:`tangentialcomplex2014`.
+
+Inconsistencies
+^^^^^^^^^^^^^^^
+Inconsistencies between the stars can occur. An inconsistency occurs when a
+simplex is not in the star of all its vertices.
+
+Let us take the same example.
+
+.. figure:: img/tc_example_07_before.png
+    :alt: Before
+    :figclass: align-center
+    Before
+
+Let us slightly move the tangent subspace :math:`T_q`
+
+.. figure:: img/tc_example_07_after.png
+    :alt: After
+    :figclass: align-center
+    After
+
+Now, the star of :math:`Q` contains :math:`QP`, but the star of :math:`P` does
+not contain :math:`QP`. We have an inconsistency.
+
+.. figure:: img/tc_example_08.png
+    :alt: After
+    :figclass: align-center
+    After
+
+One way to solve inconsistencies is to randomly perturb the positions of the
+points involved in an inconsistency. In the current implementation, this
+perturbation is done in the tangent subspace of each point. The maximum
+perturbation radius is given as a parameter to the constructor.
+
+In most cases, we recommend to provide a point set where the minimum distance
+between any two points is not too small. This can be achieved using the
+functions provided by the Subsampling module. Then, a good value to start with
+for the maximum perturbation radius would be around half the minimum distance
+between any two points. The Example with perturbation below shows an example of
+such a process.
+
+In most cases, this process is able to dramatically reduce the number of
+inconsistencies, but is not guaranteed to succeed.
+
+Output
+^^^^^^
+The result of the computation is exported as a Simplex_tree. It is the union of
+the stars of all the input points. A vertex in the Simplex Tree is the index of
+the point in the range provided by the user. The point corresponding to a
+vertex can also be obtained through the Tangential_complex::get_point function.
+Note that even if the positions of the points are perturbed, their original
+positions are kept (e.g. Tangential_complex::get_point returns the original
+position of the point).
+
+The result can be obtained after the computation of the Tangential complex
+itself and/or after the perturbation process.
+
+
+Simple example
+--------------
+
+This example builds the Tangential complex of point set. Note that the
+dimension of the kernel here is dynamic, which is slower, but more flexible:
+the intrinsic and ambient dimensions does not have to be known at compile-time.
+
+testcode::
+
+   import gudhi
+   tc = gudhi.TangentialComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]])
+
+The output is:
+
+testoutput::
+
+   Tangential complex is of dimension 2 - 25 simplices - 7 vertices.
+
+
+Example with perturbation
+-------------------------
+
+This example builds the Tangential complex of a point set, then tries to solve
+inconsistencies by perturbing the positions of points involved in inconsistent
+simplices.
+
+testcode::
+
+   import gudhi
+   tc = gudhi.TangentialComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]])
+
+The output is:
+
+testoutput::
diff --git a/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py
index d74eba57..40d06f93 100755
--- a/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py
+++ b/src/cython/example/alpha_complex_diagram_persistence_from_off_file_example.py
@@ -52,15 +52,17 @@ with open(args.file, 'r') as f:
         message = "AlphaComplex with max_edge_length=" + repr(args.max_alpha_square)
         print(message)
         
-        alpha_complex = gudhi.AlphaComplex(off_file=args.file, max_alpha_square=args.max_alpha_square)
+        alpha_complex = gudhi.AlphaComplex(off_file=args.file)
+        simplex_tree = gudhi.SimplexTree()
+        alpha_complex.create_simplex_tree(simplex_tree, max_alpha_square=args.max_alpha_square)
     
-        message = "Number of simplices=" + repr(alpha_complex.num_simplices())
+        message = "Number of simplices=" + repr(simplex_tree.num_simplices())
         print(message)
         
-        diag = alpha_complex.persistence()
+        diag = simplex_tree.persistence()
     
         print("betti_numbers()=")
-        print(alpha_complex.betti_numbers())
+        print(simplex_tree.betti_numbers())
     
         if args.no_diagram == False:
             gudhi.diagram_persistence(diag)
diff --git a/src/cython/example/alpha_complex_from_points_example.py b/src/cython/example/alpha_complex_from_points_example.py
index da98f3d4..45319f7f 100755
--- a/src/cython/example/alpha_complex_from_points_example.py
+++ b/src/cython/example/alpha_complex_from_points_example.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
-import gudhi
+from gudhi import AlphaComplex, SimplexTree
 
 """This file is part of the Gudhi Library. The Gudhi library
    (Geometric Understanding in Higher Dimensions) is a generic C++
@@ -30,38 +30,39 @@ __license__ = "GPL v3"
 
 print("#####################################################################")
 print("AlphaComplex creation from points")
-alpha_complex = gudhi.AlphaComplex(points=[[0, 0], [1, 0], [0, 1], [1, 1]],
-                                   max_alpha_square=60.0)
+alpha_complex = AlphaComplex(points=[[0, 0], [1, 0], [0, 1], [1, 1]])
+simplex_tree = SimplexTree()
+alpha_complex.create_simplex_tree(simplex_tree, max_alpha_square=60.0)
 
-if alpha_complex.find([0, 1]):
+if simplex_tree.find([0, 1]):
     print("[0, 1] Found !!")
 else:
     print("[0, 1] Not found...")
 
-if alpha_complex.find([4]):
+if simplex_tree.find([4]):
     print("[4] Found !!")
 else:
     print("[4] Not found...")
 
-if alpha_complex.insert([0, 1, 2], filtration=4.0):
+if simplex_tree.insert([0, 1, 2], filtration=4.0):
     print("[0, 1, 2] Inserted !!")
 else:
     print("[0, 1, 2] Not inserted...")
 
-if alpha_complex.insert([0, 1, 4], filtration=4.0):
+if simplex_tree.insert([0, 1, 4], filtration=4.0):
     print("[0, 1, 4] Inserted !!")
 else:
     print("[0, 1, 4] Not inserted...")
 
-if alpha_complex.find([4]):
+if simplex_tree.find([4]):
     print("[4] Found !!")
 else:
     print("[4] Not found...")
 
-print("dimension=", alpha_complex.dimension())
-print("filtered_tree=", alpha_complex.get_filtered_tree())
-print("star([0])=", alpha_complex.get_star_tree([0]))
-print("coface([0], 1)=", alpha_complex.get_coface_tree([0], 1))
+print("dimension=", simplex_tree.dimension())
+print("filtered_tree=", simplex_tree.get_filtered_tree())
+print("star([0])=", simplex_tree.get_star_tree([0]))
+print("coface([0], 1)=", simplex_tree.get_coface_tree([0], 1))
 
 print("point[0]=", alpha_complex.get_point(0))
 print("point[5]=", alpha_complex.get_point(5))
diff --git a/src/cython/include/Alpha_complex_interface.h b/src/cython/include/Alpha_complex_interface.h
index 13412994..07f5b5b4 100644
--- a/src/cython/include/Alpha_complex_interface.h
+++ b/src/cython/include/Alpha_complex_interface.h
@@ -28,6 +28,7 @@
 #include <CGAL/Epick_d.h>
 
 #include "Persistent_cohomology_interface.h"
+#include "Simplex_tree_interface.h"
 
 #include <vector>
 #include <utility>  // std::pair
@@ -49,18 +50,14 @@ class Alpha_complex_interface {
   typedef typename Simplex_tree<>::Simplex_key Simplex_key;
 
  public:
-  Alpha_complex_interface(std::vector<std::vector<double>>&points, double max_alpha_square)
+  Alpha_complex_interface(std::vector<std::vector<double>>&points)
   : pcoh_(nullptr) {
     alpha_complex_ = new Alpha_complex<Dynamic_kernel>(points);
-    alpha_complex_->create_complex(simplex_tree_, max_alpha_square);
-    simplex_tree_.initialize_filtration();
   }
 
-  Alpha_complex_interface(std::string off_file_name, double max_alpha_square, bool from_file = true)
+  Alpha_complex_interface(std::string off_file_name, bool from_file = true)
   : pcoh_(nullptr) {
     alpha_complex_ = new Alpha_complex<Dynamic_kernel>(off_file_name);
-    alpha_complex_->create_complex(simplex_tree_, max_alpha_square);
-    simplex_tree_.initialize_filtration();
   }
 
   bool find_simplex(const Simplex& vh) {
@@ -194,6 +191,11 @@ class Alpha_complex_interface {
     return persistent_betti_numbers;
   }
   
+  void create_simplex_tree(Simplex_tree_interface<>& simplex_tree, double max_alpha_square) {
+    alpha_complex_->create_complex(simplex_tree, max_alpha_square);
+    simplex_tree.initialize_filtration();
+  }
+
  private:
   Simplex_tree<> simplex_tree_;
   Persistent_cohomology_interface<Simplex_tree<>>* pcoh_;
@@ -202,7 +204,6 @@ class Alpha_complex_interface {
 
 }  // namespace alpha_complex
 
-} // namespace Gudhi
+}  // namespace Gudhi
 
 #endif  // ALPHA_COMPLEX_INTERFACE_H
-
diff --git a/src/cython/include/Simplex_tree_interface.h b/src/cython/include/Simplex_tree_interface.h
index ad7be1a8..06ec6d40 100644
--- a/src/cython/include/Simplex_tree_interface.h
+++ b/src/cython/include/Simplex_tree_interface.h
@@ -36,6 +36,7 @@ namespace Gudhi {
 
 template<typename SimplexTreeOptions = Simplex_tree_options_full_featured>
 class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
+ public:
   typedef typename Simplex_tree<SimplexTreeOptions>::Simplex_handle Simplex_handle;
   typedef typename std::pair<Simplex_handle, bool> Insertion_result;
   using Simplex = std::vector<Vertex_handle>;
@@ -50,6 +51,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
 
   bool insert_simplex_and_subfaces(const Simplex& complex, Filtration_value filtration = 0) {
     Insertion_result result = Simplex_tree<SimplexTreeOptions>::insert_simplex_and_subfaces(complex, filtration);
+    Simplex_tree<SimplexTreeOptions>::initialize_filtration();
     return (result.second);
   }
 
@@ -58,7 +60,8 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
   }
 
   void remove_maximal_simplex(const Simplex& complex) {
-    return Simplex_tree<SimplexTreeOptions>::remove_maximal_simplex(Simplex_tree<SimplexTreeOptions>::find(complex));
+    Simplex_tree<SimplexTreeOptions>::remove_maximal_simplex(Simplex_tree<SimplexTreeOptions>::find(complex));
+    Simplex_tree<SimplexTreeOptions>::initialize_filtration();
   }
 
   Complex_tree get_filtered_tree() {
@@ -66,8 +69,10 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     for (auto f_simplex : Simplex_tree<SimplexTreeOptions>::filtration_simplex_range()) {
       Simplex simplex;
       for (auto vertex : Simplex_tree<SimplexTreeOptions>::simplex_vertex_range(f_simplex)) {
+        std::cout << " " << vertex;
         simplex.insert(simplex.begin(), vertex);
       }
+      std::cout << std::endl;
       filtered_tree.push_back(std::make_pair(simplex, Simplex_tree<SimplexTreeOptions>::filtration(f_simplex)));
     }
     return filtered_tree;
diff --git a/src/cython/test/test_alpha_complex.py b/src/cython/test/test_alpha_complex.py
index 13c6f7e8..3731ccda 100755
--- a/src/cython/test/test_alpha_complex.py
+++ b/src/cython/test/test_alpha_complex.py
@@ -1,4 +1,4 @@
-from gudhi import AlphaComplex
+from gudhi import AlphaComplex, SimplexTree
 
 """This file is part of the Gudhi Library. The Gudhi library
    (Geometric Understanding in Higher Dimensions) is a generic C++
@@ -28,30 +28,30 @@ __license__ = "GPL v3"
 
 
 def test_empty_alpha():
-    alpha_complex = AlphaComplex()
+    alpha_complex = AlphaComplex(points=[[0,0]])
     assert alpha_complex.__is_defined() == True
-    assert alpha_complex.__is_persistence_defined() == False
-    assert alpha_complex.num_simplices() == 0
-    assert alpha_complex.num_vertices() == 0
 
 def test_infinite_alpha():
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
     alpha_complex = AlphaComplex(points=point_list)
     assert alpha_complex.__is_defined() == True
-    assert alpha_complex.__is_persistence_defined() == False
 
-    assert alpha_complex.num_simplices() == 11
-    assert alpha_complex.num_vertices() == 4
+    simplex_tree = SimplexTree()
+    alpha_complex.create_simplex_tree(simplex_tree)
+    assert simplex_tree.__is_persistence_defined() == False
+
+    assert simplex_tree.num_simplices() == 11
+    assert simplex_tree.num_vertices() == 4
  
-    assert alpha_complex.get_filtered_tree() == \
+    assert simplex_tree.get_filtered_tree() == \
            [([0], 0.0), ([1], 0.0), ([2], 0.0), ([3], 0.0),
             ([0, 1], 0.25), ([0, 2], 0.25), ([1, 3], 0.25),
             ([2, 3], 0.25), ([1, 2], 0.5), ([0, 1, 2], 0.5),
             ([1, 2, 3], 0.5)]
-    assert alpha_complex.get_star_tree([0]) == \
+    assert simplex_tree.get_star_tree([0]) == \
            [([0], 0.0), ([0, 1], 0.25), ([0, 1, 2], 0.5),
            ([0, 2], 0.25)]
-    assert alpha_complex.get_coface_tree([0], 1) == \
+    assert simplex_tree.get_coface_tree([0], 1) == \
            [([0, 1], 0.25), ([0, 2], 0.25)]
  
     assert point_list[0] == alpha_complex.get_point(0)
@@ -63,11 +63,13 @@ def test_infinite_alpha():
 
 def test_filtered_alpha():
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    filtered_alpha = AlphaComplex(points=point_list,
-                                  max_alpha_square=0.25)
+    filtered_alpha = AlphaComplex(points=point_list)
+
+    simplex_tree = SimplexTree()
+    filtered_alpha.create_simplex_tree(simplex_tree, max_alpha_square=0.25)
 
-    assert filtered_alpha.num_simplices() == 8
-    assert filtered_alpha.num_vertices() == 4
+    assert simplex_tree.num_simplices() == 8
+    assert simplex_tree.num_vertices() == 4
 
     assert point_list[0] == filtered_alpha.get_point(0)
     assert point_list[1] == filtered_alpha.get_point(1)
@@ -76,11 +78,11 @@ def test_filtered_alpha():
     assert filtered_alpha.get_point(4) == []
     assert filtered_alpha.get_point(125) == []
 
-    assert filtered_alpha.get_filtered_tree() == \
+    assert simplex_tree.get_filtered_tree() == \
            [([0], 0.0), ([1], 0.0), ([2], 0.0), ([3], 0.0),
             ([0, 1], 0.25), ([0, 2], 0.25), ([1, 3], 0.25),
             ([2, 3], 0.25)]
-    assert filtered_alpha.get_star_tree([0]) == \
+    assert simplex_tree.get_star_tree([0]) == \
            [([0], 0.0), ([0, 1], 0.25), ([0, 2], 0.25)]
-    assert filtered_alpha.get_coface_tree([0], 1) == \
+    assert simplex_tree.get_coface_tree([0], 1) == \
            [([0, 1], 0.25), ([0, 2], 0.25)]