Commit a7c24fd2 authored by Reimar Heinrich Leike's avatar Reimar Heinrich Leike

Merge branch 'NIFTy_5' into parallizing_mirrored_samples

parents 017b391f e50e65f4
......@@ -10,11 +10,11 @@ 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 \
python3-mpi4py python3-matplotlib \
# more optional NIFTy dependencies
&& pip3 install pyfftw \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft.git \
&& pip3 install jupyter \
&& rm -rf /var/lib/apt/lists/*
......
......@@ -47,9 +47,9 @@ Installation
- [Python 3](https://www.python.org/) (3.5.x or later)
- [SciPy](https://www.scipy.org/)
- [pypocketfft](https://gitlab.mpcdf.mpg.de/mtr/pypocketfft)
Optional dependencies:
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW) for faster Fourier transforms
- [pyHealpix](https://gitlab.mpcdf.mpg.de/ift/pyHealpix) (for harmonic
transforms involving domains on the sphere)
- [nifty_gridder](https://gitlab.mpcdf.mpg.de/ift/nifty_gridder) (for radio
......@@ -73,28 +73,12 @@ NIFTy5 and its mandatory dependencies can be installed via:
sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_5
pip3 install --user git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft
Plotting support is added via:
sudo apt-get install python3-matplotlib
NIFTy uses Numpy's FFT implementation by default. For large problems FFTW may be
used because of its higher performance. It can be installed via:
sudo apt-get install libfftw3-dev
pip3 install --user pyfftw
To enable FFTW usage in NIFTy, call
nifty5.fft.enable_fftw()
at the beginning of your code.
(Note: If you encounter problems related to `pyFFTW`, make sure that you are
using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
shipping an incorrectly configured `pyFFTW` package, which does not cooperate
with the installed `FFTW3` libraries.)
Support for spherical harmonic transforms is added via:
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
......
......@@ -5,13 +5,11 @@ import numpy as np
import nifty5 as ift
ift.fft.enable_fftw()
np.random.seed(40)
N0s, a0s, b0s, c0s = [], [], [], []
N1s, a1s, b1s, c1s = [], [], [], []
for ii in range(10, 23):
for ii in range(10, 26):
nu = 1024
nv = 1024
N = int(2**ii)
......@@ -29,34 +27,20 @@ for ii in range(10, 23):
img = ift.from_global_data(uvspace, img)
t0 = time()
GM = ift.GridderMaker(uvspace, eps=1e-7)
idx = GM.getReordering(uv)
uv = uv[idx]
vis = vis[idx]
GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
vis = ift.from_global_data(visspace, vis)
op = GM.getFull(uv).adjoint
op = GM.getFull().adjoint
t1 = time()
op(img).to_global_data()
t2 = time()
op.adjoint(vis).to_global_data()
t3 = time()
print(t2-t1, t3-t2)
N0s.append(N)
a0s.append(t1 - t0)
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
......@@ -68,10 +52,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()
......@@ -79,9 +62,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))
......
......@@ -109,7 +109,8 @@ if __name__ == '__main__':
minimizer = ift.NewtonCG(ic_newton)
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
likelihood = ift.GaussianEnergy(mean=data,
inverse_covariance=N.inverse)(signal_response)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
initial_mean = ift.MultiField.full(H.domain, 0.)
......
......@@ -2,6 +2,12 @@ NIFTy-related publications
==========================
::
@article{asclnifty5,
title={NIFTy5: Numerical Information Field Theory v5},
author={Arras, Philipp and Baltac, Mihai and Ensslin, Torsten A and Frank, Philipp and Hutschenreuter, Sebastian and Knollmueller, Jakob and Leike, Reimar and Newrzella, Max-Niklas and Platz, Lukas and Reinecke, Martin and others},
journal={Astrophysics Source Code Library},
year={2019}
}
@software{nifty,
author = {{Martin Reinecke, Theo Steininger, Marco Selig}},
......
......@@ -9,28 +9,12 @@ NIFTy5 and its mandatory dependencies can be installed via::
sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_5
pip3 install --user git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft
Plotting support is added via::
sudo apt-get install python3-matplotlib
NIFTy uses Numpy's FFT implementation by default. For large problems FFTW may be
used because of its higher performance. It can be installed via::
sudo apt-get install libfftw3-dev
pip3 install --user pyfftw
To enable FFTW usage in NIFTy, call::
nifty5.fft.enable_fftw()
at the beginning of your code.
(Note: If you encounter problems related to `pyFFTW`, make sure that you are
using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
shipping an incorrectly configured `pyFFTW` package, which does not cooperate
with the installed `FFTW3` libraries.)
Support for spherical harmonic transforms is added via::
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
......
......@@ -45,7 +45,8 @@ from .operators.block_diagonal_operator import BlockDiagonalOperator
from .operators.outer_product_operator import OuterProduct
from .operators.simple_linear_operators import (
VdotOperator, ConjugationOperator, Realizer,
FieldAdapter, ducktape, GeometryRemover, NullOperator)
FieldAdapter, ducktape, GeometryRemover, NullOperator,
MatrixProductOperator)
from .operators.value_inserter import ValueInserter
from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
......@@ -86,7 +87,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
......
......@@ -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
......
......@@ -80,8 +80,8 @@ class LogRGSpace(StructuredDomain):
return np.array(self._t_0)
def __repr__(self):
return ("LogRGSpace(shape={}, harmonic={})".format(
self.shape, self.harmonic))
return ("LogRGSpace(shape={}, bindistances={}, t_0={}, harmonic={})".format(
self.shape, self.bindistances, self.t_0, self.harmonic))
def get_default_codomain(self):
"""Returns a :class:`LogRGSpace` object representing the (position or
......
......@@ -17,8 +17,10 @@
import numpy as np
from .domain_tuple import DomainTuple
from .field import Field
from .linearization import Linearization
from .multi_domain import MultiDomain
from .operators.linear_operator import LinearOperator
from .sugar import from_random
......@@ -70,12 +72,20 @@ def _full_implementation(op, domain_dtype, target_dtype, atol, rtol,
def _check_linearity(op, domain_dtype, atol, rtol):
fld1 = from_random("normal", op.domain, dtype=domain_dtype)
fld2 = from_random("normal", op.domain, dtype=domain_dtype)
alpha = np.random.random()
alpha = np.random.random() # FIXME: this can break badly with MPI!
val1 = op(alpha*fld1+fld2)
val2 = alpha*op(fld1)+op(fld2)
_assert_allclose(val1, val2, atol=atol, rtol=rtol)
def _domain_check(op):
for dd in [op.domain, op.target]:
if not isinstance(dd, (DomainTuple, MultiDomain)):
raise TypeError(
'The domain and the target of an operator need to',
'be instances of either DomainTuple or MultiDomain.')
def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
atol=0, rtol=1e-7, only_r_linear=False):
"""
......@@ -109,6 +119,7 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
"""
if not isinstance(op, LinearOperator):
raise TypeError('This test tests only linear operators.')
_domain_check(op)
_check_linearity(op, domain_dtype, atol, rtol)
_full_implementation(op, domain_dtype, target_dtype, atol, rtol,
only_r_linear)
......@@ -162,6 +173,7 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100):
tol : float
Tolerance for the check.
"""
_domain_check(op)
for _ in range(ntries):
lin = op(Linearization.make_var(loc))
loc2, lin2 = _get_acceptable_location(op, loc, lin)
......
......@@ -15,152 +15,28 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from .utilities import iscomplextype
import numpy as np
import pypocketfft
_nthreads = 1
_use_fftw = False
_fftw_prepped = False
_fft_extra_args = {}
def nthreads():
return _nthreads
def enable_fftw():
global _use_fftw
_use_fftw = True
def disable_fftw():
global _use_fftw
_use_fftw = False
def _init_pyfftw():
global _fft_extra_args, _fftw_prepped
if not _fftw_prepped:
import pyfftw
from pyfftw.interfaces.numpy_fft import fftn, rfftn, ifftn
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(1000.)
# Optional extra arguments for the FFT calls
# if exact reproducibility is needed,
# set "planner_effort" to "FFTW_ESTIMATE"
import os
nthreads = int(os.getenv("OMP_NUM_THREADS", "1"))
_fft_extra_args = dict(planner_effort='FFTW_ESTIMATE',
threads=nthreads)
_fftw_prepped = True
def set_nthreads(nthr):
global _nthreads
_nthreads = nthr
def fftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import fftn
_init_pyfftw()
return fftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.fftn(a, axes=axes)
def rfftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import rfftn
_init_pyfftw()
return rfftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.rfftn(a, axes=axes)
return pypocketfft.c2c(a, axes=axes, nthreads=_nthreads)
def ifftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import ifftn
_init_pyfftw()
return ifftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.ifftn(a, axes=axes)
return pypocketfft.c2c(a, axes=axes, inorm=2, forward=False,
nthreads=_nthreads)
def hartley(a, axes=None):
# Check if the axes provided are valid given the shape
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes do not match array shape")
if iscomplextype(a.dtype):
raise TypeError("Hartley transform requires real-valued arrays.")
tmp = rfftn(a, axes=axes)
def _fill_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = (slice(None),)*lastaxis + (slice(0, ntmplast),)
np.add(tmp.real, tmp.imag, out=res[slice1])
def _fill_upper_half(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
slice1 = tuple(slice1)
slice2 = tuple(slice2)
np.subtract(tmp[slice2].real, tmp[slice2].imag, out=res[slice1])
for i, ax in enumerate(axes[:-1]):
dim1 = (slice(None),)*ax + (slice(0, 1),)
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half(tmp[dim1], res[dim1], axes2)
_fill_upper_half(tmp, res, axes)
return res
return _fill_array(tmp, np.empty_like(a), axes)
# Do a real-to-complex forward FFT and return the _full_ output array
def my_fftn_r2c(a, axes=None):
# Check if the axes provided are valid given the shape
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes do not match array shape")
if iscomplextype(a.dtype):
raise TypeError("Transform requires real-valued input arrays.")
tmp = rfftn(a, axes=axes)
def _fill_complex_array(tmp, res, axes):
if axes is None:
axes = tuple(range(tmp.ndim))
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
res[tuple(slice1)] = tmp
def _fill_upper_half_complex(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
# np.conjugate(tmp[slice2], out=res[slice1])
res[tuple(slice1)] = np.conjugate(tmp[tuple(slice2)])
for i, ax in enumerate(axes[:-1]):
dim1 = tuple([slice(None)]*ax + [slice(0, 1)])
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half_complex(tmp[dim1], res[dim1], axes2)
_fill_upper_half_complex(tmp, res, axes)
return res
return _fill_complex_array(tmp, np.empty_like(a, dtype=tmp.dtype), axes)
def my_fftn(a, axes=None):
return fftn(a, axes=axes)
return pypocketfft.genuine_hartley(a, axes=axes, nthreads=_nthreads)
......@@ -15,114 +15,96 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..fft import hartley
from ..operators.linear_operator import LinearOperator
from ..sugar import from_global_data, makeDomain
import numpy as np
class GridderMaker(object):
def __init__(self, domain, eps=1e-15):
domain = makeDomain(domain)
if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or
not len(domain.shape) == 2):
raise ValueError("need domain with exactly one 2D RGSpace")
nu, nv = domain.shape
if nu % 2 != 0 or nv % 2 != 0:
raise ValueError("dimensions must be even")
rat = 3 if eps < 1e-11 else 2
nu2, nv2 = rat*nu, rat*nv
nspread = int(-np.log(eps)/(np.pi*(rat-1)/(rat-.5)) + .5) + 1
nu2 = max([nu2, 2*nspread])
nv2 = max([nv2, 2*nspread])
r2lamb = rat*rat*nspread/(rat*(rat-.5))
oversampled_domain = RGSpace(
[nu2, nv2], distances=[1, 1], harmonic=False)
self._nspread = nspread
self._r2lamb = r2lamb
self._rest = _RestOperator(domain, oversampled_domain, r2lamb)
def getReordering(self, uv):
from nifty_gridder import peanoindex
nu2, nv2 = self._rest._domain.shape
return peanoindex(uv, nu2, nv2)
def getGridder(self, uv):
return RadioGridder(self._rest.domain, self._nspread, self._r2lamb, uv)
def __init__(self, dirty_domain, uv, eps=2e-13):
import nifty_gridder
dirty_domain = makeDomain(dirty_domain)
if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace)
or not len(dirty_domain.shape) == 2):
raise ValueError("need dirty_domain with exactly one 2D RGSpace")
if uv.ndim != 2:
raise ValueError("uv must be a 2D array")
if uv.shape[1] != 2:
raise ValueError("second dimension of uv must have length 2")
dstx, dsty = dirty_domain[0].distances
# wasteful hack to adjust to shape required by nifty_gridder
uvw = np.empty((uv.shape[0], 3), dtype=np.float64)
uvw[:, 0:2] = uv
uvw[:, 2] = 0.
# Scale uv such that 0<uv<=1 which is assumed by nifty_gridder
uvw[:, 0] = uvw[:, 0]*dstx
uvw[:, 1] = uvw[:, 1]*dsty
speedOfLight = 299792458.
bl = nifty_gridder.Baselines(uvw, np.array([speedOfLight]))
nxdirty, nydirty = dirty_domain.shape
gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.)
nu, nv = gconf.Nu(), gconf.Nv()
self._idx = nifty_gridder.getIndices(
bl, gconf, np.zeros((uv.shape[0], 1), dtype=np.bool))
self._bl = bl
du, dv = 1./(nu*dstx), 1./(nv*dsty)
grid_domain = RGSpace([nu, nv], distances=[du, dv], harmonic=True)
self._rest = _RestOperator(dirty_domain, grid_domain, gconf)
self._gridder = RadioGridder(grid_domain, bl, gconf, self._idx)
def getGridder(self):
return self._gridder
def getRest(self):
return self._rest
def getFull(self, uv):
return self.getRest() @ self.getGridder(uv)
def getFull(self):
return self.getRest() @ self._gridder
def ms2vis(self, x):
return self._bl.ms2vis(x, self._idx)
class _RestOperator(LinearOperator):
def __init__(self, domain, oversampled_domain, r2lamb):
self._domain = makeDomain(oversampled_domain)
self._target = domain
nu, nv = domain.shape
nu2, nv2 = oversampled_domain.shape
# compute deconvolution operator
rng = np.arange(nu)
k = np.minimum(rng, nu-rng)
c = np.pi*r2lamb/nu2**2
self._deconv_u = np.roll(np.exp(c*k**2), -nu//2).reshape((-1, 1))
rng = np.arange(nv)
k = np.minimum(rng, nv-rng)
c = np.pi*r2lamb/nv2**2
self._deconv_v = np.roll(
np.exp(c*k**2)/r2lamb, -nv//2).reshape((1, -1))
class _RestOperator(LinearOperator):
def __init__(self, dirty_domain, grid_domain, gconf):
self._domain = makeDomain(grid_domain)
self._target = makeDomain(dirty_domain)
self._gconf = gconf
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
nu, nv = self._target.shape
res = x.to_global_data()
if mode == self.TIMES:
res = hartley(res)
res = np.roll(res, (nu//2, nv//2), axis=(0, 1))
res = res[:nu, :nv]
res *= self._deconv_u
res *= self._deconv_v
res = self._gconf.grid2dirty(res)
else:
res = res*self._deconv_u
res *= self._deconv_v
nu2, nv2 = self._domain.shape
res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant',
constant_values=0)
res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1))
res = hartley(res)
res = self._gconf.dirty2grid(res)
return from_global_data(self._tgt(mode), res)
class RadioGridder(LinearOperator):
def __init__(self, target, nspread, r2lamb, uv):
def __init__(self, grid_domain, bl, gconf, idx):
self._domain = DomainTuple.make(
UnstructuredDomain((uv.shape[0],)))
self._target = DomainTuple.make(target)
UnstructuredDomain((bl.Nrows())))
self._target = DomainTuple.make(grid_domain)
self._bl = bl
self._gconf = gconf
self._idx = idx
self._capability = self.TIMES | self.ADJOINT_TIMES
self._nspread, self._r2lamb = int(nspread), float(r2lamb)
self._uv = uv # FIXME: should we write-protect this?
def apply(self, x, mode):
from nifty_gridder import (to_grid, to_grid_post,
from_grid, from_grid_pre)
import nifty_gridder
self._check_input(x, mode)
nu2, nv2 = self._target.shape
x = x.to_global_data()
if mode == self.TIMES:
res = to_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
res = to_grid_post(res)
x = self._bl.ms2vis(x.to_global_data().reshape((-1, 1)), self._idx)
res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x)
else:
x = from_grid_pre(x)
res = from_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx,
x.to_global_data())
res = self._bl.vis2ms(res, self._idx).reshape((-1,))
return from_global_data(self._tgt(mode), res)
# 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) 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)
......@@ -28,12 +28,15 @@ class Adder(Operator):
field : Field or MultiField
The field by which the input is shifted.
"""
def __init__(self, field):
def __init__(self, field, neg=False):
if not isinstance(field, (Field, MultiField)):
raise TypeError
self._field = field
self._domain = self._target = field.domain
self._neg = bool(neg)
def apply(self, x):
self._check_input(x)
if self._neg:
return x - self._field
return x + self._field
......@@ -105,7 +105,7 @@ class DiagonalOperator(EndomorphicOperator):
res._spaces = None
else:
res._spaces = tuple(set(self._spaces) | set(spc))
res._ldiag = ldiag
res._ldiag = np.array(ldiag)
res._fill_rest()
return res<