diff options
Diffstat (limited to 'src/python')
| -rw-r--r-- | src/python/CMakeLists.txt | 4 | ||||
| -rw-r--r-- | src/python/gudhi/differentiation/__init__.py | 3 | ||||
| -rw-r--r-- | src/python/gudhi/differentiation/tensorflow.py | 234 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/CubicalLayer.py | 66 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/LowerStarSimplexTreeLayer.py | 77 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/RipsLayer.py | 75 | ||||
| -rw-r--r-- | src/python/gudhi/tensorflow/__init__.py | 5 | ||||
| -rw-r--r-- | src/python/test/test_diff.py | 2 |
8 files changed, 226 insertions, 240 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 6c17223e..49ec3fea 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -67,7 +67,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_strong_witness_complex', ") # Modules that should not be auto-imported in __init__.py set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'representations', ") - set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'differentiation', ") + set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'tensorflow', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'point_cloud', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'weighted_rips_complex', ") @@ -271,7 +271,7 @@ if(PYTHONINTERP_FOUND) file(COPY "gudhi/persistence_graphical_tools.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/") file(COPY "gudhi/wasserstein" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") - file(COPY "gudhi/differentiation" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") + file(COPY "gudhi/tensorflow" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/clustering" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi" FILES_MATCHING PATTERN "*.py") file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") diff --git a/src/python/gudhi/differentiation/__init__.py b/src/python/gudhi/differentiation/__init__.py deleted file mode 100644 index 3b7790e4..00000000 --- a/src/python/gudhi/differentiation/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from .tensorflow import * - -__all__ = ["LowerStarSimplexTreeLayer", "RipsLayer", "CubicalLayer"] diff --git a/src/python/gudhi/differentiation/tensorflow.py b/src/python/gudhi/differentiation/tensorflow.py deleted file mode 100644 index 15d5811e..00000000 --- a/src/python/gudhi/differentiation/tensorflow.py +++ /dev/null @@ -1,234 +0,0 @@ -import numpy as np -import tensorflow as tf -from ..rips_complex import RipsComplex -from ..cubical_complex import CubicalComplex - -# In this file, we write functions based on the Gudhi library that compute persistence diagrams associated to -# different filtrations (lower star, Rips, cubical), as well as the corresponding positive and negative -# simplices. We also wrap these functions into Tensorflow models. - - - -######################################### -# Lower star filtration on simplex tree # -######################################### - -# The parameters of the model are the vertex function values of the simplex tree. - -def _LowerStarSimplexTree(simplextree, filtration, dimension): - # Parameters: simplextree (simplex tree on which to compute persistence) - # filtration (function values on the vertices of st), - # dimension (homology dimension), - - for s,_ in simplextree.get_filtration(): - simplextree.assign_filtration(s, -1e10) - - # 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 - dgm = simplextree.persistence() - - # Get vertex pairs for optimization. First, get all simplex pairs - pairs = simplextree.persistence_pairs() - - # Then, loop over all simplex pairs - indices, pers = [], [] - for s1, s2 in pairs: - # Select pairs with good homological dimension and finite lifetime - if len(s1) == dimension+1 and len(s2) > 0: - # Get IDs of the vertices corresponding to the filtration values of the simplices - l1, l2 = np.array(s1), np.array(s2) - i1 = l1[np.argmax(filtration[l1])] - i2 = l2[np.argmax(filtration[l2])] - indices.append(i1) - indices.append(i2) - # Compute lifetime - pers.append(simplextree.filtration(s2)-simplextree.filtration(s1)) - - # Sort vertex pairs wrt lifetime - perm = np.argsort(pers) - indices = np.reshape(indices, [-1,2])[perm][::-1,:].flatten() - - return np.array(indices, dtype=np.int32) - -class LowerStarSimplexTreeLayer(tf.keras.layers.Layer): - """ - TensorFlow layer for computing lower-star persistence out of a simplex tree - - Attributes: - simplextree (gudhi.SimplexTree()): underlying simplex tree - dimension (int): homology dimension - """ - def __init__(self, simplextree, dimension=0, **kwargs): - super().__init__(dynamic=True, **kwargs) - self.dimension = dimension - self.simplextree = simplextree - - def build(self): - super.build() - - 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 - """ - # Don't try to compute gradients for the vertex pairs - indices = tf.stop_gradient(_LowerStarSimplexTree(self.simplextree, filtration.numpy(), self.dimension)) - # Get persistence diagram - self.dgm = tf.reshape(tf.gather(filtration, indices), [-1,2]) - return self.dgm - - - - - - - - - -############################ -# Vietoris-Rips filtration # -############################ - -# The parameters of the model are the point coordinates. - -def _Rips(DX, max_edge, dimension): - # Parameters: DX (distance matrix), - # max_edge (maximum edge length for Rips filtration), - # dimension (homology dimension) - - # Compute the persistence pairs with Gudhi - rc = RipsComplex(distance_matrix=DX, max_edge_length=max_edge) - st = rc.create_simplex_tree(max_dimension=dimension+1) - dgm = st.persistence() - pairs = st.persistence_pairs() - - # Retrieve vertices v_a and v_b by picking the ones achieving the maximal - # distance among all pairwise distances between the simplex vertices - indices, pers = [], [] - for s1, s2 in pairs: - if len(s1) == dimension+1 and len(s2) > 0: - l1, l2 = np.array(s1), np.array(s2) - i1 = [l1[v] for v in np.unravel_index(np.argmax(DX[l1,:][:,l1]),[len(l1), len(l1)])] - i2 = [l2[v] for v in np.unravel_index(np.argmax(DX[l2,:][:,l2]),[len(l2), len(l2)])] - indices.append(i1) - indices.append(i2) - pers.append(st.filtration(s2)-st.filtration(s1)) - - # Sort points with distance-to-diagonal - perm = np.argsort(pers) - indices = np.reshape(indices, [-1,4])[perm][::-1,:].flatten() - - return np.array(indices, dtype=np.int32) - -class RipsLayer(tf.keras.layers.Layer): - """ - TensorFlow layer for computing Rips persistence out of a point cloud - - Attributes: - maximum_edge_length (float): maximum edge length for the Rips complex - dimension (int): homology dimension - """ - def __init__(self, maximum_edge_length=12, dimension=1, **kwargs): - super().__init__(dynamic=True, **kwargs) - self.max_edge = maximum_edge_length - self.dimension = dimension - - def build(self): - super.build() - - 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] - """ - # Compute distance matrix - DX = tf.math.sqrt(tf.reduce_sum((tf.expand_dims(X, 1)-tf.expand_dims(X, 0))**2, 2)) - # Compute vertices associated to positive and negative simplices - # Don't compute gradient for this operation - indices = tf.stop_gradient(_Rips(DX.numpy(), self.max_edge, self.dimension)) - # Get persistence diagram by simply picking the corresponding entries in the distance matrix - if self.dimension > 0: - dgm = tf.reshape(tf.gather_nd(DX, tf.reshape(indices, [-1,2])), [-1,2]) - else: - indices = tf.reshape(indices, [-1,2])[1::2,:] - dgm = tf.concat([tf.zeros([indices.shape[0],1]), tf.reshape(tf.gather_nd(DX, indices), [-1,1])], axis=1) - return dgm - - - - - - - - - -###################### -# Cubical filtration # -###################### - -# The parameters of the model are the pixel values. - -def _Cubical(X, dimension): - # Parameters: X (image), - # dimension (homology dimension) - - # Compute the persistence pairs with Gudhi - cc = CubicalComplex(dimensions=X.shape, top_dimensional_cells=X.flatten()) - cc.persistence() - try: - cof = cc.cofaces_of_persistence_pairs()[0][dimension] - except IndexError: - cof = np.array([]) - - if len(cof) > 0: - # Sort points with distance-to-diagonal - Xs = X.shape - pers = [X[np.unravel_index(cof[idx,1], Xs)] - X[np.unravel_index(cof[idx,0], Xs)] for idx in range(len(cof))] - perm = np.argsort(pers) - cof = cof[perm[::-1]] - - # Retrieve and ouput image indices/pixels corresponding to positive and negative simplices - D = len(Xs) if len(cof) > 0 else 1 - ocof = np.array([0 for _ in range(D*2*cof.shape[0])]) - count = 0 - for idx in range(0,2*cof.shape[0],2): - ocof[D*idx:D*(idx+1)] = np.unravel_index(cof[count,0], Xs) - ocof[D*(idx+1):D*(idx+2)] = np.unravel_index(cof[count,1], Xs) - count += 1 - return np.array(ocof, dtype=np.int32) - -class CubicalLayer(tf.keras.layers.Layer): - """ - TensorFlow layer for computing cubical persistence out of a cubical complex - - Attributes: - dimension (int): homology dimension - """ - def __init__(self, dimension=1, **kwargs): - super().__init__(dynamic=True, **kwargs) - self.dimension = dimension - - def build(self): - super.build() - - 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 - """ - # Compute pixels associated to positive and negative simplices - # Don't compute gradient for this operation - indices = tf.stop_gradient(_Cubical(X.numpy(), self.dimension)) - # Get persistence diagram by simply picking the corresponding entries in the image - dgm = tf.reshape(tf.gather_nd(X, tf.reshape(indices, [-1,len(X.shape)])), [-1,2]) - return dgm diff --git a/src/python/gudhi/tensorflow/CubicalLayer.py b/src/python/gudhi/tensorflow/CubicalLayer.py new file mode 100644 index 00000000..e36adec5 --- /dev/null +++ b/src/python/gudhi/tensorflow/CubicalLayer.py @@ -0,0 +1,66 @@ +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(X, dimension): + # Parameters: X (image), + # dimension (homology dimension) + + # Compute the persistence pairs with Gudhi + cc = CubicalComplex(dimensions=X.shape, top_dimensional_cells=X.flatten()) + cc.persistence() + try: + cof = cc.cofaces_of_persistence_pairs()[0][dimension] + except IndexError: + cof = np.array([]) + + if len(cof) > 0: + # Sort points with distance-to-diagonal + Xs = X.shape + pers = [X[np.unravel_index(cof[idx,1], Xs)] - X[np.unravel_index(cof[idx,0], Xs)] for idx in range(len(cof))] + perm = np.argsort(pers) + cof = cof[perm[::-1]] + + # Retrieve and ouput image indices/pixels corresponding to positive and negative simplices + D = len(Xs) if len(cof) > 0 else 1 + ocof = np.array([0 for _ in range(D*2*cof.shape[0])]) + count = 0 + for idx in range(0,2*cof.shape[0],2): + ocof[D*idx:D*(idx+1)] = np.unravel_index(cof[count,0], Xs) + ocof[D*(idx+1):D*(idx+2)] = np.unravel_index(cof[count,1], Xs) + count += 1 + return np.array(ocof, dtype=np.int32) + +class CubicalLayer(tf.keras.layers.Layer): + """ + TensorFlow layer for computing cubical persistence out of a cubical complex + + Attributes: + dimension (int): homology dimension + """ + def __init__(self, dimension=1, **kwargs): + super().__init__(dynamic=True, **kwargs) + self.dimension = dimension + + def build(self): + super.build() + + 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 + """ + # Compute pixels associated to positive and negative simplices + # Don't compute gradient for this operation + indices = tf.stop_gradient(_Cubical(X.numpy(), self.dimension)) + # Get persistence diagram by simply picking the corresponding entries in the image + dgm = tf.reshape(tf.gather_nd(X, tf.reshape(indices, [-1,len(X.shape)])), [-1,2]) + return dgm diff --git a/src/python/gudhi/tensorflow/LowerStarSimplexTreeLayer.py b/src/python/gudhi/tensorflow/LowerStarSimplexTreeLayer.py new file mode 100644 index 00000000..fc963d2f --- /dev/null +++ b/src/python/gudhi/tensorflow/LowerStarSimplexTreeLayer.py @@ -0,0 +1,77 @@ +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, dimension): + # Parameters: simplextree (simplex tree on which to compute persistence) + # filtration (function values on the vertices of st), + # dimension (homology dimension), + + for s,_ in simplextree.get_filtration(): + simplextree.assign_filtration(s, -1e10) + + # 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 + dgm = simplextree.persistence() + + # Get vertex pairs for optimization. First, get all simplex pairs + pairs = simplextree.persistence_pairs() + + # Then, loop over all simplex pairs + indices, pers = [], [] + for s1, s2 in pairs: + # Select pairs with good homological dimension and finite lifetime + if len(s1) == dimension+1 and len(s2) > 0: + # Get IDs of the vertices corresponding to the filtration values of the simplices + l1, l2 = np.array(s1), np.array(s2) + i1 = l1[np.argmax(filtration[l1])] + i2 = l2[np.argmax(filtration[l2])] + indices.append(i1) + indices.append(i2) + # Compute lifetime + pers.append(simplextree.filtration(s2)-simplextree.filtration(s1)) + + # Sort vertex pairs wrt lifetime + perm = np.argsort(pers) + indices = np.reshape(indices, [-1,2])[perm][::-1,:].flatten() + + return np.array(indices, dtype=np.int32) + +class LowerStarSimplexTreeLayer(tf.keras.layers.Layer): + """ + TensorFlow layer for computing lower-star persistence out of a simplex tree + + Attributes: + simplextree (gudhi.SimplexTree()): underlying simplex tree + dimension (int): homology dimension + """ + def __init__(self, simplextree, dimension=0, **kwargs): + super().__init__(dynamic=True, **kwargs) + self.dimension = dimension + self.simplextree = simplextree + + def build(self): + super.build() + + 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 + """ + # Don't try to compute gradients for the vertex pairs + indices = tf.stop_gradient(_LowerStarSimplexTree(self.simplextree, filtration.numpy(), self.dimension)) + # Get persistence diagram + self.dgm = tf.reshape(tf.gather(filtration, indices), [-1,2]) + return self.dgm + diff --git a/src/python/gudhi/tensorflow/RipsLayer.py b/src/python/gudhi/tensorflow/RipsLayer.py new file mode 100644 index 00000000..373e021e --- /dev/null +++ b/src/python/gudhi/tensorflow/RipsLayer.py @@ -0,0 +1,75 @@ +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, dimension): + # Parameters: DX (distance matrix), + # max_edge (maximum edge length for Rips filtration), + # dimension (homology dimension) + + # Compute the persistence pairs with Gudhi + rc = RipsComplex(distance_matrix=DX, max_edge_length=max_edge) + st = rc.create_simplex_tree(max_dimension=dimension+1) + dgm = st.persistence() + pairs = st.persistence_pairs() + + # Retrieve vertices v_a and v_b by picking the ones achieving the maximal + # distance among all pairwise distances between the simplex vertices + indices, pers = [], [] + for s1, s2 in pairs: + if len(s1) == dimension+1 and len(s2) > 0: + l1, l2 = np.array(s1), np.array(s2) + i1 = [l1[v] for v in np.unravel_index(np.argmax(DX[l1,:][:,l1]),[len(l1), len(l1)])] + i2 = [l2[v] for v in np.unravel_index(np.argmax(DX[l2,:][:,l2]),[len(l2), len(l2)])] + indices.append(i1) + indices.append(i2) + pers.append(st.filtration(s2)-st.filtration(s1)) + + # Sort points with distance-to-diagonal + perm = np.argsort(pers) + indices = np.reshape(indices, [-1,4])[perm][::-1,:].flatten() + + return np.array(indices, dtype=np.int32) + +class RipsLayer(tf.keras.layers.Layer): + """ + TensorFlow layer for computing Rips persistence out of a point cloud + + Attributes: + maximum_edge_length (float): maximum edge length for the Rips complex + dimension (int): homology dimension + """ + def __init__(self, maximum_edge_length=12, dimension=1, **kwargs): + super().__init__(dynamic=True, **kwargs) + self.max_edge = maximum_edge_length + self.dimension = dimension + + def build(self): + super.build() + + 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] + """ + # Compute distance matrix + DX = tf.math.sqrt(tf.reduce_sum((tf.expand_dims(X, 1)-tf.expand_dims(X, 0))**2, 2)) + # Compute vertices associated to positive and negative simplices + # Don't compute gradient for this operation + indices = tf.stop_gradient(_Rips(DX.numpy(), self.max_edge, self.dimension)) + # Get persistence diagram by simply picking the corresponding entries in the distance matrix + if self.dimension > 0: + dgm = tf.reshape(tf.gather_nd(DX, tf.reshape(indices, [-1,2])), [-1,2]) + else: + indices = tf.reshape(indices, [-1,2])[1::2,:] + dgm = tf.concat([tf.zeros([indices.shape[0],1]), tf.reshape(tf.gather_nd(DX, indices), [-1,1])], axis=1) + return dgm + diff --git a/src/python/gudhi/tensorflow/__init__.py b/src/python/gudhi/tensorflow/__init__.py new file mode 100644 index 00000000..47335a25 --- /dev/null +++ b/src/python/gudhi/tensorflow/__init__.py @@ -0,0 +1,5 @@ +from .CubicalLayer import CubicalLayer +from .LowerStarSimplexTreeLayer import LowerStarSimplexTreeLayer +from .RipsLayer import RipsLayer + +__all__ = ["LowerStarSimplexTreeLayer", "RipsLayer", "CubicalLayer"] diff --git a/src/python/test/test_diff.py b/src/python/test/test_diff.py index 129b9f03..73a03697 100644 --- a/src/python/test/test_diff.py +++ b/src/python/test/test_diff.py @@ -1,4 +1,4 @@ -from gudhi.differentiation import * +from gudhi.tensorflow import * import numpy as np import tensorflow as tf import gudhi as gd |