%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);
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
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')
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]))
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))
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))
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)))
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))
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, '.');