PyGSP

PyGSP

# PyGSP (0.5.1)
# matplotlib (3.1.2)
# networkx (2.4)
# numpy (1.17.4)

from pygsp import graphs, filters
import matplotlib.pyplot as plt
G = graphs.Logo()
G.estimate_lmax()
g = filters.Heat(G, tau=100)
import numpy as np
DELTAS = [20, 30, 1090]
s = np.zeros(G.N)
s[DELTAS] = 1
s = g.filter(s)
G.plot_signal(s, highlight=DELTAS, backend='matplotlib')
plt.show()
###########################################################
import numpy as np
import matplotlib.pyplot as plt
from pygsp import graphs, filters, plotting

plotting.BACKEND = 'matplotlib'
plt.rcParams['figure.figsize'] = (10, 5)

rs = np.random.RandomState(42) # Reproducible results.
W = rs.uniform(size=(30, 30)) # Full graph.
W[W < 0.93] = 0 # Sparse graph.
W = W + W.T # Symmetric graph.
np.fill_diagonal(W, 0) # No self-loops.
G = graphs.Graph(W)
print('{} nodes, {} edges'.format(G.N, G.Ne))

connected_flg = G.is_connected()
directed_flg = G.is_directed()

# We can retrieve our weight matrix, which is stored in a sparse format.
retrieve_flg = (G.W == W).all()
W_type = type(G.W)

# The graph Laplacian (combinatorial by default).
# Laplacian
G.L.shape#(30, 30)

# We can also compute and get the graph Fourier basis
# the Fourier basis
G.compute_fourier_basis()
G.U.shape# (30, 30)

# the graph differential operator, useful to e.g. compute the gradient or smoothness of a signal.
# the differential operator
G.compute_differential_operator()
G.D.shape#(60, 30)

# To be able to plot a graph, we need to embed its nodes in a 2D or 3D space.
# Let’s set some coordinates with pygsp.graphs.Graph.set_coordinates() and plot our graph.
G.set_coordinates('ring2D')
G.plot()
plt.show()
##############################################################
# As in classical signal processing, the Fourier transform
# plays a central role in graph signal processing. Getting
# the Fourier basis is however computationally intensive as
# it needs to fully diagonalize the Laplacian. While it can
# be used to filter signals on graphs, a better alternative
# is to use one of the fast approximations (see pygsp.filters.Filter.filter()).
# Let’s plot the second and third eigenvectors (the first is constant).
G = graphs.Logo()
G.compute_fourier_basis()
fig, axes = plt.subplots(1, 2, figsize=(10, 3))
for i, ax in enumerate(axes):
G.plot_signal(G.U[:, i+1], vertex_size=30, ax=ax)
_ = ax.set_title('Eigenvector {}'.format(i+2))
ax.set_axis_off()
fig.tight_layout()
plt.show()

G2 = graphs.Ring(N=50)
G2.compute_fourier_basis()
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
G2.plot_signal(G2.U[:, 4], ax=axes[0])
G2.set_coordinates('line1D')
G2.plot_signal(G2.U[:, 1:4], ax=axes[1])
fig.tight_layout()
plt.show()
###############################################################
# Filters
# To filter signals on graphs, we need to define filters.
# They are represented in the toolbox by the pygsp.filters.Filter class.
# Filters are usually defined in the spectral domain.
# let’s define and plot that low-pass filter:
tau = 1
def g(x):
return 1. / (1. + tau * x)
g = filters.Filter(G, g)
fig, ax = plt.subplots()
g.plot(plot_eigenvalues=True, ax=ax)
_ = ax.set_title('Filter frequency response')
plt.show()
###############################################################
# Let’s create a graph signal and add some random noise.
# Graph signal: each letter gets a different value + additive noise.
s = np.zeros(G.N)
s[G.info['idx_g']-1] = -1
s[G.info['idx_s']-1] = 0
s[G.info['idx_p']-1] = 1
G.plot()
plt.show()
s += rs.uniform(-0.5, 0.5, size=G.N)
# We can now try to denoise that signal by filtering it with the above defined low-pass filter.
s2 = g.filter(s)
fig, axes = plt.subplots(1, 2, figsize=(10, 3))
G.plot_signal(s, vertex_size=30, ax=axes[0])
_ = axes[0].set_title('Noisy signal')
axes[0].set_axis_off()
G.plot_signal(s2, vertex_size=30, ax=axes[1])
_ = axes[1].set_title('Cleaned signal')
axes[1].set_axis_off()
fig.tight_layout()
plt.show()
# While the noise is largely removed thanks to the filter, some energy is diffused between the letters.
###############################################################
###############################################################
# Next contents will show you how to easily construct a wavelet frame,
# a kind of filter bank, and apply it to a signal.
# This tutorial will walk you into computing the wavelet coefficients
# of a graph, visualizing filters in the vertex domain, and using the
# wavelets to estimate the curvature of a 3D shape.
###############################################################
# spectral graph wavelets
# 显示3d必须导入 from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from pygsp import graphs, filters, plotting, optimization, utils
from mpl_toolkits.mplot3d import Axes3D

G = graphs.Bunny()
# estimate the largest eigenvalue G.lmax
lamida_max_L = G.estimate_lmax()

