.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/dowker/plot_optimize.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_dowker_plot_optimize.py: .. _dowker_optimize_cycles_gallery: Optimal cycles ========================================= In this example we will - generate a point cloud - generate a collection of subsets (this can be viewed as a *cover*, and also as a *hypergraph*) - compute the homology of the associated Dowker complex - analyze cycle representatives - plot cycle representatives .. GENERATED FROM PYTHON SOURCE LINES 18-28 .. code-block:: Python import oat_python as oat import plotly.graph_objects as go import plotly.express as px import numpy as np import sklearn import sklearn.metrics from sklearn.neighbors import NearestNeighbors import itertools .. GENERATED FROM PYTHON SOURCE LINES 29-31 Generate a point cloud ------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 33-34 Set parameter values .. GENERATED FROM PYTHON SOURCE LINES 34-38 .. code-block:: Python n_points = 60 maxdis = None maxdim = 1 .. GENERATED FROM PYTHON SOURCE LINES 39-40 Generate a point cloud consisting of three noisy circles .. GENERATED FROM PYTHON SOURCE LINES 40-51 .. code-block:: Python points = [] for seed in range(3): circle = oat.point_cloud.annulus(n_points=n_points, inner_radius=1, outer_radius=2.5, random_seed=seed) circle[:,0] += 2 * np.cos( seed * 2 * np.pi / 4) circle[:,1] += 2 * np.sin( seed * 2 * np.pi / 4) points.append( circle ) points = np.concatenate(points) dismat = sklearn.metrics.pairwise_distances(points) dismat = ( dismat + dismat.T ) / 2 .. GENERATED FROM PYTHON SOURCE LINES 52-53 Plot the point cloud .. GENERATED FROM PYTHON SOURCE LINES 53-58 .. code-block:: Python trace = go.Scatter(x=points[:,0],y=points[:,1], mode="markers") fig = go.Figure([trace]) fig.update_layout( title=dict(text="Circle with noise"), height=650 ) fig .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 59-61 Choose a cover ------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 63-96 .. code-block:: Python radius_neighbor = 1.2; # hyperedges will be neighborhoods of vertices, with this radius radius_net = 1.0; # we'll make epsilon net with this value of epsilon # COMPUTE AN EPSILON NET net, _ = oat.dissimilarity.farthest_point_sampling( # the algorithm uses farthest point sampling metric_space = dict( point_cloud = points ), stopping_condition = dict( epsilon = radius_net ) ) # PLOT THE COVER data = [] data.append( go.Scatter(x=points[:,0],y=points[:,1], mode="markers", name="Cloud", showlegend=True) ) data.append( go.Scatter(x=points[net,0],y=points[net,1], mode="markers", marker=dict(symbol="triangle-up", size=10), name="Net", showlegend=True) ) for counter, v in enumerate(net): point = points[v] trace = oat.plot.ball_2d( point[0], point[1], radius=radius_neighbor, n_points=100 ) trace.update( opacity=0.2, name=f"Cover {counter}" ) data.append(trace) fig = go.Figure(data) fig.update_layout( title=dict(text="Hyperedges"), height=650, yaxis=dict( scaleanchor="x", # Anchor y-axis scaling to x-axis scaleratio=1.0 # Set the aspect ratio (1.0 for square aspect) ), ) fig .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 97-105 Compute homology ------------------------------------------------------- We'll compute the homology of - the dual hypergraph; that is, the hypergraph where vertices are balls, and for each vertex `v` we have a hyperedge that contains every ball to which `v` belongs - equivalently, the witness complex where every point is a witness and net points are landmarks The homology solver only accepts hypergraphs represented by a list of lists of integers, currently. If your hypergraph has a different format (e.g., if vertices are strings), then you can use some built-in tools to help translate back and forth between this format and list-of-list format; see the [rbs_reduced](rbs_reduced.ipynb) notebook for examples. .. GENERATED FROM PYTHON SOURCE LINES 107-108 Format the cover as a family of sorted lists: .. GENERATED FROM PYTHON SOURCE LINES 110-116 .. code-block:: Python # data structure holding the whole cloud net_wrapper = NearestNeighbors(n_neighbors=1, algorithm='ball_tree').fit( points[net] ) # witness complex where all points are witnesses, and net points are landmarks cover = [ sorted(list(x)) for x in net_wrapper.radius_neighbors( points, radius=radius_neighbor, return_distance=False, ) ] .. GENERATED FROM PYTHON SOURCE LINES 117-118 Factor the boundary matrix to compute homology: .. GENERATED FROM PYTHON SOURCE LINES 120-125 .. code-block:: Python decomposition = oat.core.dowker.BoundaryMatrixDecompositionDowker( dowker_simplices = cover, max_homology_dimension = 2 ) .. GENERATED FROM PYTHON SOURCE LINES 126-127 Inspect the Betti numbers .. GENERATED FROM PYTHON SOURCE LINES 129-131 .. code-block:: Python decomposition.betti_numbers() .. rst-class:: sphx-glr-script-out .. code-block:: none [1, 4, 0] .. GENERATED FROM PYTHON SOURCE LINES 132-133 Place a cycle basis for homology into a dataframe: .. GENERATED FROM PYTHON SOURCE LINES 135-138 .. code-block:: Python homology = decomposition.homology() homology .. raw:: html

	dimension	cycle_representative	cycle_representative_nonzero_coefficient_count	unique_simplex_id
