.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/vietoris_rips/plot_rdv.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_vietoris_rips_plot_rdv.py: .. _vietoris_rips_rdv: R=DV and other matrix decompositions ======================================== `R = DV factorization `_ is a standard tool to compute persistent homology. OAT obtains R=DV factorizations by way of Umatch decompositions. This tutorial shows how to obtain - Umatch decompositions - R=DV factorizations - R=DV factorizations of the anti-transpose of D (for persistent cohomology) - RU=D factorizations .. GENERATED FROM PYTHON SOURCE LINES 18-27 .. code-block:: Python import oat_python as oat from fractions import Fraction import matplotlib.pyplot as plt import numpy as np import pandas as pd import plotly.graph_objects as go from scipy.sparse import csr_matrix .. GENERATED FROM PYTHON SOURCE LINES 28-30 Setup ----------------- .. GENERATED FROM PYTHON SOURCE LINES 32-33 Define a function to check equality of sparse matrices .. GENERATED FROM PYTHON SOURCE LINES 33-44 .. code-block:: Python def sparse_matrices_equal(A, B): """Check if two scipy sparse matrices are equal.""" # Check shape first if A.shape != B.shape: return False # Subtract and check if all elements are zero diff = (A != B).nnz == 0 return diff .. GENERATED FROM PYTHON SOURCE LINES 45-46 Generate a point cloud .. GENERATED FROM PYTHON SOURCE LINES 46-49 .. code-block:: Python points = oat.point_cloud.sphere_or_slice_spiral(n_points=8, noise_scale=0.07, random_seed=0) .. GENERATED FROM PYTHON SOURCE LINES 50-51 Plot the point cloud .. GENERATED FROM PYTHON SOURCE LINES 51-56 .. code-block:: Python trace = go.Scatter3d(x=points[:,0],y=points[:,1],z=points[:,2], mode="markers", marker=dict(opacity=1, size=4, color=points[:,2], colorscale="Aggrnyl")) fig = go.Figure(data=trace) fig.update_layout(title=dict(text="Point cloud"), height=800,width=850) fig .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 57-59 Compute persistent homology ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 61-62 We compute persistent homology by factoring the boundary matrix. The following cell generates a sparse distance matrix and feeds it to the persistent homology solver. The result is a decomposition boundary matrix. We will extract information from this matrix in the following cells. .. GENERATED FROM PYTHON SOURCE LINES 64-65 Prepare a dissimilarity matrix .. GENERATED FROM PYTHON SOURCE LINES 65-80 .. code-block:: Python # Calculate the minimum enclosing radius. There is a theorem that states all # homology vanishes above this filtration parameter, so we can exclude # simplices whose diameter is larger than this value from the complex. enclosing_radius = oat.dissimilarity.enclosing_radius_for_points(points) # Format the dissimilarity matrix. # This produces a dissimilarity matrix where all values over # enclosing_radius + 0.0000000001 are removed. Adding 0.0000000001 avoids # problems due to numerical error. dissimilarity_matrix = oat.dissimilarity.sparse_matrix_for_points( points = points, max_dissimilarity = enclosing_radius + 0.00000000001, ) .. GENERATED FROM PYTHON SOURCE LINES 81-82 Create and decompose the boundary matrix .. GENERATED FROM PYTHON SOURCE LINES 82-87 .. code-block:: Python decomposition = oat.core.vietoris_rips.BoundaryMatrixDecomposition( dissimilarity_matrix = dissimilarity_matrix, ) .. GENERATED FROM PYTHON SOURCE LINES 88-89 Compute the boundary of a chain .. GENERATED FROM PYTHON SOURCE LINES 89-101 .. code-block:: Python chain = pd.DataFrame([ { "simplex":(0,1,2), "coefficient":Fraction(1, 1) }, { "simplex":(4,5,6), "coefficient":Fraction(1, 1) }, ]) boundary_matrix = decomposition.boundary_matrix_oracle() boundary = boundary_matrix.boundary_for_chain( chain ) oat.plot.display_dataframes_side_by_side( chain, boundary, titles= ["Chain", "Boundary for this chain"] ) .. raw:: html

Chain

	simplex	coefficient
0	(0, 1, 2)	1
1	(4, 5, 6)	1

Boundary for this chain

	simplex	filtration	coefficient
0	(0, 1)	0.573701	1
1	(0, 2)	0.974061	-1
2	(1, 2)	1.041700	1
3	(5, 6)	1.121170	1
4	(4, 5)	1.476147	1
5	(4, 6)	1.503094	-1

.. GENERATED FROM PYTHON SOURCE LINES 102-103 Sort the simplices in the boundary by lexicographic order .. GENERATED FROM PYTHON SOURCE LINES 103-107 .. code-block:: Python boundary.sort_values(inplace=True, by="simplex") boundary .. raw:: html

	simplex	filtration	coefficient
0	(0, 1)	0.573701	1
1	(0, 2)	0.974061	-1
2	(1, 2)	1.041700	1
4	(4, 5)	1.476147	1
5	(4, 6)	1.503094	-1
3	(5, 6)	1.121170	1

