Commit c7e93629 authored by Philipp Arras's avatar Philipp Arras
Browse files

Tweaks

parent 1d6e850e
Pipeline #93689 failed with stages
in 6 minutes and 15 seconds
......@@ -84,8 +84,8 @@ class FinuFFT(LinearOperator):
self._eps = float(eps)
dst = np.array(self._target[0].distances)
pos = (2*np.pi*pos*dst) % (2*np.pi)
self._eps = float(eps/10)
if pos.ndim > 1:
self._eps = float(eps)
if pos.shape[0] > 1:
self._pos = [pos[:, k] for k in range(pos.shape[1])]
s = 'nufft' + str(pos.shape[1]) + 'd'
else:
......@@ -99,7 +99,6 @@ class FinuFFT(LinearOperator):
if mode == self.TIMES:
res = self._f(*self._pos, c=x.val_rw(),
n_modes=self._target[0].shape, eps=self._eps).real
# TODO is this .real needed?
if mode == self.ADJOINT_TIMES:
else:
res = self._fadj(*self._pos, f=x.val, eps=self._eps)
return makeField(self._tgt(mode), res)
......@@ -59,8 +59,6 @@ def test_gridding(nxdirty, nydirty, N, eps):
vis[i]*np.exp(2j*np.pi*(x*uv[i, 0]*dstx + y*uv[i, 1]*dsty))).real
ift.myassert(_l2error(dft, pynu) < eps)
def test_cartesian():
nx, ny = 32, 42
dstx, dsty = 0.3, 0.2
......@@ -86,7 +84,6 @@ def test_cartesian():
np.testing.assert_allclose(res, res1*vol)
@pmp('eps', [1e-2, 1e-6, 2e-13])
@pmp('nxdirty', [32, 128])
@pmp('nydirty', [32, 48, 128])
......@@ -105,7 +102,7 @@ def test_build(nxdirty, nydirty, N, eps):
@pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13])
@pmp('nxdirty', [32, 128])
@pmp('N', [1, 10 , 100])
@pmp('N', [1, 10, 100])
def test_finu1d(nxdirty, N, eps):
pos = ift.random.current_rng().random((N)) - 0.5
vis = (ift.random.current_rng().standard_normal(N)
......@@ -117,7 +114,7 @@ def test_finu1d(nxdirty, N, eps):
# Nifty
dom = ift.RGSpace((nxdirty), distances=0.2)
dstx = dom.distances
pos = pos /dstx
pos = pos / dstx
Op = ift.FinuFFT(dom, pos=pos, eps=eps)
vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
pynu = Op(vis2).val
......@@ -125,12 +122,9 @@ def test_finu1d(nxdirty, N, eps):
x = -nxdirty/2 + np.arange(nxdirty)
dft = pynu*0
print(x)
print(dstx)
for i in range(N):
dft += (
vis[i]*np.exp(2j*np.pi*(x*pos[i]*dstx))).real
ift.myassert(_l2error(dft, pynu) < eps)
dft += (vis[i]*np.exp(2j*np.pi*(x*pos[i]*dstx))).real
ift.myassert(_l2error(dft, pynu) < eps*10)
@pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13])
@pmp('nxdirty', [32, 128])
......@@ -160,62 +154,48 @@ def test_finu2d(nxdirty, nydirty, N, eps):
for i in range(N):
dft += (
vis[i]*np.exp(2j*np.pi*(x*uv[i, 0]*dstx + y*uv[i, 1]*dsty))).real
ift.myassert(_l2error(dft, pynu) < eps)
ift.myassert(_l2error(dft, pynu) < eps*10)
@pmp('eps', [1e-2, 1e-4, 1e-7, 1e-11])
@pmp('nxdirty', [32, 128])
@pmp('nydirty', [32, 48])
@pmp('nzdirty', [32, 54])
@pmp('N', [1, 10])
def test_finu3d(nxdirty, nydirty, nzdirty, N, eps):
pos = ift.random.current_rng().random((N, 3)) - 0.5
vis = (ift.random.current_rng().standard_normal(N)
+ 1j*ift.random.current_rng().standard_normal(N))
# Nifty
dom = ift.RGSpace((nxdirty, nydirty, nzdirty), distances=(0.2, 1.12, 0.7))
dstx, dsty, dstz = dom.distances
pos[:, 0] = pos[:, 0]/dstx
pos[:, 1] = pos[:, 1]/dsty
pos[:, 2] = pos[:, 2]/dstz
Op = ift.FinuFFT(dom, pos=pos, eps=eps)
vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
# @pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 2e-13])
# @pmp('nxdirty', [32, 128])
# @pmp('nydirty', [32, 48, 128])
# @pmp('nzdirty', [32, 54]) #FIXME crashes with 55
# @pmp('N', [1, 10, 100])
# def test_finu3d(nxdirty, nydirty, nzdirty, N, eps):
# pos = ift.random.current_rng().random((N, 3)) - 0.5
# vis = (ift.random.current_rng().standard_normal(N)
# + 1j*ift.random.current_rng().standard_normal(N))
# # if N > 2:
# # pos[-1] = 0
# # pos[-2] = 1e-5
# # Nifty
# dom = ift.RGSpace((nxdirty, nydirty, nzdirty), distances=(0.2, 1.12, 0.7))
# dstx, dsty, dstz = dom.distances
# pos[:, 0] = pos[:, 0]/dstx
# pos[:, 1] = pos[:, 1]/dsty
# pos[:, 2] = pos[:, 2]/dstz
# Op = ift.FinuFFT(dom, pos=pos, eps=eps)
# vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
# pynu = Op(vis2).val
# # DFT
# x, y, z= np.meshgrid(
# *[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty, nzdirty]], indexing='ij')
# dft = pynu*0.
# for i in range(N):
# dft += (
# vis[i]*np.exp(2j*np.pi*(x*pos[i, 0]*dstx + y*pos[i, 1]*dsty + z*pos[i, 2]*dstz))).real
# ift.myassert(_l2error(dft, pynu) < eps)
#NOTE Could we switch to FFT here?
#NOTE simplify to one test
#NOTE test cartesian
pynu = Op(vis2).val
# DFT
x, y, z = np.meshgrid(
*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty, nzdirty]], indexing='ij')
dft = pynu*0.
for i in range(N):
dft += (
vis[i]*np.exp(2j*np.pi*(x*pos[i, 0]*dstx + y*pos[i, 1]*dsty + z*pos[i, 2]*dstz))).real
ift.myassert(_l2error(dft, pynu) < eps*10)
@pmp('eps', [1e-2, 1e-6, 2e-13])
@pmp('space', [ift.RGSpace(128),
ift.RGSpace(32, 64),
ift.RGSpace([32, 64]),
ift.RGSpace([4, 27, 32])])
@pmp('N', [ 10, 100])
@pmp('N', [1, 10, 100])
def test_build_finufft(space, N, eps):
# if len(space.shape) > 1:
pos = ift.random.current_rng().random((N, len(space.shape))) - 0.5
# else:
# pos = ift.random.current_rng().random(N) - 0.5
RF = ift.FinuFFT(space, pos=pos, eps=eps)
# Consistency checks
flt = np.float64
cmplx = np.complex128
# We set rtol=eps here, because the gridder operator only guarantees
# adjointness to this accuracy.
ift.extra.check_linear_operator(RF, cmplx, flt, only_r_linear=True, rtol=eps)
#FIXME somethings wrong with N = 1
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment