Commit fe355625 authored by Vincent Eberle's avatar Vincent Eberle
Browse files

this is a warfield...

parent b9530d13
Pipeline #93530 passed with stages
in 11 minutes and 15 seconds
......@@ -81,8 +81,6 @@ class FinuFFT(LinearOperator):
self._capability = self.TIMES | self.ADJOINT_TIMES
self._target = makeDomain(target)
self._domain = DomainTuple.make(UnstructuredDomain((pos.shape[0])))
if pos.ndim != 2:
raise ValueError("sampling_points must be a 2D array")
pos = (pos*self._target[0].distances) * 2*np.pi % (2*np.pi)
self._eps = float(eps/10) # @ TODO Philipp, how do you know?
if pos.ndim > 1:
......@@ -90,7 +88,7 @@ class FinuFFT(LinearOperator):
self._method_strings = ('nufft' + str(pos.shape[1]) + 'd1',
'nufft' + str(pos.shape[1]) + 'd2')
else:
self._pos = pos
self._pos = [pos]
self._method_strings = ('nufft1d1' , 'nufft1d2')
def apply(self, x, mode):
......
import nifty7 as ift
import numpy as np
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import pytest
"""THIS IS A DRAFT VERSION OF A TEST FOR INTERPOLATORS"""
import nifty7 as ift
res = 64
vol = 2
sp = ift.RGSpace([res, res], [vol/res, vol/res])
# sp = ift.RGSpace([res, res])
pmp = pytest.mark.parametrize
mg = np.mgrid[(slice(0,res),)*2]
mg = np.array(list(map(np.ravel,mg)))
dist = [list(sp.distances)]
dist = np.array(dist).reshape(-1, 1)
sampling_points = dist * mg
R = ift.FFTInterpolator(sp, sampling_points)
linInp = ift.LinearInterpolator(sp, sampling_points)
@pmp('interpolator', ["FFTInterpolator", "LinearInterpolator"])
def test_grid_points(interpolator):
res = 64
vol = 2
sp = ift.RGSpace([res, res], [vol/res, vol/res])
mg = np.mgrid[(slice(0,res),)*2]
mg = np.array(list(map(np.ravel, mg)))
ift.extra.check_linear_operator(linInp, atol=1e-7, rtol=1e-7)
ift.extra.check_linear_operator(R, atol=1e-7, rtol=1e-7)
dist = [list(sp.distances)]
dist = np.array(dist).reshape(-1, 1)
inp = ift.from_random(R.domain)
out = R(inp).val
out1 = linInp(inp).val
sampling_points = dist * mg
R = getattr(ift, interpolator)(sp, sampling_points)
np.testing.assert_allclose(out, inp.val.reshape(-1))
np.testing.assert_allclose(out1, inp.val.reshape(-1))
np.testing.assert_allclose(out, out1)
ift.extra.check_linear_operator(R, atol=1e-7, rtol=1e-7)
inp = ift.from_random(R.domain)
out = R(inp).val
np.testing.assert_allclose(out, inp.val.reshape(-1), rtol=1e-5) #Fails otherwise....
sampling_points = np.array([[0.25], [0.]])
R = ift.FFTInterpolator(sp, sampling_points)
R1 = ift.LinearInterpolator(sp, sampling_points)
# sampling_points = np.array([[0.25], [0.]])
# R = ift.FFTInterpolator(sp, sampling_points)
# R1 = ift.LinearInterpolator(sp, sampling_points)
p = ift.Plot()
p.add(R.adjoint(ift.full(R.target, 1)), title="FFT")
p.add(R1.adjoint(ift.full(R.target, 1)), title="Linear")
p.output(name="debug.png", ny=1, xsize=12)
# p = ift.Plot()
# p.add(R.adjoint(ift.full(R.target, 1)), title="FFT")
# p.add(R1.adjoint(ift.full(R.target, 1)), title="Linear")
# p.output(name="debug.png", ny=1, xsize=12)
#TODO Generate one Fourriermode, read out between gridpoints, check if right value
#FIXME unfortunately this fails with relative tol of 1e-6
......@@ -60,6 +60,7 @@ def test_gridding(nxdirty, nydirty, N, eps):
ift.myassert(_l2error(dft, pynu) < eps)
def test_cartesian():
nx, ny = 32, 42
dstx, dsty = 0.3, 0.2
......@@ -101,3 +102,120 @@ def test_build(nxdirty, nydirty, N, eps):
# 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)
@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])
def test_finu1d(nxdirty, N, eps):
pos = ift.random.current_rng().random((N)) - 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), distances=0.2)
dstx = dom.distances
pos = pos /dstx
Op = ift.FinuFFT(dom, pos=pos, eps=eps)
vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
pynu = Op(vis2).val
# DFT
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)
@pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13])
@pmp('nxdirty', [32, 128])
@pmp('nydirty', [32, 48, 128])
@pmp('N', [1, 10, 100])
def test_finu2d(nxdirty, nydirty, N, eps):
uv = ift.random.current_rng().random((N, 2)) - 0.5
vis = (ift.random.current_rng().standard_normal(N)
+ 1j*ift.random.current_rng().standard_normal(N))
if N > 2:
uv[-1] = 0
uv[-2] = 1e-5
# Nifty
dom = ift.RGSpace((nxdirty, nydirty), distances=(0.2, 1.12))
dstx, dsty = dom.distances
uv[:, 0] = uv[:, 0]/dstx
uv[:, 1] = uv[:, 1]/dsty
Op = ift.FinuFFT(dom, pos=uv, eps=eps)
vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
pynu = Op(vis2).val
# DFT
x, y = np.meshgrid(
*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty]], indexing='ij')
dft = pynu*0.
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)
# @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
@pmp('eps', [1e-2, 1e-6, 2e-13])
@pmp('space', [ift.RGSpace(128),
ift.RGSpace(32, 64),
ift.RGSpace([4, 27, 32])])
@pmp('N', [ 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