From c4b0aeb20d920ba366a656a9aee7afe78871c9c7 Mon Sep 17 00:00:00 2001 From: Rémi Flamary Date: Tue, 25 Jun 2019 14:26:16 +0200 Subject: add fgw examples in doc --- docs/source/auto_examples/plot_barycenter_fgw.rst | 268 ++++++++++++++++++++++ 1 file changed, 268 insertions(+) create mode 100644 docs/source/auto_examples/plot_barycenter_fgw.rst (limited to 'docs/source/auto_examples/plot_barycenter_fgw.rst') diff --git a/docs/source/auto_examples/plot_barycenter_fgw.rst b/docs/source/auto_examples/plot_barycenter_fgw.rst new file mode 100644 index 0000000..2c44a65 --- /dev/null +++ b/docs/source/auto_examples/plot_barycenter_fgw.rst @@ -0,0 +1,268 @@ + + +.. _sphx_glr_auto_examples_plot_barycenter_fgw.py: + + +================================= +Plot graphs' barycenter using FGW +================================= + +This example illustrates the computation barycenter of labeled graphs using FGW + +Requires networkx >=2 + +.. [18] Vayer Titouan, Chapel Laetitia, Flamary R{'e}mi, Tavenard Romain + and Courty Nicolas + "Optimal Transport for structured data with application on graphs" + International Conference on Machine Learning (ICML). 2019. + + + + +.. code-block:: python + + + # Author: Titouan Vayer + # + # License: MIT License + + #%% load libraries + import numpy as np + import matplotlib.pyplot as plt + import networkx as nx + import math + from scipy.sparse.csgraph import shortest_path + import matplotlib.colors as mcol + from matplotlib import cm + from ot.gromov import fgw_barycenters + #%% Graph functions + + + def find_thresh(C, inf=0.5, sup=3, step=10): + """ Trick to find the adequate thresholds from where value of the C matrix are considered close enough to say that nodes are connected + Tthe threshold is found by a linesearch between values "inf" and "sup" with "step" thresholds tested. + The optimal threshold is the one which minimizes the reconstruction error between the shortest_path matrix coming from the thresholded adjency matrix + and the original matrix. + Parameters + ---------- + C : ndarray, shape (n_nodes,n_nodes) + The structure matrix to threshold + inf : float + The beginning of the linesearch + sup : float + The end of the linesearch + step : integer + Number of thresholds tested + """ + dist = [] + search = np.linspace(inf, sup, step) + for thresh in search: + Cprime = sp_to_adjency(C, 0, thresh) + SC = shortest_path(Cprime, method='D') + SC[SC == float('inf')] = 100 + dist.append(np.linalg.norm(SC - C)) + return search[np.argmin(dist)], dist + + + def sp_to_adjency(C, threshinf=0.2, threshsup=1.8): + """ Thresholds the structure matrix in order to compute an adjency matrix. + All values between threshinf and threshsup are considered representing connected nodes and set to 1. Else are set to 0 + Parameters + ---------- + C : ndarray, shape (n_nodes,n_nodes) + The structure matrix to threshold + threshinf : float + The minimum value of distance from which the new value is set to 1 + threshsup : float + The maximum value of distance from which the new value is set to 1 + Returns + ------- + C : ndarray, shape (n_nodes,n_nodes) + The threshold matrix. Each element is in {0,1} + """ + H = np.zeros_like(C) + np.fill_diagonal(H, np.diagonal(C)) + C = C - H + C = np.minimum(np.maximum(C, threshinf), threshsup) + C[C == threshsup] = 0 + C[C != 0] = 1 + + return C + + + def build_noisy_circular_graph(N=20, mu=0, sigma=0.3, with_noise=False, structure_noise=False, p=None): + """ Create a noisy circular graph + """ + g = nx.Graph() + g.add_nodes_from(list(range(N))) + for i in range(N): + noise = float(np.random.normal(mu, sigma, 1)) + if with_noise: + g.add_node(i, attr_name=math.sin((2 * i * math.pi / N)) + noise) + else: + g.add_node(i, attr_name=math.sin(2 * i * math.pi / N)) + g.add_edge(i, i + 1) + if structure_noise: + randomint = np.random.randint(0, p) + if randomint == 0: + if i <= N - 3: + g.add_edge(i, i + 2) + if i == N - 2: + g.add_edge(i, 0) + if i == N - 1: + g.add_edge(i, 1) + g.add_edge(N, 0) + noise = float(np.random.normal(mu, sigma, 1)) + if with_noise: + g.add_node(N, attr_name=math.sin((2 * N * math.pi / N)) + noise) + else: + g.add_node(N, attr_name=math.sin(2 * N * math.pi / N)) + return g + + + def graph_colors(nx_graph, vmin=0, vmax=7): + cnorm = mcol.Normalize(vmin=vmin, vmax=vmax) + cpick = cm.ScalarMappable(norm=cnorm, cmap='viridis') + cpick.set_array([]) + val_map = {} + for k, v in nx.get_node_attributes(nx_graph, 'attr_name').items(): + val_map[k] = cpick.to_rgba(v) + colors = [] + for node in nx_graph.nodes(): + colors.append(val_map[node]) + return colors + + + + + + + +Generate data +------------- + + + +.. code-block:: python + + + #%% circular dataset + # We build a dataset of noisy circular graphs. + # Noise is added on the structures by random connections and on the features by gaussian noise. + + + np.random.seed(30) + X0 = [] + for k in range(9): + X0.append(build_noisy_circular_graph(np.random.randint(15, 25), with_noise=True, structure_noise=True, p=3)) + + + + + + + +Plot data +--------- + + + +.. code-block:: python + + + #%% Plot graphs + + plt.figure(figsize=(8, 10)) + for i in range(len(X0)): + plt.subplot(3, 3, i + 1) + g = X0[i] + pos = nx.kamada_kawai_layout(g) + nx.draw(g, pos=pos, node_color=graph_colors(g, vmin=-1, vmax=1), with_labels=False, node_size=100) + plt.suptitle('Dataset of noisy graphs. Color indicates the label', fontsize=20) + plt.show() + + + + +.. image:: /auto_examples/images/sphx_glr_plot_barycenter_fgw_001.png + :align: center + + + + +Barycenter computation +---------------------- + + + +.. code-block:: python + + + #%% We compute the barycenter using FGW. Structure matrices are computed using the shortest_path distance in the graph + # Features distances are the euclidean distances + Cs = [shortest_path(nx.adjacency_matrix(x)) for x in X0] + ps = [np.ones(len(x.nodes())) / len(x.nodes()) for x in X0] + Ys = [np.array([v for (k, v) in nx.get_node_attributes(x, 'attr_name').items()]).reshape(-1, 1) for x in X0] + lambdas = np.array([np.ones(len(Ys)) / len(Ys)]).ravel() + sizebary = 15 # we choose a barycenter with 15 nodes + + A, C, log = fgw_barycenters(sizebary, Ys, Cs, ps, lambdas, alpha=0.95, log=True) + + + + + + + +Plot Barycenter +------------------------- + + + +.. code-block:: python + + + #%% Create the barycenter + bary = nx.from_numpy_matrix(sp_to_adjency(C, threshinf=0, threshsup=find_thresh(C, sup=100, step=100)[0])) + for i, v in enumerate(A.ravel()): + bary.add_node(i, attr_name=v) + + #%% + pos = nx.kamada_kawai_layout(bary) + nx.draw(bary, pos=pos, node_color=graph_colors(bary, vmin=-1, vmax=1), with_labels=False) + plt.suptitle('Barycenter', fontsize=20) + plt.show() + + + +.. image:: /auto_examples/images/sphx_glr_plot_barycenter_fgw_002.png + :align: center + + + + +**Total running time of the script:** ( 0 minutes 2.065 seconds) + + + +.. only :: html + + .. container:: sphx-glr-footer + + + .. container:: sphx-glr-download + + :download:`Download Python source code: plot_barycenter_fgw.py ` + + + + .. container:: sphx-glr-download + + :download:`Download Jupyter notebook: plot_barycenter_fgw.ipynb ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ -- cgit v1.2.3