0	0	simplex coefficient 0 (0,) 1	1	(0,)
1	1	simplex coefficient 0 (0, 15) -...	12	(7, 15)
2	1	simplex coefficient 0 (0, 10) -...	16	(7, 16)
3	1	simplex coefficient 0 (0, 15) -1 ...	7	(8, 15)
4	1	simplex coefficient 0 (0, 9) 1 ...	3	(9, 25)

.. GENERATED FROM PYTHON SOURCE LINES 139-140 Inspect a cycle representative .. GENERATED FROM PYTHON SOURCE LINES 142-144 .. code-block:: Python homology["cycle_representative"][1] .. raw:: html

	simplex	coefficient
0	(0, 15)	-1
1	(0, 17)	1
2	(2, 13)	1
3	(2, 14)	-1
4	(3, 8)	1
5	(3, 11)	-1
6	(4, 11)	1
7	(4, 17)	-1
8	(6, 8)	-1
9	(6, 14)	1
10	(7, 13)	-1
11	(7, 15)	1

.. GENERATED FROM PYTHON SOURCE LINES 145-147 Plot cycles ------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 149-176 .. code-block:: Python # DATA FOR THE POINT CLOUD data = [] data.append( go.Scatter(x=points[:,0],y=points[:,1], mode="markers", name="Cloud", showlegend=True) ) data.append( go.Scatter(x=points[net,0],y=points[net,1], mode="markers+text", marker=dict(symbol="triangle-up", size=10), name="Net", showlegend=True, text=[str(x) for x in range(len(net))], textposition='bottom center', ) ) # DATA FOR THE CYCLE colors = px.colors.qualitative.Plotly # specifies a (discrete) sequence of colors, represented by a list of strings for rownum, row in homology.iterrows(): if row["dimension"] != 1: continue trace_cycle = oat.plot.trace_2d_for_edges( edges = row["cycle_representative"].simplex, points = points[net], name = f"Cycle {rownum}", ) data.append(trace_cycle) # add the trace to our list # PLOT fig = go.Figure( data ) fig.update_layout(title="Coefficients still appear in the hover data", height=650 ) fig .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 177-179 Optimize ------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 181-185 Within homology class ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This method will find a minimal cycle representative within the same homology class .. GENERATED FROM PYTHON SOURCE LINES 187-193 .. code-block:: Python optimal = decomposition.optimize_cycle( unique_simplex_id=[7,15], problem_type="preserve homology class" ) optimal .. raw:: html

	cost	nnz	chain
type of chain
initial_cycle	12.0	12	simplex coefficient 0 (7, 15) ...
optimal_cycle	10.0	10	simplex coefficient 0 (15, 26) -...
surface_between_cycles	NaN	4	simplex coefficient 0 (7, 15, 26) ...
difference in essential chains	NaN	6	simplex coefficient 0 (15, 26) -...
Ax + z - y	NaN	6	simplex coefficient 0 (15, 26) -...

.. GENERATED FROM PYTHON SOURCE LINES 194-195 Inspect the initial cycle .. GENERATED FROM PYTHON SOURCE LINES 195-198 .. code-block:: Python optimal["chain"]["initial_cycle"] .. raw:: html

	simplex	coefficient
0	(7, 15)	1
1	(7, 13)	-1
2	(6, 14)	1
3	(6, 8)	-1
4	(4, 17)	-1
5	(4, 11)	1
6	(3, 11)	-1
7	(3, 8)	1
8	(2, 14)	-1
9	(2, 13)	1
10	(0, 17)	1
11	(0, 15)	-1

