Skip to content

Commit 5059ae4

Browse files
committed
First attempt at testElement for FFB2D
1 parent 7faf923 commit 5059ae4

File tree

1 file changed

+47
-0
lines changed

1 file changed

+47
-0
lines changed

tests/test_FFBbasis2D.py

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
from unittest import TestCase
44

55
import numpy as np
6+
from scipy.special import jv
67

78
from aspire.basis import FFBBasis2D
89
from aspire.image import Image
@@ -263,6 +264,52 @@ def testFFBBasis2DExpand(self):
263264
)
264265
)
265266

267+
def testElement(self):
268+
ell = 1
269+
sgn = -1
270+
k = 2
271+
272+
indices = self.basis.indices()
273+
ells = indices["ells"]
274+
sgns = indices["sgns"]
275+
ks = indices["ks"]
276+
277+
g2d = grid_2d(self.L, dtype=self.dtype)
278+
mask = g2d["r"] < 1
279+
280+
r0 = self.basis.r0[k, ell]
281+
282+
# TODO: Figure out where these factors of 1 / 2 are coming from.
283+
r = g2d["r"] * self.L / 4
284+
285+
im = np.zeros((self.L, self.L), dtype=self.dtype)
286+
im[mask] = (
287+
(-1) ** k
288+
* np.sqrt(np.pi)
289+
* r0
290+
* jv(ell, 2 * np.pi * r[mask])
291+
/ ((2 * np.pi * r[mask]) ** 2 - r0**2)
292+
)
293+
294+
if sgn == 1:
295+
im *= np.sqrt(2) * np.cos(ell * g2d["phi"])
296+
else:
297+
im *= np.sqrt(2) * np.sin(ell * g2d["phi"])
298+
299+
coef_ref = np.zeros(self.basis.count, dtype=self.dtype)
300+
coef_ref[(ells == ell) & (sgns == sgn) & (ks == k)] = 1
301+
302+
im_ref = self.basis.evaluate(coef_ref).asnumpy()[0]
303+
304+
coef = self.basis.expand(im)
305+
306+
# NOTE: These tolerances are expected to be rather loose since the
307+
# above expression for `im` is derived from the analytical formulation
308+
# (eq. 6 in Zhao and Singer, 2013) and does not take into account
309+
# discretization and other approximations.
310+
self.assertTrue(np.allclose(im, im_ref, atol=1e-1))
311+
self.assertTrue(np.allclose(coef, coef_ref, atol=1e-1))
312+
266313
def testRotate(self):
267314
# Now low res (8x8) had problems;
268315
# better with odd (7x7), but still not good.

0 commit comments

Comments
 (0)