Commit 85a520fd authored by Martin Reinecke's avatar Martin Reinecke
Browse files

remove NFFT

parent 35dc4f9b
Pipeline #47410 passed with stages
in 16 minutes and 28 seconds
......@@ -10,7 +10,7 @@ RUN apt-get update && apt-get install -y \
# Testing dependencies
python3-pytest-cov jupyter \
# Optional NIFTy dependencies
libfftw3-dev python3-mpi4py python3-matplotlib python3-pynfft \
libfftw3-dev python3-mpi4py python3-matplotlib \
# more optional NIFTy dependencies
&& pip3 install pyfftw \
&& pip3 install git+ \
......@@ -45,18 +45,6 @@ for ii in range(10, 23):
b0s.append(t2 - t1)
c0s.append(t3 - t2)
t0 = time()
op = ift.NFFT(uvspace, uv)
t1 = time()
t2 = time()
t3 = time()
a1s.append(t1 - t0)
b1s.append(t2 - t1)
c1s.append(t3 - t2)
print('Measure rest operator')
sc = ift.StatCalculator()
op = GM.getRest().adjoint
......@@ -68,7 +56,6 @@ t_fft = sc.mean
print('FFT shape', res.shape)
plt.scatter(N0s, a0s, label='Gridder mr')
plt.scatter(N1s, a1s, marker='^', label='NFFT')
# no idea why this is necessary, but if it is omitted, the range is wrong
plt.ylim(min(a0s+a1s), max(a0s+a1s))
......@@ -79,9 +66,7 @@ plt.savefig('bench0.png')
plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
plt.scatter(N1s, b1s, color='r', marker='^', label='NFFT times')
plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
plt.scatter(N1s, c1s, color='r', label='NFFT adjoint times')
plt.axhline(sc.mean, label='FFT')
plt.axhline(sc.mean + np.sqrt(sc.var))
plt.axhline(sc.mean - np.sqrt(sc.var))
......@@ -86,7 +86,6 @@ from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.correlated_fields import CorrelatedField, MfCorrelatedField
from .library.adjust_variances import (make_adjust_variances_hamiltonian,
from .library.nfft import NFFT
from .library.gridder import GridderMaker
from . import extra
# 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
# 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 <>.
# Copyright(C) 2018-2019 Max-Planck-Society
# Resolve is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import nifty5 as ift
class NFFT(ift.LinearOperator):
"""Performs a non-equidistant Fourier transform, i.e. a Fourier transform
followed by a degridding operation.
domain : RGSpace
Domain of the operator. It has to be two-dimensional and have shape
`(2N, 2N)`. The coordinates of the lower left pixel of the dirty image
are `(-N,-N)`, and of the upper right pixel `(N-1,N-1)`.
uv : numpy.ndarray
2D numpy array of type float64 and shape (M,2), where M is the number
of measurements. uv[i,0] and uv[i,1] contain the u and v coordinates
of measurement #i, respectively. All coordinates must lie in the range
`[-0.5; 0,5[`.
def __init__(self, domain, uv):
from pynfft.nfft import NFFT
npix = domain.shape[0]
assert npix == domain.shape[1]
assert len(domain.shape) == 2
assert type(npix) == int, "npix must be integer"
assert npix > 0 and (
npix % 2) == 0, "npix must be an even, positive integer"
assert isinstance(uv, np.ndarray), "uv must be a Numpy array"
assert uv.dtype == np.float64, "uv must be an array of float64"
assert uv.ndim == 2, "uv must be a 2D array"
assert uv.shape[0] > 0, "at least one point needed"
assert uv.shape[1] == 2, "the second dimension of uv must be 2"
assert np.all(uv >= -0.5) and np.all(uv <= 0.5),\
"all coordinates must lie between -0.5 and 0.5"
self._domain = ift.DomainTuple.make(domain)
self._target = ift.DomainTuple.make(
self._capability = self.TIMES | self.ADJOINT_TIMES
self.npt = uv.shape[0]
self.plan = NFFT(self.domain.shape, self.npt, m=6)
self.plan.x = uv
def apply(self, x, mode):
self._check_input(x, mode)
if mode == self.TIMES:
self.plan.f_hat = x.to_global_data()
res = self.plan.trafo().copy()
self.plan.f = x.to_global_data()
res = self.plan.adjoint().copy()
return ift.Field.from_global_data(self._tgt(mode), res)
......@@ -293,10 +293,3 @@ def testValueInserter(sp, seed):
ind.append(np.random.randint(0, ss-1))
op = ift.ValueInserter(sp, ind)
def testNFFT():
dom = ift.RGSpace(2*(16,))
uv = np.array([[.2, .4], [-.22, .452]])
op = ift.NFFT(dom, uv)
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