-
Notifications
You must be signed in to change notification settings - Fork 0
/
Laplacian.py
69 lines (56 loc) · 2.06 KB
/
Laplacian.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
""" Approximate Laplace matrix via heat kernel. """
import numexpr as ne
import numpy as np
def computeLaplaceMatrix(sqdist, t, L):
""" Compute heat approximation to Laplacian matrix using logarithms. """
# FIXME Is the scaling ok in here?
lt = np.log(2.0 / t)
ne.evaluate('sqdist / (-2.0 * t)', out=L)
# sqdist is nonnegative, with 0 on the diagonal
# so the largest element of each row of L is 0
# no logsumexp needed
density = ne.evaluate("sum(exp(L), axis=1)")[:, None]
ne.evaluate("log(density)", out=density)
# sum in rows must be 1, except for 2/t factor
ne.evaluate('exp(L - density + lt)', out=L)
# fix diagonal to account for -f(x)?
L[np.diag_indices(len(L))] -= 2.0 / t
Laplacian = computeLaplaceMatrix
#
# tests based on computeLaplaceMatrix2
#
import unittest
class LaplaceTests (unittest.TestCase):
""" Correctness and speed tests. """
def row_sums(self, f):
""" Test if rows add up to 0. """
import data
threshold = 1E-10
print
for d in data.tests('small'):
L = np.random.rand(*d.shape)
f(d, 0.1, L)
error = np.max(np.abs(np.sum(L, axis=1)))
print "Absolute error: ", error
self.assertLess(error, threshold)
def test_rows_numpy(self):
""" Check sums of rows with numpy. """
self.row_sums(computeLaplaceMatrix)
def speed(self, f):
""" Test speed on larger data sets. """
import data
from tools import test_speed
d = data.closefarsimplices(500, 0.1, 5)[0]
print "\nPoints: 1000"
L = np.random.rand(*d.shape)
test_speed(f, d, 0.1, L)
d = data.closefarsimplices(1000, 0.1, 5)[0]
L = np.random.rand(*d.shape)
print "Points: 2000"
test_speed(f, d, 0.1, L)
def test_speed_numpy(self):
""" Speed with numpy. """
self.speed(computeLaplaceMatrix)
if __name__ == "__main__":
suite = unittest.TestLoader().loadTestsFromTestCase(LaplaceTests)
unittest.TextTestRunner(verbosity=2).run(suite)