.. GENERATED FROM PYTHON SOURCE LINES 108-109 Compute a coboundary .. GENERATED FROM PYTHON SOURCE LINES 109-119 .. code-block:: Python cochain = pd.DataFrame([ { "simplex":(0,1), "coefficient":Fraction(1, 1) }, ]) coboundary = boundary_matrix.coboundary_for_cochain( cochain ) oat.plot.display_dataframes_side_by_side( cochain, coboundary, titles= ["Cochain", "Coboundary for this cochain"] ) .. raw:: html

Cochain

	simplex	coefficient
0	(0, 1)	1

Coboundary for this cochain

	simplex	filtration	coefficient
0	(0, 1, 2)	1.041700	1
1	(0, 1, 4)	1.493781	1
2	(0, 1, 3)	1.560769	1
3	(0, 1, 5)	1.643623	1
4	(0, 1, 6)	1.790893	1

.. GENERATED FROM PYTHON SOURCE LINES 120-122 A boundary matrix ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 124-133 .. code-block:: Python D = decomposition.boundary_matrix_as_csr() plt.figure(figsize=(6, 6)) plt.spy(D, markersize=2, color="orange") plt.title("D (boundary matrix)", y=-0.18) # Move title to the bottom plt.tight_layout() plt.show() .. image-sg:: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_001.png :alt: D (boundary matrix) :srcset: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 134-141 The simplex corresponding to each row and column of the boundary matrix can be obtained as follows. - This list includes all simplices of dimension <= max_homology_dimension, and all *negative* simplices of dimension max_homology_dimension + 1. - Simplices are sorted first by dimension, then by filtration value, then lexicographically. .. GENERATED FROM PYTHON SOURCE LINES 141-144 .. code-block:: Python decomposition.boundary_matrix_indices_df() .. raw:: html

	simplex	filtration
0	(0,)	0.000000
1	(1,)	0.000000
2	(2,)	0.000000
3	(3,)	0.000000
4	(4,)	0.000000
5	(5,)	0.000000
6	(6,)	0.000000
7	(7,)	0.000000
8	(6, 7)	0.551279
9	(0, 1)	0.573701
10	(0, 2)	0.974061
11	(1, 2)	1.041700
12	(5, 7)	1.049453
13	(2, 5)	1.088970
14	(5, 6)	1.121170
15	(1, 4)	1.170000
16	(3, 6)	1.181233
17	(1, 5)	1.294119
18	(0, 3)	1.299655
19	(4, 7)	1.299778
20	(2, 6)	1.329181
21	(2, 3)	1.438239
22	(4, 5)	1.476147
23	(0, 4)	1.493781
24	(4, 6)	1.503094
25	(3, 4)	1.529175
26	(3, 7)	1.545687
27	(1, 3)	1.560769
28	(0, 5)	1.643623
29	(2, 7)	1.651619
30	(1, 6)	1.719551
31	(0, 6)	1.790893

