up

.github/workflows/sphinx.yml

@@ -0,0 +1,18 @@
+name: Deploy Sphinx documentation to Pages
+
+on:
+  push:
+    branches: [main] # branch to trigger deployment
+
+jobs:
+  pages:
+    runs-on: ubuntu-20.04
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+    permissions:
+      pages: write
+      id-token: write
+    steps:
+      - id: deployment
+        uses: sphinx-notes/pages@v3

pyproject.toml

@@ -3,11 +3,12 @@ requires = ["flit_core >=2,<4"]
 build-backend = "flit_core.buildapi"
 
 [project]
-name =  "py-project-toml"
+name =  "totalpersistence"
 authors = [
-    {name = "Daniel Ciborowski", email = "[email protected]"},
+    {name = "Pablo Riera", email = "[email protected]"},
+    {name = "Manuela Cerdeiro", email = "[email protected]"}
 ]
-description = "Sample Python Project for creating a new Python Module"
+description = "Total persistence for persistence diagrams"
 readme = "README.md"
 classifiers = [
     "Development Status :: 3 - Alpha",
@@ -21,11 +22,16 @@ classifiers = [
 ]
 requires-python = ">=3.8.1"
 dynamic = ["version"]
-
-[project.optional-dependencies]
-spark = [
-    "pyspark>=3.0.0"
+dependencies = [
+    "numpy>=1.21.0",
+    "scipy>=1.7.0",
+    "matplotlib>=3.4.0",
+    "pandas>=1.3.0",
+    "ripser>=0.6.0",
+    "persim>=0.3.0"
 ]
+
+
 test = [
     "bandit[toml]==1.7.5",
     "black==23.3.0",
@@ -47,12 +53,12 @@ test = [
 ]
 
 [project.urls]
-Documentation = "https://github.com/microsoft/python-package-template/tree/main#readme"
-Source = "https://github.com/microsoft/python-package-template"
-Tracker = "https://github.com/microsoft/python-package-template/issues"
+Documentation = "https://github.com/habla-liaa/totalpersistence"
+Source = "https://github.com/habla-liaa/totalpersistence"
+Tracker = "https://github.com/habla-liaa/totalpersistence/issues"
 
 [tool.flit.module]
-name = "python_package"
+name = "totalpersistence"
 
 [tool.bandit]
 exclude_dirs = ["build","dist","tests","scripts"]

src/totalpersistence/totalpersistence.py

@@ -0,0 +1,66 @@
+import numpy as np
+from ripser import ripser
+from scipy.spatial.distance import squareform
+from .utils import (
+    lipschitz,
+    general_position_distance_matrix,
+    to_condensed_form,
+    conematrix,
+    kercoker_bars,
+    matrix_size_from_condensed,
+)
+from IPython import embed
+
+
+def totalpersistence(X, Y, f, maxdim=1, eps=0, tol=1e-11, perturb=1e-7):
+    dX = general_position_distance_matrix(X, perturb)
+    dY = general_position_distance_matrix(Y, perturb)
+
+    data, dgm, dgmX, dgmY = kercoker_via_cone(dX, dY, f, maxdim, eps, tol, perturb)
+    return data, dgm, dgmX, dgmY
+
+
+def kercoker_via_cone(dX, dY, f, maxdim=1, cone_eps=0, tol=1e-11):
+    """
+    TODO: Compute the total persistence diagram using the cone algorithm.
+
+    Parameters
+    ----------
+    dX : np.array
+        Distance matrix of the source space in condensed form.
+    dY : np.array
+        Distance matrix of the target space in condensed form.
+    f : np.array
+        Function values.
+    """
+
+    n = matrix_size_from_condensed(dX)
+    m = matrix_size_from_condensed(dY)
+
+    f = np.array(f)
+
+    # dY_ff = d(f(x_i),f(x_j)) para todo i,j
+    i, j = np.triu_indices(n, k=1)
+    f_i, f_j = f[i], f[j]
+    f_pos = to_condensed_form(f_i, f_j, m)
+    dY_ff = dY[f_pos.astype(int)]
+
+    # dY_fy = d(f(x_i),y_j) para todo i,j
+    indices = np.indices((n, m))
+    i = indices[0].flatten()
+    j = indices[1].flatten()
+    f_i = f[i]
+    DY_fy = np.zeros((n, m))
+    DY_fy[i, j] = squareform(dY)[f_i, j]
+
+    L = lipschitz(dX, dY_ff)
+    dY = dY / L
+
+    D = conematrix(squareform(dX), squareform(dY), DY_fy, cone_eps)
+
+    dgmX = ripser(squareform(dX), distance_matrix=True, maxdim=maxdim)["dgms"]
+    dgmY = ripser(squareform(dY), distance_matrix=True, maxdim=maxdim)["dgms"]
+    dgm = ripser(D, maxdim=maxdim, distance_matrix=True)["dgms"]
+
+    bars = kercoker_bars(dgm, dgmX, dgmY, cone_eps, tol)
+    return bars, dgm, dgmX, dgmY

src/totalpersistence/utils.py

@@ -0,0 +1,94 @@
+import numpy as np
+from scipy.spatial.distance import pdist, squareform
+from IPython import embed
+
+
+def findclose(x, A, tol=1e-5):
+    return ((x + tol) >= A) & ((x - tol) <= A)
+
+
+def lipschitz(dX, dY):
+    return np.max(dY / dX)
+
+
+def matrix_size_from_condensed(dX):
+    n = len(dX)
+    return int(0.5 * (np.sqrt(8 * n + 1) - 1) + 1)
+
+
+def to_condensed_form(i, j, m):
+    return m * i + j - ((i + 2) * (i + 1)) // 2.0
+
+
+def general_position_distance_matrix(X, perturb=1e-7):
+    n = len(X)
+    Xperturbation = perturb * np.random.rand((n * (n - 1) // 2))
+    dX = pdist(X) + Xperturbation
+    return dX
+
+
+def conematrix(DX, DY, DY_fy, eps):
+    n = len(DX)
+    m = len(DY)
+
+    D = np.zeros((n + m + 1, n + m + 1))
+    D[0:n, 0:n] = DX
+    D[n : n + m, n : n + m] = DY
+
+    D[0:n, n : n + m] = DY_fy
+    D[n : n + m, 0:n] = DY_fy.T
+
+    R = max(DX.max(), DY_fy.max()) + 1
+
+    D[n + m, n : n + m] = R
+    D[n : n + m, n + m] = R
+
+    D[n + m, :n] = eps
+    D[:n, n + m] = eps
+
+    return D
+
+
+def kercoker_bars(dgm, dgmX, dgmY, cone_eps, tol=1e-11):
+    """
+    Find cokernel and kernel bars in the persistence diagram.
+    TODO: optimize
+    """
+    data = []
+    for k in range(len(dgm)):
+        coker = []
+        ker = []
+        for r in dgm[k]:
+            b, d = r
+            if d > cone_eps + tol:
+                # coker
+                # b_c = b_y_i
+                # d_c = d_y_i
+                m = findclose(b, dgmY[k][:, 0], tol) & findclose(d, dgmY[k][:, 1], tol)
+                if sum(m):
+                    coker.append((b, d))
+
+                # b_c = b_y_i
+                # d_c = b_x_j
+                if any(findclose(b, dgmY[k][:, 0], tol)) and any(findclose(d, dgmX[k][:, 0], tol)):
+                    coker.append((b, d))
+
+                # ker
+                if k > 0:
+                    # b_c = b_x_i (dim-1)
+                    # d_c = d_x_i (dim-1)
+                    m = findclose(b, dgmX[k - 1][:, 0], tol) & findclose(d, dgmX[k - 1][:, 1], tol)
+                    if sum(m):
+                        ker.append((b, d))
+
+                    # b_c = d_y_i (dim-1)
+                    # d_c = d_x_j (dim-1)
+                    if any(findclose(b, dgmY[k - 1][:, 1], tol)) and any(findclose(d, dgmX[k - 1][:, 1], tol)):
+                        ker.append((b, d))
+
+        for c in coker:
+            data.append({"dim": k, "set": "coker", "b": c[0], "d": c[1]})
+        for c in ker:
+            data.append({"dim": k - 1, "set": "ker", "b": c[0], "d": c[1]})
+
+    return data

tests/test_kercoker.py

@@ -0,0 +1,38 @@
+import numpy as np
+from totalpersistence.totalpersistence import kercoker_via_cone, general_position_distance_matrix
+
+
+def test_kercoker_via_cone():
+    # Test input data
+    X = np.array([[0, 0], [0, 1], [1, 0], [1, 1], [2, 0], [2, 1]])
+    Y = np.array([[0, 2], [1, 2], [2, 2]])
+    f = np.array([0, 0, 1, 1, 2, 2])
+
+    # Compute perturbed distance matrices
+    dX = general_position_distance_matrix(X, perturb=1e-13)
+    dY = general_position_distance_matrix(Y, perturb=1e-13)
+
+    # Run the function
+    result = kercoker_via_cone(dX, dY, f, maxdim=2, cone_eps=0)[0]
+    print(result)
+
+    # Expected output
+    expected = [
+        {"dim": 0, "set": "coker", "b": 0.0, "d": float("inf")},
+        {"dim": 0, "set": "ker", "b": 0.5, "d": 1.0},
+        {"dim": 0, "set": "ker", "b": 0.5, "d": 1.0},
+        {"dim": 0, "set": "ker", "b": 0.0, "d": 1.0},
+        {"dim": 0, "set": "ker", "b": 0.0, "d": 1.0},
+        {"dim": 0, "set": "ker", "b": 0.0, "d": 1.0},
+        {"dim": 1, "set": "ker", "b": 1.0, "d": 1.4142135381698608},
+        {"dim": 1, "set": "ker", "b": 1.0, "d": 1.4142135381698608},
+    ]
+    tol = 1e-10
+
+    # Compare results
+    assert len(result) == len(expected)
+    for r, e in zip(result, expected):
+        assert r["dim"] == e["dim"]
+        assert r["set"] == e["set"]
+        assert abs(r["b"] - e["b"]) < tol
+        assert abs(r["d"] - e["d"]) < tol if e["d"] != float("inf") else r["d"] == float("inf")