%matplotlib inline
import numpy as np
from scipy import sparse
import scipy.sparse.linalg
import matplotlib as mpl
from matplotlib import pyplot as plt
import gudhi
import pygsp as pg
import sys
sys.path.append('..')
from data.s2_6_complex_to_laplacians import build_boundaries, build_laplacians
The simplicial complex is built from a Delaunay triangulation. See the CGAL documentation.
Alternatives:
scipy.spatial.Delaunay
.
Only for convex shapes?n_points = 100 # Approximate.
dim = 2
rs = np.random.RandomState(None)
points = np.concatenate([
[0.2, 1] * rs.uniform(size=(n_points//4, dim)),
[0.2, 0.8] + [0.6, 0.2] * rs.uniform(size=(n_points//6, dim)),
[0.2, 0] + [0.6, 0.2] * rs.uniform(size=(n_points//6, dim)),
[0.8, 0] + [0.2, 1] * rs.uniform(size=(n_points//4, dim)),
[0.4, 0.4] + [0.2, 0.2] * rs.uniform(size=(n_points//6, dim))
])
points = np.random.uniform(size=(n_points, dim))
#points = pg.graphs.Grid2d(10).coords
n_points = points.shape[0]
print(f'{n_points} points')
plt.scatter(*points.T);
100 points
from scipy import spatial
tri = spatial.Delaunay(points) print(f'{tri.simplices.shape[0]} triangles')
plt.triplot(points[:,0], points[:,1], tri.simplices);
ac = gudhi.AlphaComplex(points)
st = ac.create_simplex_tree()
before = st.num_simplices()
_ = st.prune_above_filtration(1e-2)
print(f'filtration: {before} => {st.num_simplices()} simplices')
assert st.num_vertices() == n_points
assert st.dimension() == dim
filtration: 573 => 498 simplices
simplices
is a list of dictionaries, with one dictionary per degree $k$.
Each dictionary maps a simplex (represented as a set of vertices) to an integer that will be its index in the boundary and Laplacian operators.
def extract_simplices(simplex_tree):
simplices = [dict() for _ in range(simplex_tree.dimension()+1)]
for simplex, _ in simplex_tree.get_skeleton(simplex_tree.dimension()):
k = len(simplex)
simplices[k-1][frozenset(simplex)] = len(simplices[k-1])
return simplices
simplices = extract_simplices(st)
for k, s in enumerate(simplices):
print(f'n_{k} = {len(s):,} {k}-simplices')
n_0 = 100 0-simplices n_1 = 251 1-simplices n_2 = 147 2-simplices
def get_positions(simplices, dim):
polygons = list()
for i, simplex in enumerate(simplices[dim].keys()):
assert simplices[dim][simplex] == i # Dictionary is ordered.
polygon = list()
for vertex in simplex:
polygon.append(points[vertex])
polygons.append(polygon)
return polygons
lines = get_positions(simplices, 1)
triangles = get_positions(simplices, 2)
def value2color(values):
values -= values.min()
values /= values.max()
return mpl.cm.viridis(values)
def plot_nodes(colors, ax=None, **kwargs):
if ax is None:
fig, ax = plt.subplots()
ax.scatter(points[:, 0], points[:, 1], c=colors, **kwargs)
return ax.figure, ax
def plot_edges(colors, ax=None, **kwargs):
if ax is None:
fig, ax = plt.subfigs()
colors = value2color(colors)
collection = mpl.collections.LineCollection(lines, colors=colors, **kwargs)
ax.add_collection(collection)
ax.autoscale()
def plot_triangles(colors, ax=None, **kwargs):
if ax is None:
fig, ax = plt.subfigs()
colors = value2color(colors)
for triangle, color in zip(triangles, colors):
triangle = plt.Polygon(triangle, color=color, **kwargs)
ax.add_patch(triangle)
ax.autoscale()
def plot_triangles_plain( ax=None, **kwargs):
if ax is None:
fig, ax = plt.subfigs()
colors = len(simplices[1])*['linen']
for triangle, color in zip(triangles, colors):
triangle = plt.Polygon(triangle, color=color, **kwargs)
ax.add_patch(triangle)
ax.autoscale()
boundaries = build_boundaries(simplices)
s1 = np.zeros(len(simplices[1]))
s1[100:150] = 1
#s1 = np.random.uniform(size=len(simplices[1]))
s0 = boundaries[0] @ s1
s2 = boundaries[1].T @ s1
fig, ax = plt.subplots()
plot_nodes(s0, ax, zorder=3)
plot_edges(s1, ax, zorder=2)
plot_triangles(s2, ax, zorder=1)
#plt.colorbar(ax.collections[0], ax=ax)
laplacians = build_laplacians(boundaries)
for k, laplacian in enumerate(laplacians):
print('{}-simplices: {:,} simplices, {:.2%} sparse'.format(k, laplacian.shape[0], laplacian.nnz/np.prod(laplacian.shape)))
assert laplacian.shape == (len(simplices[k]), len(simplices[k]))
0-simplices: 100 simplices, 6.02% sparse 1-simplices: 251 simplices, 2.59% sparse 2-simplices: 147 simplices, 2.49% sparse
def plot_eigen(laplacian, plot, n_eigenvectors=4):
values, vectors = sparse.linalg.eigsh(laplacian, n_eigenvectors, which='SM')
fix, axes = plt.subplots(1, n_eigenvectors+3, figsize=(20, 3))
axes[0].plot(values, '.')
axes[1].scatter(vectors[:, 0], vectors[:, 1])
axes[2].scatter(vectors[:, 2], vectors[:, 3])
for i in range(n_eigenvectors):
plot(vectors[:, i], axes[i+3])
axes[i+3].axis('off')
plot_eigen(laplacians[0], plot_nodes)
The edges indeed seem to be embeded in linear subspaces.
TODO:
plot_eigen(laplacians[1], plot_edges)
plot_eigen(laplacians[2], plot_triangles)
We study the following dynamical system.
$$\frac{\partial s(t)}{\partial t}=-L_k s(t)$$values, vectors = [], []
for j , laplacian in enumerate(laplacians):
w, v = np.linalg.eigh(laplacians[j].toarray())
values.append(w)
vectors.append(v)
class Heat:
def __init__(self, t, lmax):
self.t = t
self.lmax = lmax
def __call__(self, x):
return np.exp(-self.t * x / self.lmax)
def plot(self):
x = np.linspace(0, self.lmax)
plt.plot(x, self(x))
dim=0
heat = Heat(t=10, lmax=values[dim][-1])
heat.plot()
s = np.zeros(len(simplices[0]))
s[50:70] = 1
#s = np.random.uniform(size=len(simplices[0]))
def filter(d, signal):
return vectors[d] @ (heat(values[d]) * (vectors[d].T @ signal))
s = filter(0, s)
print(s.shape)
fig, ax = plot_nodes(s)
plt.colorbar(ax.collections[0], ax=ax)
#plot_edges(np.zeros(len(simplices[1])), ax=ax, zorder=0)
print(np.linalg.norm(s, 1))
(100,) 20.000000620898522
dim=1
iterations=1000
heat = Heat(t=iterations, lmax=values[dim][-1])
heat.plot()
s = np.zeros(len(simplices[dim]))
#s[100:150] =1
s = np.random.uniform(size=len(simplices[dim]))
init_cond=np.copy(s)
def filter(d, signal):
return vectors[d] @ (heat(values[d]) * (vectors[d].T @ signal))
s = filter(dim, s)
filtred_s=np.copy(s)
print(s.shape)
fig, ax = plt.subplots()
plot_edges(s, ax,norm=False, zorder=2)
plot_triangles_plain(ax, zorder=1)
#plt.colorbar(ax.collections[0], ax=ax)
print(np.linalg.norm(s, 1))
(251,) 139.46565826729983
As studied in Control Using Higher Order Laplacians in Network Topologies, Chapter IV.
if np.isclose(values[1][0],10**-10):
epsilon=0.0001
harm_dist=np.linalg.norm(vectors[1][0]-filtred_s, 2)
c=~np.isclose(values[1],10**-8)
t_epsilon=-1/values[1][c][0]*np.log(epsilon/(np.linalg.norm(init_cond, 2)*len(init_cond)))
print('The solution is epsilon-harmonic for all t > {}'.format(t_epsilon))
print('After {} iterations the distance from one harmonic representative is {}'.format(iterations,harm_dist))
else:
print('The solution converges to zero. After {} iterations the norm is {}'.format(iterations,np.linalg.norm(filtred_s, 2)))
The solution is epsilon-harmonic for all t > 183.92344235057902 After 1000 iterations the distance from one harmonic representative is 1.368369778518699
graph = pg.graphs.Grid2d(10)
g = pg.filters.Heat(graph, 100)
s = np.zeros(graph.N)
s[50:70] = 1
s = g.filter(s, method='exact')
graph.plot_signal(s)
print(np.linalg.norm(s, 1))
2020-11-04 19:48:30,224:[WARNING](pygsp.graphs.graph._check_fourier_properties): The eigenvalues vector G.e is not available, we need to compute the Fourier basis. Explicitly call G.compute_fourier_basis() once beforehand to suppress the warning.
19.999999999999993
Simplicial complex and ground truths from Control Using Higher Order Laplacians in Network Topologies, figure 3.
simplices = [
{
frozenset([0]): 0,
frozenset([1]): 1,
frozenset([2]): 2,
frozenset([3]): 3,
frozenset([4]): 4,
},
{
frozenset([0, 1]): 0,
frozenset([0, 2]): 1,
frozenset([1, 2]): 2,
frozenset([1, 3]): 3,
frozenset([2, 3]): 4,
frozenset([2, 4]): 5,
frozenset([3, 4]): 6,
},
{
frozenset([1, 2, 3]): 0,
},
]
boundaries = build_boundaries(simplices)
laplacians = build_laplacians(boundaries)
_, vectors = sparse.linalg.eigsh(laplacians[1], 2, which='SM')
print(vectors)
values, vectors = np.linalg.eigh(laplacians[1].toarray())
plt.plot(values, '.');
[[ 0.4999086 0.36199972] [-0.49990845 -0.36199963] [ 0.432163 0.06092054] [ 0.06774557 0.30107912] [-0.36441737 0.24015856] [ 0.2966719 -0.54123765] [-0.29667178 0.5412376 ]]