From af146a2e48c16855355ac599cbc617250727d244 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 25 Nov 2016 16:00:19 +0000 Subject: 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 --- src/cython/cython/alpha_complex.pyx | 259 ++------------------- src/cython/cython/simplex_tree.pyx | 2 +- src/cython/doc/alpha_complex_user.rst | 39 ++-- .../doc/persistence_graphical_tools_user.rst | 5 +- src/cython/doc/tangential_complex_ref.rst | 10 + src/cython/doc/tangential_complex_sum.rst | 16 ++ src/cython/doc/tangential_complex_user.rst | 144 ++++++++++++ ...ex_diagram_persistence_from_off_file_example.py | 10 +- .../example/alpha_complex_from_points_example.py | 25 +- src/cython/include/Alpha_complex_interface.h | 17 +- src/cython/include/Simplex_tree_interface.h | 7 +- src/cython/test/test_alpha_complex.py | 38 +-- 12 files changed, 266 insertions(+), 306 deletions(-) create mode 100644 src/cython/doc/tangential_complex_ref.rst create mode 100644 src/cython/doc/tangential_complex_sum.rst create mode 100644 src/cython/doc/tangential_complex_user.rst (limited to 'src') 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( 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(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, - 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(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, 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 ≥ 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 #include "Persistent_cohomology_interface.h" +#include "Simplex_tree_interface.h" #include #include // std::pair @@ -49,18 +50,14 @@ class Alpha_complex_interface { typedef typename Simplex_tree<>::Simplex_key Simplex_key; public: - Alpha_complex_interface(std::vector>&points, double max_alpha_square) + Alpha_complex_interface(std::vector>&points) : pcoh_(nullptr) { alpha_complex_ = new Alpha_complex(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(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>* 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 class Simplex_tree_interface : public Simplex_tree { + public: typedef typename Simplex_tree::Simplex_handle Simplex_handle; typedef typename std::pair Insertion_result; using Simplex = std::vector; @@ -50,6 +51,7 @@ class Simplex_tree_interface : public Simplex_tree { bool insert_simplex_and_subfaces(const Simplex& complex, Filtration_value filtration = 0) { Insertion_result result = Simplex_tree::insert_simplex_and_subfaces(complex, filtration); + Simplex_tree::initialize_filtration(); return (result.second); } @@ -58,7 +60,8 @@ class Simplex_tree_interface : public Simplex_tree { } void remove_maximal_simplex(const Simplex& complex) { - return Simplex_tree::remove_maximal_simplex(Simplex_tree::find(complex)); + Simplex_tree::remove_maximal_simplex(Simplex_tree::find(complex)); + Simplex_tree::initialize_filtration(); } Complex_tree get_filtered_tree() { @@ -66,8 +69,10 @@ class Simplex_tree_interface : public Simplex_tree { for (auto f_simplex : Simplex_tree::filtration_simplex_range()) { Simplex simplex; for (auto vertex : Simplex_tree::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::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)] -- cgit v1.2.3