# Simple filtering: heat diffusion
taus = [10, 25, 50]
g = filters.Heat(G, taus)
# create a signal as a Kronecker delta located on one vertex,
# e.g. the vertex 20. That signal is our heat source.
s = np.zeros(G.N)
DELTA = 20
s[DELTA] = 1
# We can now simulate heat diffusion by filtering our signal s with each of our heat kernels.
s = g.filter(s, method='chebyshev')
# finally plot the filtered signal showing heat diffusion at different scales.
fig = plt.figure(figsize=(10, 3))
num_fitter = g.Nf
for i in range(num_fitter):
# 只要Axes3D导入存在就行 from mpl_toolkits.mplot3d import Axes3D
ax = fig.add_subplot(1, num_fitter, i + 1, projection='3d')
G.plot_signal(s[:, i], colorbar=True, ax=ax)
title = r'Heat diffusion, $	au={}$'.format(taus[i])
_ = ax.set_title(title)
ax.set_axis_off()
fig.tight_layout()
plt.show()

# We can visualize the atoms as we did with the heat kernel,
# by filtering a Kronecker delta placed at one specific vertex.
s = g.localize(DELTA)
fig = plt.figure(figsize=(10, 2.5))
for i in range(num_fitter):
ax = fig.add_subplot(1, 3, i + 1, projection='3d')
G.plot_signal(s[:, i], ax=ax)
_ = ax.set_title('Wavelet {}'.format(i + 1))
ax.set_axis_off()
fig.tight_layout()
plt.show()
###################################################################
# Curvature estimation
# As a last and more applied example, let us try to estimate
# the curvature of the underlying 3D model by only using spectral
# filtering on the nearest-neighbor graph formed by its point cloud.

# let us try to estimate the curvature of the underlying 3D model
# by only using spectral filtering on the nearest-neighbor graph formed by its point cloud.
# Doing so gives us a 3-dimensional signal:
'''
s = G.coords
s = g.filter(s)
# The curvature is then estimated by taking the ℓ1 or ℓ2 norm across the 3D position.
s = np.linalg.norm(s, ord=2, axis=1)
fig = plt.figure(figsize=(10, 7))
for i in range(4):
ax = fig.add_subplot(2, 2, i + 1, projection='3d')
G.plot_signal(s[:, i], ax=ax)
title = 'Curvature estimation (scale {})'.format(i + 1)
_ = ax.set_title(title)
ax.set_axis_off()
fig.tight_layout()
plt.show()
'''
###################################################################
# The pygsp.filters.Filter.localize() method can be used to
# visualize a filter in the vertex domain instead of doing it manually.
# localize(i, **kwargs)# Localize the kernels at a node (to visualize them).
# i: Index of the node where to localize the kernel.
import matplotlib
N = 20
DELTA = N//2 * (N+1)
G = graphs.Grid2d(N)
G.estimate_lmax()
g = filters.Heat(G, 100)
s = g.localize(DELTA)
G.plot_signal(s, highlight=DELTA)
plt.show()
#############################################
# Visualizing wavelets atoms
# Let’s now replace the Heat filter by a filter bank of wavelets.
# We can create a filter bank using one of the predefined filters,
# such as pygsp.filters.MexicanHat to design a set of Mexican hat wavelets.
g = filters.MexicanHat(G, Nf=6) # Nf = 6 filters in the filter bank.
fig, ax = plt.subplots(figsize=(10, 5))
g.plot(ax=ax)
_ = ax.set_title('Filter bank of mexican hat wavelets')
plt.show()

# A better coverage could be obtained by adapting the filter bank with
# pygsp.filters.WarpedTranslates or by using another filter bank like pygsp.filters.Itersine.
# pygsp.filters.Itersine(G, Nf=6, overlap=2.0, **kwargs)
# Create an itersine half overlap filter bank of Nf filters. Going from 0 to lambda_max.

# Filter bank’s representation in Fourier and time (ring graph) domains.
import matplotlib.pyplot as plt
G = graphs.Ring(N=20)
G.estimate_lmax()
G.set_coordinates('line1D')
g = filters.HalfCosine(G)
s = g.localize(G.N // 2)
fig, axes = plt.subplots(1, 2)
g.plot(ax=axes[0])
G.plot_signal(s, ax=axes[1])
plt.show()

# class Meyer(Filter):# Use of this kernel for SGWT proposed by
# Nora Leonardi and Dimitri Van De Ville in :cite:`leonardi2011wavelet`.
import matplotlib.pyplot as plt
G = graphs.Ring(N=20)
G.estimate_lmax()
G.set_coordinates('line1D')
g = filters.Meyer(G)
s = g.localize(G.N // 2)
fig, axes = plt.subplots(1, 2)
g.plot(ax=axes[0])
G.plot_signal(s, ax=axes[1])
plt.show()

#########################
import numpy as np
from pygsp import graphs, plotting

# Create a random sensor graph
G = graphs.Sensor(N=256, distribute=True, seed=42)
G.compute_fourier_basis()
# Create label signal
label_signal = np.copysign(np.ones(G.N), G.U[:, 3])
G.plot_signal(label_signal)
plt.show()
# The up figure shows a plot of the original label signal, that we wish to recover, on the graph.

rs = np.random.RandomState(42)
# Create the mask
M = rs.rand(G.N)
M = (M > 0.6).astype(float) # Probability of having no label on a vertex.
# Applying the mask to the data
sigma = 0.1
subsampled_noisy_label_signal = M * (label_signal + sigma * rs.standard_normal(G.N))
G.plot_signal(subsampled_noisy_label_signal)
plt.show()
# The up figure shows the label signal on the graph after the application of the subsampling mask and the addition of noise. The label of more than half of the vertices has been set to 00.