During extension algorithm one has to compute lover_link for every vertex in the complex. So let us implement search for the lower link first. It requires quite a lot of code: first we find a star, then link and finally lower link for the given simplex.
In [90]:
from itertools import combinations, chain
def simplex_closure(a):
"""Returns the generator that iterating over all subsimplices (of all dimensions) in the closure
of the simplex a. The simplex a is also included.
"""
return chain.from_iterable([combinations(a, l) for l in range(1, len(a) + 1)])
def closure(K):
"""Add all missing subsimplices to K in order to make it a simplicial complex."""
return list({s for a in K for s in simplex_closure(a)})
def contained(a, b):
"""Returns True is a is a subsimplex of b, False otherwise."""
return all((v in b for v in a))
def star(s, cx):
"""Return the set of all simplices in the cx that contais simplex s.
"""
return {p for p in cx if contained(s, p)}
def intersection(s1, s2):
"""Return the intersection of s1 and s2."""
return list(set(s1).intersection(s2))
def link(s, cx):
"""Returns link of the simplex s in the complex cx.
"""
# Link consists of all simplices from the closed star that have
# empty intersection with s.
return [c for c in closure(star(s, cx)) if not intersection(s, c)]
def simplex_value(s, f, aggregate):
"""Return the value of f on vertices of s
aggregated by the aggregate function.
"""
return aggregate([f[v] for v in s])
def lower_link(s, cx, f):
"""Return the lower link of the simplex s in the complex cx.
The dictionary f is the mapping from vertices (integers)
to the values on vertices.
"""
sval = simplex_value(s, f, min)
return [s for s in link(s, cx)
if simplex_value(s, f, max) < sval]
Let us test the above function on the simple example: full triangle with values 0, 1 and 2 on the vertices labeled with 1, 2 and 3.
In [87]:
K = closure([(1, 2, 3)])
f = {1: 0, 2: 1, 3: 2}
for v in (1, 2, 3):
print"{0}: {1}".format((v,), lower_link((v,), K, f))
Now let us implement an extension algorithm. We are leaving out the cancelling step for clarity.
In [91]:
def join(a, b):
"""Return the join of 2 simplices a and b."""
return tuple(sorted(set(a).union(b)))
def extend(K, f):
"""Extend the field to the complex K.
Function on vertices is given in f.
Returns the pair V, C, where V is the dictionary containing discrete gradient vector field
and C is the list of all critical cells.
"""
V = dict()
C = []
for v in (s for s in K if len(s)==1):
# Add your own code
pass
return V, C
Let us test the algorithm on the example from the previous step (full triangle).
In [92]:
K = closure([(1, 2, 3)])
f = {1: 0, 2: 1, 3: 2}
extend(K, f)
Out[92]:
In [100]:
K = closure([(1, 2, 3), (2, 3, 4)])
f = {1: 0, 2: 1, 3: 2, 4: 0}
extend(K, f)
Out[100]:
In [101]:
K = closure([(1, 2, 3), (2, 3, 4)])
f = {1: 0, 2: 1, 3: 2, 4: 3}
extend(K, f)
Out[101]:
In [ ]: