diff --git a/Dockerfile b/Dockerfile index b121048922d1beb08d418c54585f21150a3ee5c5..210872a95720c723a02a35cbaf07c881be84842d 100644 --- a/Dockerfile +++ b/Dockerfile @@ -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+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 8658f8f32fafed19cd07ce5b80d5e49956df4d88..5dd68c1cd157605a28445edec33500db46dfaeba 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -8,7 +8,6 @@ import nifty5 as ift np.random.seed(40) N0s, a0s, b0s, c0s = [], [], [], [] -N1s, a1s, b1s, c1s = [], [], [], [] for ii in range(10, 23): nu = 1024 @@ -44,18 +43,6 @@ for ii in range(10, 23): b0s.append(t2 - t1) c0s.append(t3 - t2) - t0 = time() - op = ift.NFFT(uvspace, uv) - t1 = time() - op(img).to_global_data() - t2 = time() - op.adjoint(vis).to_global_data() - t3 = time() - N1s.append(N) - a1s.append(t1 - t0) - b1s.append(t2 - t1) - c1s.append(t3 - t2) - print('Measure rest operator') sc = ift.StatCalculator() op = GM.getRest().adjoint @@ -67,10 +54,9 @@ t_fft = sc.mean print('FFT shape', res.shape) plt.scatter(N0s, a0s, label='Gridder mr') -plt.scatter(N1s, a1s, marker='^', label='NFFT') plt.legend() # no idea why this is necessary, but if it is omitted, the range is wrong -plt.ylim(min(a0s+a1s), max(a0s+a1s)) +plt.ylim(min(a0s), max(a0s)) plt.ylabel('time [s]') plt.title('Initialization') plt.loglog() @@ -78,9 +64,7 @@ plt.savefig('bench0.png') plt.close() 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)) diff --git a/nifty5/__init__.py b/nifty5/__init__.py index dd1ccb298edf7886a01c37635a3873a78dc8f1c2..eee17e5bbc2a618113f5da0a202d09581e3a270b 100644 --- a/nifty5/__init__.py +++ b/nifty5/__init__.py @@ -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, do_adjust_variances) -from .library.nfft import NFFT from .library.gridder import GridderMaker from . import extra diff --git a/nifty5/domains/gl_space.py b/nifty5/domains/gl_space.py index d6091d5388129f14515d1f5c2e0ce4659d688e25..85a0dba9927bfb62bfad30766076ebbbb641b016 100644 --- a/nifty5/domains/gl_space.py +++ b/nifty5/domains/gl_space.py @@ -103,7 +103,9 @@ class GLSpace(StructuredDomain): The partner domain """ from ..domains.lm_space import LMSpace - return LMSpace(lmax=self._nlat-1, mmax=self._nlon//2) + mmax = self._nlon//2 + lmax = max(mmax, self._nlat-1) + return LMSpace(lmax=lmax, mmax=mmax) def check_codomain(self, codomain): """Raises `TypeError` if `codomain` is not a matching partner domain diff --git a/nifty5/library/nfft.py b/nifty5/library/nfft.py deleted file mode 100644 index c96866d6ee948c0fbda23a7720c53b792059e5d9..0000000000000000000000000000000000000000 --- a/nifty5/library/nfft.py +++ /dev/null @@ -1,73 +0,0 @@ -# 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 . -# -# 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. - - Parameters - ---------- - 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( - ift.UnstructuredDomain(uv.shape[0])) - 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 - self.plan.precompute() - - 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() - else: - self.plan.f = x.to_global_data() - res = self.plan.adjoint().copy() - return ift.Field.from_global_data(self._tgt(mode), res) diff --git a/nifty5/operators/energy_operators.py b/nifty5/operators/energy_operators.py index 70f79de170baa88c79e2831269ccf55a55061a43..17a97d8040ee44e643c73c59219a4370a6281987 100644 --- a/nifty5/operators/energy_operators.py +++ b/nifty5/operators/energy_operators.py @@ -191,7 +191,12 @@ class PoissonianEnergy(EnergyOperator): def apply(self, x): self._check_input(x) - res = x.sum() - x.log().vdot(self._d) + res = x.sum() + tmp = res.val.local_data if isinstance(res, Linearization) else res + # if we have no infinity here, we can continue with the calculation; + # otherwise we know that the result must also be infinity + if not np.isinf(tmp): + res = res - x.log().vdot(self._d) if not isinstance(x, Linearization): return Field.scalar(res) if not x.want_metric: diff --git a/nifty5/operators/simple_linear_operators.py b/nifty5/operators/simple_linear_operators.py index 096ebb2a24fb03d0a99bff835d79a9dd63a50b00..1f3d11d2d5c978bfe66033f72a8727e4c190e8a6 100644 --- a/nifty5/operators/simple_linear_operators.py +++ b/nifty5/operators/simple_linear_operators.py @@ -147,7 +147,16 @@ class FieldAdapter(LinearOperator): return MultiField(self._tgt(mode), (x,)) def __repr__(self): - return 'FieldAdapter' + s = 'FieldAdapter' + dom = isinstance(self._domain, MultiDomain) + tgt = isinstance(self._target, MultiDomain) + if dom and tgt: + s += ' {} <- {}'.format(self._target.keys(), self._domain.keys()) + elif dom: + s += ' <- {}'.format(self._domain.keys()) + elif tgt: + s += ' {} <-'.format(self._target.keys()) + return s class _SlowFieldAdapter(LinearOperator): diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py index 97e020763c1803636d9d1f7ddf0ee1ee69e7b97b..56d2534ae535598649e02037a313dce4f9e1e217 100644 --- a/test/test_operators/test_adjoint.py +++ b/test/test_operators/test_adjoint.py @@ -293,10 +293,3 @@ def testValueInserter(sp, seed): ind.append(np.random.randint(0, ss-1)) op = ift.ValueInserter(sp, ind) ift.extra.consistency_check(op) - - -def testNFFT(): - dom = ift.RGSpace(2*(16,)) - uv = np.array([[.2, .4], [-.22, .452]]) - op = ift.NFFT(dom, uv) - ift.extra.consistency_check(op)