diff options
Diffstat (limited to 'src/python/gudhi')
| -rw-r--r-- | src/python/gudhi/tensorflow/__init__.py | 5 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/cubical_layer.py | 76 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py | 84 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/rips_layer.py | 91 |
4 files changed, 256 insertions, 0 deletions
diff --git a/src/python/gudhi/tensorflow/__init__.py b/src/python/gudhi/tensorflow/__init__.py new file mode 100644 index 00000000..1599cf52 --- /dev/null +++ b/src/python/gudhi/tensorflow/__init__.py @@ -0,0 +1,5 @@ +from .cubical_layer import CubicalLayer +from .lower_star_simplex_tree_layer import LowerStarSimplexTreeLayer +from .rips_layer import RipsLayer + +__all__ = ["LowerStarSimplexTreeLayer", "RipsLayer", "CubicalLayer"] diff --git a/src/python/gudhi/tensorflow/cubical_layer.py b/src/python/gudhi/tensorflow/cubical_layer.py new file mode 100644 index 00000000..369b0e54 --- /dev/null +++ b/src/python/gudhi/tensorflow/cubical_layer.py @@ -0,0 +1,76 @@ +import numpy as np +import tensorflow as tf +from ..cubical_complex import CubicalComplex + +###################### +# Cubical filtration # +###################### + +# The parameters of the model are the pixel values. + +def _Cubical(Xflat, Xdim, dimensions): + # Parameters: Xflat (flattened image), + # Xdim (shape of non-flattened image) + # dimensions (homology dimensions) + + # Compute the persistence pairs with Gudhi + # We reverse the dimensions because CubicalComplex uses Fortran ordering + cc = CubicalComplex(dimensions=Xdim[::-1], top_dimensional_cells=Xflat) + cc.compute_persistence() + + # Retrieve and ouput image indices/pixels corresponding to positive and negative simplices + cof_pp = cc.cofaces_of_persistence_pairs() + + L_cofs = [] + for dim in dimensions: + + try: + cof = cof_pp[0][dim] + except IndexError: + cof = np.array([]) + + L_cofs.append(np.array(cof, dtype=np.int32)) + + return L_cofs + +class CubicalLayer(tf.keras.layers.Layer): + """ + TensorFlow layer for computing the persistent homology of a cubical complex + """ + def __init__(self, dimensions, min_persistence=None, **kwargs): + """ + Constructor for the CubicalLayer class + + Parameters: + dimensions (List[int]): homology dimensions + min_persistence (List[float]): minimum distance-to-diagonal of the points in the output persistence diagrams (default None, in which case 0. is used for all dimensions) + """ + super().__init__(dynamic=True, **kwargs) + self.dimensions = dimensions + self.min_persistence = min_persistence if min_persistence != None else [0.] * len(self.dimensions) + assert len(self.min_persistence) == len(self.dimensions) + + def call(self, X): + """ + Compute persistence diagram associated to a cubical complex filtered by some pixel values + + Parameters: + X (TensorFlow variable): pixel values of the cubical complex + + Returns: + dgms (list of TensorFlow variables): list of cubical persistence diagrams of length self.dimensions, where each element contains a finite persistence diagram of shape [num_finite_points, 2] + """ + # Compute pixels associated to positive and negative simplices + # Don't compute gradient for this operation + Xflat = tf.reshape(X, [-1]) + Xdim = X.shape + indices_list = _Cubical(Xflat.numpy(), Xdim, self.dimensions) + # Get persistence diagram by simply picking the corresponding entries in the image + self.dgms = [tf.reshape(tf.gather(Xflat, indices), [-1,2]) for indices in indices_list] + for idx_dim in range(len(self.min_persistence)): + min_pers = self.min_persistence[idx_dim] + if min_pers >= 0: + finite_dgm = self.dgms[idx_dim] + persistent_indices = tf.where(tf.math.abs(finite_dgm[:,1]-finite_dgm[:,0]) > min_pers) + self.dgms[idx_dim] = tf.gather(finite_dgm, indices=persistent_indices) + return self.dgms diff --git a/src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py b/src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py new file mode 100644 index 00000000..cf7df6de --- /dev/null +++ b/src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py @@ -0,0 +1,84 @@ +import numpy as np +import tensorflow as tf + +######################################### +# Lower star filtration on simplex tree # +######################################### + +# The parameters of the model are the vertex function values of the simplex tree. + +def _LowerStarSimplexTree(simplextree, filtration, dimensions): + # Parameters: simplextree (simplex tree on which to compute persistence) + # filtration (function values on the vertices of st), + # dimensions (homology dimensions), + + simplextree.reset_filtration(-np.inf, 0) + + # Assign new filtration values + for i in range(simplextree.num_vertices()): + simplextree.assign_filtration([i], filtration[i]) + simplextree.make_filtration_non_decreasing() + + # Compute persistence diagram + simplextree.compute_persistence() + + # Get vertex pairs for optimization. First, get all simplex pairs + pairs = simplextree.lower_star_persistence_generators() + + L_indices = [] + for dimension in dimensions: + + finite_pairs = pairs[0][dimension] if len(pairs[0]) >= dimension+1 else np.empty(shape=[0,2]) + essential_pairs = pairs[1][dimension] if len(pairs[1]) >= dimension+1 else np.empty(shape=[0,1]) + + finite_indices = np.array(finite_pairs.flatten(), dtype=np.int32) + essential_indices = np.array(essential_pairs.flatten(), dtype=np.int32) + + L_indices.append((finite_indices, essential_indices)) + + return L_indices + +class LowerStarSimplexTreeLayer(tf.keras.layers.Layer): + """ + TensorFlow layer for computing lower-star persistence out of a simplex tree + """ + def __init__(self, simplextree, dimensions, min_persistence=None, **kwargs): + """ + Constructor for the LowerStarSimplexTreeLayer class + + Parameters: + simplextree (gudhi.SimplexTree): underlying simplex tree. Its vertices MUST be named with integers from 0 to n = number of vertices. Note that its filtration values are modified in each call of the class. + dimensions (List[int]): homology dimensions + min_persistence (List[float]): minimum distance-to-diagonal of the points in the output persistence diagrams (default None, in which case 0. is used for all dimensions) + """ + super().__init__(dynamic=True, **kwargs) + self.dimensions = dimensions + self.simplextree = simplextree + self.min_persistence = min_persistence if min_persistence != None else [0. for _ in range(len(self.dimensions))] + assert len(self.min_persistence) == len(self.dimensions) + + def call(self, filtration): + """ + Compute lower-star persistence diagram associated to a function defined on the vertices of the simplex tree + + Parameters: + F (TensorFlow variable): filter function values over the vertices of the simplex tree. The ith entry of F corresponds to vertex i in self.simplextree + + Returns: + dgms (list of tuple of TensorFlow variables): list of lower-star persistence diagrams of length self.dimensions, where each element of the list is a tuple that contains the finite and essential persistence diagrams of shapes [num_finite_points, 2] and [num_essential_points, 1] respectively + """ + # Don't try to compute gradients for the vertex pairs + indices = _LowerStarSimplexTree(self.simplextree, filtration.numpy(), self.dimensions) + # Get persistence diagrams + self.dgms = [] + for idx_dim, dimension in enumerate(self.dimensions): + finite_dgm = tf.reshape(tf.gather(filtration, indices[idx_dim][0]), [-1,2]) + essential_dgm = tf.reshape(tf.gather(filtration, indices[idx_dim][1]), [-1,1]) + min_pers = self.min_persistence[idx_dim] + if min_pers >= 0: + persistent_indices = tf.where(tf.math.abs(finite_dgm[:,1]-finite_dgm[:,0]) > min_pers) + self.dgms.append((tf.reshape(tf.gather(finite_dgm, indices=persistent_indices),[-1,2]), essential_dgm)) + else: + self.dgms.append((finite_dgm, essential_dgm)) + return self.dgms + diff --git a/src/python/gudhi/tensorflow/rips_layer.py b/src/python/gudhi/tensorflow/rips_layer.py new file mode 100644 index 00000000..7b5edfa3 --- /dev/null +++ b/src/python/gudhi/tensorflow/rips_layer.py @@ -0,0 +1,91 @@ +import numpy as np +import tensorflow as tf +from ..rips_complex import RipsComplex + +############################ +# Vietoris-Rips filtration # +############################ + +# The parameters of the model are the point coordinates. + +def _Rips(DX, max_edge, dimensions): + # Parameters: DX (distance matrix), + # max_edge (maximum edge length for Rips filtration), + # dimensions (homology dimensions) + + # Compute the persistence pairs with Gudhi + rc = RipsComplex(distance_matrix=DX, max_edge_length=max_edge) + st = rc.create_simplex_tree(max_dimension=max(dimensions)+1) + st.compute_persistence() + pairs = st.flag_persistence_generators() + + L_indices = [] + for dimension in dimensions: + + if dimension == 0: + finite_pairs = pairs[0] + essential_pairs = pairs[2] + else: + finite_pairs = pairs[1][dimension-1] if len(pairs[1]) >= dimension else np.empty(shape=[0,4]) + essential_pairs = pairs[3][dimension-1] if len(pairs[3]) >= dimension else np.empty(shape=[0,2]) + + finite_indices = np.array(finite_pairs.flatten(), dtype=np.int32) + essential_indices = np.array(essential_pairs.flatten(), dtype=np.int32) + + L_indices.append((finite_indices, essential_indices)) + + return L_indices + +class RipsLayer(tf.keras.layers.Layer): + """ + TensorFlow layer for computing Rips persistence out of a point cloud + """ + def __init__(self, dimensions, maximum_edge_length=np.inf, min_persistence=None, **kwargs): + """ + Constructor for the RipsLayer class + + Parameters: + maximum_edge_length (float): maximum edge length for the Rips complex + dimensions (List[int]): homology dimensions + min_persistence (List[float]): minimum distance-to-diagonal of the points in the output persistence diagrams (default None, in which case 0. is used for all dimensions) + """ + super().__init__(dynamic=True, **kwargs) + self.max_edge = maximum_edge_length + self.dimensions = dimensions + self.min_persistence = min_persistence if min_persistence != None else [0. for _ in range(len(self.dimensions))] + assert len(self.min_persistence) == len(self.dimensions) + + def call(self, X): + """ + Compute Rips persistence diagram associated to a point cloud + + Parameters: + X (TensorFlow variable): point cloud of shape [number of points, number of dimensions] + + Returns: + dgms (list of tuple of TensorFlow variables): list of Rips persistence diagrams of length self.dimensions, where each element of the list is a tuple that contains the finite and essential persistence diagrams of shapes [num_finite_points, 2] and [num_essential_points, 1] respectively + """ + # Compute distance matrix + DX = tf.norm(tf.expand_dims(X, 1)-tf.expand_dims(X, 0), axis=2) + # Compute vertices associated to positive and negative simplices + # Don't compute gradient for this operation + indices = _Rips(DX.numpy(), self.max_edge, self.dimensions) + # Get persistence diagrams by simply picking the corresponding entries in the distance matrix + self.dgms = [] + for idx_dim, dimension in enumerate(self.dimensions): + cur_idx = indices[idx_dim] + if dimension > 0: + finite_dgm = tf.reshape(tf.gather_nd(DX, tf.reshape(cur_idx[0], [-1,2])), [-1,2]) + essential_dgm = tf.reshape(tf.gather_nd(DX, tf.reshape(cur_idx[1], [-1,2])), [-1,1]) + else: + reshaped_cur_idx = tf.reshape(cur_idx[0], [-1,3]) + finite_dgm = tf.concat([tf.zeros([reshaped_cur_idx.shape[0],1]), tf.reshape(tf.gather_nd(DX, reshaped_cur_idx[:,1:]), [-1,1])], axis=1) + essential_dgm = tf.zeros([cur_idx[1].shape[0],1]) + min_pers = self.min_persistence[idx_dim] + if min_pers >= 0: + persistent_indices = tf.where(tf.math.abs(finite_dgm[:,1]-finite_dgm[:,0]) > min_pers) + self.dgms.append((tf.reshape(tf.gather(finite_dgm, indices=persistent_indices),[-1,2]), essential_dgm)) + else: + self.dgms.append((finite_dgm, essential_dgm)) + return self.dgms + |