.. GENERATED FROM PYTHON SOURCE LINES 199-200 Inspect the optimal cycle .. GENERATED FROM PYTHON SOURCE LINES 200-203 .. code-block:: Python optimal["chain"]["optimal_cycle"] .. raw:: html

	simplex	coefficient
0	(15, 26)	-1
1	(14, 26)	1
2	(6, 14)	1
3	(6, 8)	-1
4	(4, 17)	-1
5	(4, 11)	1
6	(3, 11)	-1
7	(3, 8)	1
8	(0, 17)	1
9	(0, 15)	-1

.. GENERATED FROM PYTHON SOURCE LINES 204-205 Compare cycles visually .. GENERATED FROM PYTHON SOURCE LINES 207-238 .. code-block:: Python # DATA FOR THE POINT CLOUD data = [] data.append( go.Scatter(x=points[:,0],y=points[:,1], mode="markers", name="Cloud", showlegend=True) ) data.append( go.Scatter(x=points[net,0],y=points[net,1], mode="markers+text", marker=dict(symbol="triangle-up", size=10), name="Net", showlegend=True, text=[str(x) for x in range(len(net))], textposition='bottom center', ) ) # DATA FOR THE CYCLE colors = px.colors.qualitative.Plotly # specifies a (discrete) sequence of colors, represented by a list of strings for cycle_num, cycle_name in enumerate(["optimal_cycle", "initial_cycle"]): cycle_color = colors[ cycle_num % len(colors) ] for entrynum, entry in optimal["chain"][cycle_name].iterrows(): edge = entry["simplex"] coefficient = entry["coefficient"] trace = oat.plot.trace_2d_for_edge( edge=edge, points=points[net] ) trace.update( name=cycle_name, text=f"simplex {edge}
linear coefficent {coefficient}", opacity=0.7, line=dict(color=cycle_color,)) # customize appearance trace.update( legendgroup=cycle_num) # group edges that belong to the same cycle trace.update( showlegend = entrynum==0 ) data.append(trace) # append to the data group # PLOT fig = go.Figure( data ) fig.update_layout(title="Initial and optimal cycles", height=650 ) fig .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 239-243 Across homology classes ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The following method allows us to optimize the cycle representative by adding not only boundaries, but cycle represntatives for other homology classes. .. GENERATED FROM PYTHON SOURCE LINES 245-253 .. code-block:: Python optimal = decomposition.optimize_cycle( unique_simplex_id = (7,15), problem_type = "preserve homology basis (once)" ) optimal .. raw:: html

	cost	nnz	chain
type of chain
initial_cycle	12.0	12	simplex coefficient 0 (7, 15) ...
optimal_cycle	4.0	4	simplex coefficient 0 (15, 26) -...
surface_between_cycles	NaN	7	simplex coefficient 0 (8, 15, 22) ...
difference in essential chains	NaN	16	simplex coefficient 0 (15, 26) ...
Ax + z - y	NaN	11	simplex coefficient 0 (15, 26) ...

.. GENERATED FROM PYTHON SOURCE LINES 254-255 Compare cycles visually .. GENERATED FROM PYTHON SOURCE LINES 257-289 .. code-block:: Python # DATA FOR THE POINT CLOUD data = [] data.append( go.Scatter(x=points[:,0],y=points[:,1], mode="markers", name="Cloud", showlegend=True) ) data.append( go.Scatter(x=points[net,0],y=points[net,1], mode="markers+text", marker=dict(symbol="triangle-up", size=10), name="Net", showlegend=True, text=[str(x) for x in range(len(net))], textposition='bottom center', ) ) # DATA FOR THE CYCLE colors = px.colors.qualitative.Plotly # specifies a (discrete) sequence of colors, represented by a list of strings for cycle_num, cycle_name in enumerate(["optimal_cycle", "initial_cycle"]): cycle_color = colors[ cycle_num % len(colors) ] for entrynum, entry in optimal["chain"][cycle_name].iterrows(): edge = entry["simplex"] coefficient = entry["coefficient"] trace = oat.plot.trace_2d_for_edge( edge=edge, points=points[net] ) trace.update( name=cycle_name, text=f"simplex {edge}
linear coefficent {coefficient}", opacity=0.7, line=dict(color=cycle_color,)) # customize appearance trace.update( legendgroup=cycle_num) # group edges that belong to the same cycle trace.update( showlegend = entrynum==0 ) data.append(trace) # append to the data group # PLOT fig = go.Figure( data ) fig.update_layout(title="Initial and optimal cycles", height=650 ) fig .. raw:: html

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