Add gudhi backend

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- (#33) Add gudhi backend.
+
 ### Fixed
 
 - (#30) Update `.readthedocs.yaml` and `docs/conf.py` to current RTDs requirements.

cechmate/__init__.py

@@ -1,7 +1,7 @@
-"""Cechmate init file."""
-
 # from .filtrations import *
 # from .solver import *
 # from .utils import *
 
-from ._version import __version__ as __version__
+from ._version import __version__
+
+__all__ = ["__version__"]

cechmate/filtrations/__init__.py

@@ -1,5 +1,7 @@
-from .alpha import Alpha as Alpha
-from .rips import Rips as Rips
-from .cech import Cech as Cech
-from .extended import Extended as Extended
+from .alpha import Alpha
+from .cech import Cech
+from .rips import Rips
+# from .extended import *
 # from .miniball import get_boundary
+
+__all__ = ["Alpha", "Cech", "Rips"]

cechmate/filtrations/alpha.py

@@ -1,15 +1,21 @@
+from collections import defaultdict
 import itertools
 import time
 import warnings
+import logging
 
 import numpy as np
 from numpy import linalg
 from scipy import spatial
 
-from .base import BaseFiltration
+from cechmate.filtrations.base import BaseFiltration
+import gudhi
+from numpy.typing import NDArray
 
 __all__ = ["Alpha"]
 
+logger = logging.getLogger(__name__)
+
 
 class Alpha(BaseFiltration):
     """Construct an Alpha filtration from the given data.
@@ -39,7 +45,10 @@ def build(self, X):
         X: Nxd array
             Array of N Euclidean vectors in d dimensions
         """
-
+        warnings.warn(
+            "This function is deprecated and will be removed in a future release. Use fit instead.",
+            DeprecationWarning,
+        )
         if X.shape[0] < X.shape[1]:
             warnings.warn(
                 "The input point cloud has more columns than rows; "
@@ -61,6 +70,10 @@ def build(self, X):
                 "Finished spatial.Delaunay triangulation (Elapsed Time %.3g)"
                 % (time.time() - tic)
             )
+            logger.info(
+                "Finished spatial.Delaunay triangulation (Elapsed Time %.3g)"
+                % (time.time() - tic)
+            )
             print("Building alpha filtration...")
             tic = time.time()
 
@@ -110,13 +123,132 @@ def build(self, X):
                 % (time.time() - tic)
             )
 
+        logger.info(
+            "Finished building alpha filtration (Elapsed Time %.3g)"
+            % (time.time() - tic)
+        )
+
+        # NOTE: The following line makes this method have special return type. The 0-simplices
+        # are indexed on lists, whereas all other simplices are indexed on tuples.
         simplices = [([i], 0) for i in range(X.shape[0])]
         simplices.extend(filtration.items())
 
         self.simplices_ = simplices
 
         return simplices
 
+    def fit(self, X) -> list[tuple[tuple[np.int32], np.float64]]:
+        """
+        Do the Alpha filtration of a Euclidean point set (requires scipy)
+
+        Parameters
+        ===========
+        X: Nxd array
+            Array of N Euclidean vectors in d dimensions
+        """
+
+        if X.shape[0] < X.shape[1]:
+            warnings.warn(
+                "The input point cloud has more dimensions than data points; "
+                + "did you mean to transpose?"
+            )
+        maxdim = self.maxdim or X.shape[1] - 1
+
+        ## Step 1: Figure out the filtration
+        if self.verbose:
+            print("Doing spatial.Delaunay triangulation...")
+            tic = time.time()
+
+        logger.info("Doing spatial.Delaunay triangulation...")
+        delaunay_faces = spatial.Delaunay(X).simplices
+
+        if self.verbose:
+            print(
+                "Finished spatial.Delaunay triangulation (Elapsed Time %.3g)"
+                % (time.time() - tic)
+            )
+            logger.info(
+                "Finished spatial.Delaunay triangulation (Elapsed Time %.3g)"
+                % (time.time() - tic)
+            )
+            print("Building alpha filtration...")
+            tic = time.time()
+
+        logger.info("Building alpha filtration...")
+
+        filtration = defaultdict(lambda: float("inf"))
+        circumcenter_cache = {}
+
+        filtration = {(i,): np.float64(0.0) for i in range(X.shape[0])}
+
+        for dim in range(maxdim + 2, 1, -1):
+            for s in range(delaunay_faces.shape[0]):
+                simplex = delaunay_faces[s, :]
+                for sigma in itertools.combinations(simplex, dim):
+                    sigma = tuple(sorted(sigma))
+
+                    if sigma not in filtration:
+                        if sigma not in circumcenter_cache:
+                            circumcenter_cache[sigma] = self._get_circumcenter(
+                                X[sigma, :]
+                            )
+                        _, rSqr = circumcenter_cache[sigma]
+                        if np.isfinite(rSqr):
+                            filtration[sigma] = rSqr
+
+                    if sigma in filtration:
+                        for i in range(dim):  # Propagate alpha filtration value
+                            tau = sigma[0:i] + sigma[i + 1 : :]
+                            if tau not in filtration:
+                                if len(tau) > 1:
+                                    if tau not in circumcenter_cache:
+                                        circumcenter_cache[tau] = (
+                                            self._get_circumcenter(X[tau, :])
+                                        )
+                                    xtau, rtauSqr = circumcenter_cache[tau]
+                                    if np.sum((X[sigma[i], :] - xtau) ** 2) < rtauSqr:
+                                        filtration[tau] = filtration[sigma]
+                            else:
+                                filtration[tau] = min(
+                                    filtration[tau], filtration[sigma]
+                                )
+
+        # Convert from squared radii to radii
+        for sigma in filtration:
+            filtration[sigma] = np.sqrt(filtration[sigma])
+
+        ## Step 2: Take care of numerical artifacts that may result
+        ## in simplices with greater filtration values than their co-faces
+        simplices_bydim = [set([]) for _ in range(maxdim + 2)]
+        for simplex in filtration.keys():
+            simplices_bydim[len(simplex) - 1].add(simplex)
+
+        simplices_bydim = simplices_bydim[2::]
+        simplices_bydim.reverse()
+
+        for simplices_dim in simplices_bydim:
+            for sigma in simplices_dim:
+                for i in range(len(sigma)):
+                    tau = sigma[0:i] + sigma[i + 1 : :]
+                    if filtration[tau] > filtration[sigma]:
+                        filtration[tau] = filtration[sigma]
+
+        if self.verbose:
+            print(
+                "Finished building alpha filtration (Elapsed Time %.3g)"
+                % (time.time() - tic)
+            )
+        logger.info(
+            "Finished building alpha filtration (Elapsed Time %.3g)"
+            % (time.time() - tic)
+        )
+
+        simplices = list(filtration.items())
+
+        self.simplices_ = simplices
+
+        return simplices
+
     def _get_circumcenter(self, X):
         """
         Compute the circumcenter and circumradius of a simplex
@@ -189,3 +321,30 @@ def minor(A, j):
                 x = x.dot(V.T) + muV
             return (x, rSqr)
         return (np.inf, np.inf)  # SC2 (Points not in general position)
+
+    def transform(self, simplices=None, ripser_format=True) -> list[NDArray]:
+        """
+        Compute persistent homology.
+        """
+        simplices_ = simplices or self.simplices_
+
+        simplex_tree = gudhi.SimplexTree()
+        for simplex, filtration_value in simplices_:
+            simplex_tree.insert(simplex, filtration_value)
+
+        persistence = simplex_tree.persistence()
+
+        if not ripser_format:
+            return persistence
+
+        # convert to ripser.py format
+        ripser_output = []
+        for dim, (birth, death) in persistence:
+            while len(ripser_output) <= dim:
+                ripser_output.append([])
+            if death == float("inf"):
+                death = -1
+            ripser_output[dim].append(np.array([birth, death]))
+        ripser_output = [np.array(dgm) for dgm in ripser_output]
+
+        return ripser_output

cechmate/filtrations/base.py

@@ -1,13 +1,25 @@
 """All filtrations should have a base interface."""
 
+from typing import Literal
+import warnings
+from numpy.typing import NDArray
+import numpy as np
+
 # from ..solver import phat_diagrams
 
 
 class BaseFiltration:
     """Base filtration that implements constructor and `diagrams` method."""
 
-    def __init__(self, maxdim=None, verbose=True):
-        """Default constructor
+    def __init__(
+        self,
+        maxdim=None,
+        verbose=True,
+        solver: Literal["phat", "gudhi", "ripser"] = "gudhi",
+    ):
+        """Base filtration class.
+
+        Not all solvers are implemented for all filtrations. See each filtration class for details.
 
         Parameters
         ----------
@@ -16,11 +28,13 @@
             Maximum dimension of homology to compute
         verbose: boolean
             If True, then print logging statements.
-
+        solver: Literal["phat", "gudhi", "ripser"], default="gudhi"
+            Solver to use for persistent homology.
         """
 
         self.maxdim = maxdim
         self.verbose = verbose
+        self.solver = solver
 
         self.simplices_ = None
         self.diagrams_ = None
@@ -42,7 +56,44 @@
             the persistence diagram for Hk
 
         """
+        warnings.warn(
+            "This function is deprecated and will be removed in a future release. Use transform instead.",
+            DeprecationWarning,
+        )
         simplices = simplices or self.simplices_
-        # self.diagrams_ = phat_diagrams(simplices, show_inf)
+        # TODO: Update this call.
+        self.diagrams_ = phat_diagrams(simplices, show_inf)
 
         return self.diagrams_
+
+    def transform(self, simplices=None, ripser_format=True) -> list[NDArray]:
+        """
+        Compute persistent homology.
+        """
+        import gudhi
+
+        simplices_ = simplices or self.simplices_
+
+        if simplices_ is None:
+            raise ValueError("No simplices to transform.")
+
+        simplex_tree = gudhi.SimplexTree()
+        for simplex, filtration_value in simplices_:
+            simplex_tree.insert(simplex, filtration_value)
+
+        persistence = simplex_tree.persistence()
+
+        if not ripser_format:
+            return persistence
+
+        # convert to ripser.py format
+        ripser_output = []
+        for dim, (birth, death) in persistence:
+            while len(ripser_output) <= dim:
+                ripser_output.append([])
+            if death == float("inf"):
+                death = -1
+            ripser_output[dim].append(np.array([birth, death]))
+        ripser_output = [np.array(dgm) for dgm in ripser_output]
+
+        return ripser_output

cechmate/filtrations/cech.py

@@ -1,4 +1,5 @@
 import itertools
+import warnings
 import numpy as np
 
 from .base import BaseFiltration
@@ -35,7 +36,10 @@ def build(self, X):
         simplices:
             Cech filtration for the data X
         """
-
+        warnings.warn(
+            "This function is deprecated and will be removed in a future release. Use fit instead.",
+            DeprecationWarning,
+        )
         N = X.shape[0]
         xr = np.arange(N)
         xrl = xr.tolist()
@@ -57,3 +61,38 @@ def build(self, X):
         self.simplices_ = simplices
 
         return simplices
+
+    def fit(self, X) -> list[tuple[list[int], int]]:
+        """Compute the Cech filtration of a Euclidean point set for simplices up to order :code:`self.max_dim`.
+
+        Parameters
+        ===========
+
+        X: Nxd array
+            N Euclidean vectors in d dimensions
+
+        Returns
+        ==========
+
+        simplices:
+            Cech filtration for the data X
+        """
+
+        N = X.shape[0]
+        xr = np.arange(N)
+        xrl = xr.tolist()
+        maxdim = self.maxdim or X.shape[1] - 1
+        miniball = miniball_cache(X)
+
+        # start with vertices
+        simplices = [([i], 0) for i in range(N)]
+
+        # then higher order simplices
+        for k in range(maxdim + 1):
+            for idxs in itertools.combinations(xrl, k + 2):
+                _, r2 = miniball(frozenset(idxs), frozenset([]))
+                simplices.append((list(idxs), np.sqrt(r2)))
+
+        self.simplices_ = simplices
+
+        return simplices