bindings/python/dionysus/circular/__init__.py
 author Dmitriy Morozov Sun, 24 Nov 2019 13:35:45 -0800 changeset 298 02e680f6d987 parent 267 2f02384a4d9b permissions -rw-r--r--
```
def smooth(filtration, cocycle):
from    cvxopt          import spmatrix, matrix
from    cvxopt.blas     import copy
from    lsqr            import lsqr

coefficient = []
coface_indices = []
face_indices = []
for i,s in enumerate(filtration):
if s.dimension() > 2: continue

c = 1
for sb in s.boundary:
j = filtration(sb)
coefficient.append(c)
coface_indices.append(i)
face_indices.append(j)
c *= -1

# Cocycle can be larger than D; we implicitly project it down
cocycle_max = max(zz[1] for zz in cocycle)

# D is a coboundary matrix
dimension = max(max(coface_indices), max(face_indices), cocycle_max) + 1
D = spmatrix(coefficient, coface_indices, face_indices, (dimension, dimension))

z = spmatrix([zz[0] for zz in cocycle],
[zz[1] for zz in cocycle],
[0     for zz in cocycle], (dimension, 1))

v1 = D * z
if bool(D*D):
raise Exception('D^2 is not 0')
if bool(v1):
raise Exception('Expected a cocycle as input')
z = matrix(z)

def Dfun(x,y,trans = 'N'):
if trans == 'N':
copy(D * x, y)
elif trans == 'T':
copy(D.T * x, y)
else:
assert False, "Unexpected trans parameter"

tol = 1e-10
show = False
maxit = None
solution = lsqr(Dfun, matrix(z), show = show, atol = tol, btol = tol, itnlim = maxit)

z_smooth = z - D*solution[0]

# print sum(z_smooth**2)
# assert sum((D*z_smooth)**2) < tol and sum((D.T*z_smooth)**2) < tol, "Expected a harmonic cocycle"
if not (sum((D*z_smooth)**2) < tol and sum((D.T*z_smooth)**2) < tol):
raise Exception("Expected a harmonic cocycle: %f %f" % (sum((D*z_smooth)**2), sum((D.T*z_smooth)**2)))

values = []
vertices = ((i,s) for (i,s) in enumerate(filtration) if s.dimension() == 0)
for i,s in vertices:
v = [v for v in s.vertices][0]
if v >= len(values):
values.extend((None for i in xrange(len(values), v+1)))
values[v] = solution[0][i]

return values
```