.. GENERATED FROM PYTHON SOURCE LINES 145-165 A differential Umatch decomposition ----------------------------------------- A differential Umatch decomposition is a tuple of matrices :math:`(J,M,D,J)`, where - :math:`D` is a boundary matrix - :math:`J` is a square upper triangular matrix with ones on the diagonal - :math:`M` is a generalized matching matrix See :term:`Differential Umatch Decomposition` for definitions and background. Umatch decompositions can be used to obtain - :math:`R=DV` decompositions by setting :math:`V=J` and :math:`R=DJ` - :math:`RU=D` decompositions by setting :math:`U = J^{-1}`. - :math:`R=WD` decompositions by setting :math:`W=J^{-1}` and :math:`R=J^{-1}D` - this is the calculation used to obtain persistent cohomology - it is commonly described as :math:`R=DV` decomposition of the antitranspose :math:`D^\perp` - :math:`YR=D` decompositions by setting :math:`Y = J`. .. GENERATED FROM PYTHON SOURCE LINES 167-192 .. code-block:: Python D = decomposition.boundary_matrix_as_csr() J = decomposition.differential_comb_as_csr() M = decomposition.generalized_matching_matrix_as_csr() D.data = D.data.astype(float) J.data = J.data.astype(float) M.data = M.data.astype(float) # plot sparsity patterns fig, axs = plt.subplots(1,4) fig.set_figwidth(12) axs[0].spy(J, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[1].spy(M, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[2].spy(D, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[3].spy(J, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[0].set_title("J", y=-0.2) axs[1].set_title("M", y=-0.2) axs[2].set_title("D", y=-0.2) axs[3].set_title("J", y=-0.2) plt.tight_layout() .. image-sg:: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_002.png :alt: J, M, D, J :srcset: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 193-194 Verify that that JM = DJ .. GENERATED FROM PYTHON SOURCE LINES 194-196 .. code-block:: Python sparse_matrices_equal( J@M, D@J ) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 197-198 Visualize the difference JM - DJ, just to be thorough: .. GENERATED FROM PYTHON SOURCE LINES 198-211 .. code-block:: Python fig, axs = plt.subplots(1,3) fig.set_figwidth(12) axs[0].spy(D @ J, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[1].spy(J @ M, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[2].spy( (J@M) - (D@J), precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[0].set_title("DJ", y=-0.2) axs[1].set_title("JM", y=-0.2) axs[2].set_title("JM - DJ", y=-0.2) plt.tight_layout() .. image-sg:: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_003.png :alt: DJ, JM, JM - DJ :srcset: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 212-215 # An R = DV decomposition As noted above, Umatch decompositions can be used to obtain :math:`R=DV` decompositions by setting :math:`V=J` and :math:`R=DJ`. .. GENERATED FROM PYTHON SOURCE LINES 217-218 Scipy has limited functionality for coefficients of type Fraction, so convert to float .. GENERATED FROM PYTHON SOURCE LINES 218-235 .. code-block:: Python J.data = J.data.astype(float) D.data = D.data.astype(float) # plot sparsity patterns fig, axs = plt.subplots(1,3) fig.set_figwidth(12) axs[0].spy(D @ J, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[1].spy(D, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[2].spy(J, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[0].set_title("R", y=-0.2) axs[1].set_title("D", y=-0.2) axs[2].set_title("V", y=-0.2) plt.tight_layout() .. image-sg:: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_004.png :alt: R, D, V :srcset: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 236-237 Verify that R: = DJ is in fact reducted .. GENERATED FROM PYTHON SOURCE LINES 237-257 .. code-block:: Python def bottom_nonzero_rows_are_unique(A: csr_matrix): """ Returns True if, for all distinct nonzero columns i and j of A, the bottom (last) nonzero element of each column occurs in a different row. """ # Find nonzero columns nonzero_cols = np.flatnonzero(A.getnnz(axis=0)) # For each nonzero column, find the row index of the bottom nonzero element bottom_rows = [] for col in nonzero_cols: col_data = A[:, col].tocoo() if col_data.row.size > 0: bottom_rows.append(col_data.row.max()) # Check if all bottom rows are unique return len(bottom_rows) == len(set(bottom_rows)) # Example usage: bottom_nonzero_rows_are_unique(D @J) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 258-262 An RU = D decomposition ------------------------------ As noted above, Umatch decompositions can be used to obtain :math:`RU = D` decompositions by setting :math:`U = J^{-1}` and :math:`R = DJ`. .. GENERATED FROM PYTHON SOURCE LINES 264-284 .. code-block:: Python Jinv = decomposition.differential_comb_inverse_as_csr() # Scipy has limited functionality for coefficients of type Fraction, so convert to float Jinv.data = Jinv.data.astype(float) D.data = D.data.astype(float) # plot sparsity patterns fig, axs = plt.subplots(1,3) fig.set_figwidth(12) axs[0].spy((D @ J) @ Jinv, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[1].spy(Jinv, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[2].spy(D @ J, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[0].set_title("R", y=-0.2) axs[1].set_title("U", y=-0.2) axs[2].set_title("D", y=-0.2) plt.tight_layout() .. image-sg:: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_005.png :alt: R, U, D :srcset: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 285-289 An :math:`R=DV` decomposition of :math:`D^\perp` (persistent cohomology) ------------------------------------------------------------------------------- This is equivalent to a matrix equation `R = WD` where `W` is invertible and upper triangular, and where `R` is reduced in the sense that the leading entry of every nonzero row occurs in a distinct column. .. GENERATED FROM PYTHON SOURCE LINES 291-304 .. code-block:: Python fig, axs = plt.subplots(1,3) fig.set_figwidth(12) axs[0].spy(Jinv @ D, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[1].spy(Jinv, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[2].spy(D, precision=0, marker=None, markersize=1, color="orange", aspect='equal', origin='upper') axs[0].set_title("R", y=-0.2) axs[1].set_title("W", y=-0.2) axs[2].set_title("D", y=-0.2) plt.tight_layout() .. image-sg:: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_006.png :alt: R, W, D :srcset: /auto_examples/vietoris_rips/images/sphx_glr_plot_rdv_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 305-306 Verify that :math:`R` is reduced in the correct sense. .. GENERATED FROM PYTHON SOURCE LINES 308-328 .. code-block:: Python def leading_nonzero_cols_are_unique(A: csr_matrix): """ Returns True if, for all distinct nonzero rows i and j of A, the leading (first) nonzero element of each row occurs in a different column. """ # Find nonzero rows nonzero_rows = np.flatnonzero(A.getnnz(axis=1)) # For each nonzero row, find the column index of the first nonzero element leading_cols = [] for row in nonzero_rows: row_data = A[row, :].tocoo() if row_data.col.size > 0: leading_cols.append(row_data.col.min()) # Check if all leading columns are unique return len(leading_cols) == len(set(leading_cols)) # Example usage: leading_nonzero_cols_are_unique( Jinv @ D) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.867 seconds) .. _sphx_glr_download_auto_examples_vietoris_rips_plot_rdv.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_rdv.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_rdv.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_rdv.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_