.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/vietoris_rips/plot_persistent_homology_dragon.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_persistent_homology_dragon.py: .. _vietoris_rips_dragon: Persistent Homology: Stanford Dragon ======================================== This tutorial covers persistent homology of point clouds, and filtered Vietoris-Rips complexes more generally. .. GENERATED FROM PYTHON SOURCE LINES 12-17 .. code-block:: Python import oat_python as oat import plotly.graph_objects as go import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 18-20 Load the point cloud ------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 22-24 For this tutorial we will work with a point cloud of 1000 points, sampled from the Stanford Dragon. This cell will load the points and plot them. .. GENERATED FROM PYTHON SOURCE LINES 24-44 .. code-block:: Python # Load the point cloud points = oat.point_cloud.stanford_dragon() # Plot trace = go.Scatter3d( x=points[:,0], y=points[:,1], z=points[:,2], mode="markers", marker = dict(opacity=1.0, size=3, color=points[:,1], colorscale="Peach") ) fig = go.Figure(data=trace) fig.update_layout( title=dict(text="Stanford dragon, 1000 points"), template="plotly_dark", height=800, ) fig .. raw:: html


.. GENERATED FROM PYTHON SOURCE LINES 45-52 Compute persistent homology ----------------------------------- 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 boundary matrix `D` which has been factored into a U-match decomposition `TM = DS`. This decomposition contains everything we need for persistent homology. .. GENERATED FROM PYTHON SOURCE LINES 52-70 .. code-block:: Python # compute the minimum enclosing radius; all homology vanishes above this filtration parameter enclosing = oat.dissimilarity.enclosing_radius_for_points(points) # construct a distance matrix where values over enclosing + 0.0000000001 are removed dissimilarity_matrix = oat.dissimilarity.sparse_matrix_for_points( points = points, max_dissimilarity = enclosing + 0.0000000001, # adding 0.0000000001 avoids problems due to numerical error ) # build and factor the boundary matrix decomposition = oat.core.vietoris_rips.BoundaryMatrixDecomposition( dissimilarity_matrix = dissimilarity_matrix, max_homology_dimension = 1, support_fast_column_lookup = True, ) .. GENERATED FROM PYTHON SOURCE LINES 71-72 Export a summary of the persistent homology to a data frame. .. GENERATED FROM PYTHON SOURCE LINES 72-80 .. code-block:: Python persistent_homology_dataframe \ = decomposition.persistent_homology_dataframe( return_cycle_representatives = True, return_bounding_chains = True, ) persistent_homology_dataframe .. raw:: html
dimension interval_length birth_filtration birth_simplex death_filtration death_simplex cycle_representative num_cycle_simplices bounding_chain num_bounding_simplices
id
0 0 inf 0.000000 (0,) inf None simplex filtration coefficient 0 (0,) ... 1 None NaN
1 0 0.010316 0.000000 (1,) 0.010316 (649, 679) simplex filtration coefficient 0 (0,) ... 2 simplex filtration coefficient 0 (68... 97.0
2 0 0.009326 0.000000 (2,) 0.009326 (286, 444) simplex filtration coefficient 0 (0,) ... 2 simplex filtration coefficient 0 (23... 67.0
3 0 0.005452 0.000000 (3,) 0.005452 (3, 317) simplex filtration coefficient 0 (0,) ... 2 simplex filtration coefficient 0 (355,... 4.0
4 0 0.009618 0.000000 (4,) 0.009618 (204, 437) simplex filtration coefficient 0 (1,) ... 2 simplex filtration coefficient 0 (77... 74.0
... ... ... ... ... ... ... ... ... ... ...
1306 1 0.001114 0.020244 (67, 628) 0.021358 (67, 256, 670) simplex filtration coefficient 0 (628,... 5 simplex filtration coefficient 0 (... 3.0
1307 1 0.001493 0.021049 (670, 828) 0.022542 (670, 828, 835) simplex filtration coefficient 0 (42,... 9 simplex filtration coefficient 0 ... 7.0
1308 1 0.000646 0.021578 (670, 835) 0.022224 (670, 736, 835) simplex filtration coefficient 0 (11,... 5 simplex filtration coefficient 0 ... 3.0
1309 1 0.000099 0.021842 (264, 718) 0.021941 (362, 545, 718) simplex filtration coefficient 0 (362,... 4 simplex filtration coefficient 0 ... 2.0
1310 1 0.000750 0.022640 (97, 564) 0.023391 (97, 564, 797) simplex filtration coefficient 0 (97, 7... 4 simplex filtration coefficient 0 (... 2.0

1311 rows × 10 columns



.. GENERATED FROM PYTHON SOURCE LINES 81-83 Plot the persistence diagram ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 83-89 .. code-block:: Python fig_pd = oat.plot.persistence_diagram( persistent_homology_dataframe ) fig_pd .. raw:: html


.. GENERATED FROM PYTHON SOURCE LINES 90-92 Plot the barcode ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 94-99 .. code-block:: Python fig_barcode = oat.plot.barcode( persistent_homology_dataframe ) fig_barcode .. raw:: html


.. GENERATED FROM PYTHON SOURCE LINES 100-104 Inspect homology and cycle representatives ----------------------------------------------- The `homology` object is a data frame .. GENERATED FROM PYTHON SOURCE LINES 106-108 .. code-block:: Python persistent_homology_dataframe .. raw:: html
dimension interval_length birth_filtration birth_simplex death_filtration death_simplex cycle_representative num_cycle_simplices bounding_chain num_bounding_simplices
id
0 0 inf 0.000000 (0,) inf None simplex filtration coefficient 0 (0,) ... 1 None NaN
1 0 0.010316 0.000000 (1,) 0.010316 (649, 679) simplex filtration coefficient 0 (0,) ... 2 simplex filtration coefficient 0 (68... 97.0
2 0 0.009326 0.000000 (2,) 0.009326 (286, 444) simplex filtration coefficient 0 (0,) ... 2 simplex filtration coefficient 0 (23... 67.0
3 0 0.005452 0.000000 (3,) 0.005452 (3, 317) simplex filtration coefficient 0 (0,) ... 2 simplex filtration coefficient 0 (355,... 4.0
4 0 0.009618 0.000000 (4,) 0.009618 (204, 437) simplex filtration coefficient 0 (1,) ... 2 simplex filtration coefficient 0 (77... 74.0
... ... ... ... ... ... ... ... ... ... ...
1306 1 0.001114 0.020244 (67, 628) 0.021358 (67, 256, 670) simplex filtration coefficient 0 (628,... 5 simplex filtration coefficient 0 (... 3.0
1307 1 0.001493 0.021049 (670, 828) 0.022542 (670, 828, 835) simplex filtration coefficient 0 (42,... 9 simplex filtration coefficient 0 ... 7.0
1308 1 0.000646 0.021578 (670, 835) 0.022224 (670, 736, 835) simplex filtration coefficient 0 (11,... 5 simplex filtration coefficient 0 ... 3.0
1309 1 0.000099 0.021842 (264, 718) 0.021941 (362, 545, 718) simplex filtration coefficient 0 (362,... 4 simplex filtration coefficient 0 ... 2.0
1310 1 0.000750 0.022640 (97, 564) 0.023391 (97, 564, 797) simplex filtration coefficient 0 (97, 7... 4 simplex filtration coefficient 0 (... 2.0

1311 rows × 10 columns



.. GENERATED FROM PYTHON SOURCE LINES 109-113 Inspect a cycle representative and its bounding chain --------------------------------------------------------- This dataframe is sorted by the values in the `filtration` column, with ties broken by lexicographic order on simplices. .. GENERATED FROM PYTHON SOURCE LINES 115-117 .. code-block:: Python persistent_homology_dataframe.cycle_representative[875] .. raw:: html
simplex filtration coefficient
0 (368,) 0.0 -1
1 (875,) 0.0 1


.. GENERATED FROM PYTHON SOURCE LINES 118-120 .. code-block:: Python persistent_homology_dataframe.bounding_chain[875] .. raw:: html
simplex filtration coefficient
0 (875, 876) 0.005440 -1
1 (368, 876) 0.006727 1


.. GENERATED FROM PYTHON SOURCE LINES 121-134 Plot a representative ----------------------------------------------- .. note:: Check out the :ref:`Plotting tutorials ` for more resources on plotting, especially - :ref:`styling_3d_gallery` - :ref:`cycle_representative_strategies_gallery` .. GENERATED FROM PYTHON SOURCE LINES 137-138 We'll plot a cycle representative for the following row of the persistent homology dataframe: .. GENERATED FROM PYTHON SOURCE LINES 138-140 .. code-block:: Python feature_number = 1203 .. GENERATED FROM PYTHON SOURCE LINES 141-180 .. code-block:: Python edges = persistent_homology_dataframe["cycle_representative"][feature_number]["simplex"].tolist() # the cycle triangles = persistent_homology_dataframe["bounding_chain"][feature_number]["simplex"].tolist() # the chain that bounds the cycle points = points trace_edge = oat.plot.trace_3d_for_edges( edges=edges, points=points, line=dict(color="white", width=10), name="Cycle", ) trace_triangle = oat.plot.trace_3d_for_triangles( triangles=triangles, points=points, showlegend=True, opacity=0.4, color="white", name="Bounding Chain", ) trace_points = go.Scatter3d( x=points[:,0], y=points[:,1], z=points[:,2], showlegend=True, opacity=0.5, mode="markers", marker = dict(opacity=0.8, size=3, color=points[:,1], colorscale="Peach"), name="Point Cloud", ) fig = go.Figure(data= [ trace_edge, trace_triangle, trace_points ] ) fig.update_layout( title=dict(text="Cycle representative and bounding chain"), template="plotly_dark", height=800, ) fig.update_layout() fig .. raw:: html


.. GENERATED FROM PYTHON SOURCE LINES 181-189 Compare with an optimal cycle representative ----------------------------------------------- OAT offers tools to compute optimal cycle representatives. Let's compare with an optimized version of the same cycle. .. note:: Check out the documentation for :func:`oat_python.plot.contrast_initial_and_optimal_cycles_in_3d` for details, and options to customize the plot. .. GENERATED FROM PYTHON SOURCE LINES 191-202 .. code-block:: Python fig = oat.plot.contrast_initial_and_optimal_cycles_in_3d( boundary_matrix_decomposition = decomposition, birth_simplex = persistent_homology_dataframe["birth_simplex"][feature_number], points = points, ) fig.update_layout( margin=dict(l=20, r=10, t=150, b=10), width=None, # set the width to automatic ) fig .. raw:: html


.. GENERATED FROM PYTHON SOURCE LINES 203-207 Analyze an optimal representative ----------------------------------------------- Here's how to compute the optimal cycle, and inspect its data frame. First, compute the optimal cycle: .. GENERATED FROM PYTHON SOURCE LINES 207-214 .. code-block:: Python optimal_cycle_data = decomposition.optimize_cycle( birth_simplex = persistent_homology_dataframe["birth_simplex"][feature_number], problem_type = "preserve PH basis", ) optimal_cycle_data .. raw:: html
cost num_nonzero_coefficients chain
variable
initial cycle 0.607079 75.0 simplex filtration coefficient 0 (84...
optimal_cycle 0.184934 21.0 simplex filtration coefficient 0 (23...
x NaN 200.0 simplex filtration coefficient 0 (...
surface_between_cycles NaN 522.0 simplex filtration coefficient 0...
difference_in_essential_cycles NaN 0.0 Empty DataFrame Columns: [simplex, filtration,...
time_to_formulate_the_problem 0.013879 NaN None
time_to_solve_the_problem 1.714217 NaN None


.. GENERATED FROM PYTHON SOURCE LINES 215-240 The output is a dataframe that contains the solution, as well as several other pieces of information. - ``initial_cycle``: the initial cycle representative returned by the ``decomposition`` - ``optimal_cycle``: the optimal cycle representative. This cycle has form .. math:: o = z + e + Dc where - :math:`o` is the optimal cycle, - :math:`z` is the initial cycle, - :math:`c` is a chain and :math:`Dc` is the boundary of :math:`c`, - :math:`e` is a chain in the space spanned by essential cycles (that is, cycles which never become boundaries). - Typically :math:`e` is zero, so :math:`o = z + Dc`. In fact this is will *always* true if :math:`z` is non-essential, i.e. if :math:`z` represents persistent homology class with a finite death time. If this is true, then :math:`z` and `o` will eventually become homologous, since they differ by a boundary. However, :math:`c` may have a birth time strictly later than :math:`z`, so :math:`z` and `o` may not be homologous at the birth time of :math:`z`. - ``surface_between_cycles`` is a chain :math:`c`. You can think of this chain, informally, as a surface whose boundary is the difference between the initial and optimal cycles. In particular, if :math:`e=0` then :math:`0 = z + Dc`. - ``difference_in_essential_cycles``: is the chain :math:`e` in the decomposition above. Notice that the optimal cycle has just 21 nonzero coefficients, while the initial cycle has 75! .. GENERATED FROM PYTHON SOURCE LINES 242-243 Display the data frame for the optimal cycle .. GENERATED FROM PYTHON SOURCE LINES 243-244 .. code-block:: Python optimal_cycle_data["chain"]["optimal_cycle"] .. raw:: html
simplex filtration coefficient
0 (231, 325) 0.004040 -1
1 (595, 818) 0.005324 -1
2 (564, 645) 0.005361 -1
3 (93, 616) 0.006633 1
4 (459, 964) 0.006766 -1
5 (369, 964) 0.007283 1
6 (393, 438) 0.007515 -1
7 (936, 995) 0.008283 1
8 (678, 735) 0.008574 -1
9 (111, 645) 0.008814 1
10 (438, 459) 0.009093 -1
11 (111, 250) 0.009401 -1
12 (616, 818) 0.009411 1
13 (72, 393) 0.009625 -1
14 (369, 995) 0.010594 -1
15 (346, 735) 0.010838 1
16 (346, 564) 0.011005 -1
17 (678, 936) 0.011470 1
18 (250, 595) 0.011514 -1
19 (93, 231) 0.011602 -1
20 (72, 325) 0.011787 1


.. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.604 seconds) .. _sphx_glr_download_auto_examples_vietoris_rips_plot_persistent_homology_dragon.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_persistent_homology_dragon.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_persistent_homology_dragon.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_persistent_homology_dragon.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_