Commit a439a47c authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'new_gridder' into 'NIFTy_5'

New gridder (again!)

See merge request !324
parents afbc22bc b667dc47
Pipeline #50416 passed with stages
in 19 minutes
......@@ -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-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,7 +5,6 @@ import numpy as np
import nifty5 as ift
ift.fft.enable_fftw()
np.random.seed(40)
N0s, a0s, b0s, c0s = [], [], [], []
......
......@@ -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
......
......@@ -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)
......
......@@ -17,107 +17,52 @@
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 set_nthreads(nthr):
global _nthreads
_nthreads = nthr
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
# FIXME this should not be necessary ... no one should call a complex FFT
# with a float array.
def _make_complex(a):
if a.dtype in (np.complex64, np.complex128):
return a
if a.dtype == np.float64:
return a.astype(np.complex128)
if a.dtype == np.float32:
return a.astype(np.complex64)
raise NotImplementedError
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)
return pypocketfft.fftn(_make_complex(a), axes=axes, nthreads=_nthreads)
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.rfftn(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)
# FIXME this is a temporary fix and can be done more elegantly
if axes is None:
fct = 1./a.size
else:
return np.fft.ifftn(a, axes=axes)
fct = 1./np.prod(np.take(a.shape, axes))
return pypocketfft.ifftn(_make_complex(a), axes=axes, fct=fct,
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)
return pypocketfft.hartley2(a, axes=axes, nthreads=_nthreads)
# Do a real-to-complex forward FFT and return the _full_ output array
......
......@@ -26,7 +26,8 @@ from ..sugar import from_global_data, makeDomain
class GridderMaker(object):
def __init__(self, domain, eps=1e-15):
def __init__(self, domain, eps=2e-13):
from nifty_gridder import get_w
domain = makeDomain(domain)
if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or
not len(domain.shape) == 2):
......@@ -34,20 +35,17 @@ class GridderMaker(object):
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))
nu2, nv2 = 2*nu, 2*nv
w = get_w(eps)
nsafe = (w+1)//2
nu2 = max([nu2, 2*nsafe])
nv2 = max([nv2, 2*nsafe])
oversampled_domain = RGSpace(
[nu2, nv2], distances=[1, 1], harmonic=False)
self._nspread = nspread
self._r2lamb = r2lamb
self._rest = _RestOperator(domain, oversampled_domain, r2lamb)
self._eps = eps
self._rest = _RestOperator(domain, oversampled_domain, eps)
def getReordering(self, uv):
from nifty_gridder import peanoindex
......@@ -55,7 +53,7 @@ class GridderMaker(object):
return peanoindex(uv, nu2, nv2)
def getGridder(self, uv):
return RadioGridder(self._rest.domain, self._nspread, self._r2lamb, uv)
return RadioGridder(self._rest.domain, self._eps, uv)
def getRest(self):
return self._rest
......@@ -65,22 +63,22 @@ class GridderMaker(object):
class _RestOperator(LinearOperator):
def __init__(self, domain, oversampled_domain, r2lamb):
def __init__(self, domain, oversampled_domain, eps):
from nifty_gridder import correction_factors
self._domain = makeDomain(oversampled_domain)
self._target = domain
nu, nv = domain.shape
nu2, nv2 = oversampled_domain.shape
fu = correction_factors(nu2, nu//2+1, eps)
fv = correction_factors(nv2, nv//2+1, eps)
# 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))
self._deconv_u = np.roll(fu[k], -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))
self._deconv_v = np.roll(fv[k], -nv//2).reshape((1, -1))
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
......@@ -105,24 +103,20 @@ class _RestOperator(LinearOperator):
class RadioGridder(LinearOperator):
def __init__(self, target, nspread, r2lamb, uv):
def __init__(self, target, eps, uv):
self._domain = DomainTuple.make(
UnstructuredDomain((uv.shape[0],)))
self._target = DomainTuple.make(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._nspread, self._r2lamb = int(nspread), float(r2lamb)
self._eps = float(eps)
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)
from nifty_gridder import to_grid, from_grid
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)
nu2, nv2 = self._target.shape
res = to_grid(self._uv, x.to_global_data(), nu2, nv2, self._eps)
else:
x = from_grid_pre(x)
res = from_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
res = from_grid(self._uv, x.to_global_data(), self._eps)
return from_global_data(self._tgt(mode), res)
......@@ -349,11 +349,10 @@ def _plot2D(f, ax, **kwargs):
rgb = _rgb_data(f.to_global_data())
have_rgb = True
label = kwargs.pop("label", None)
foo = kwargs.pop("norm", None)
norm = {} if foo is None else {'norm': foo}
aspect = kwargs.pop("aspect", None)
foo = kwargs.pop("aspect", None)
aspect = {} if foo is None else {'aspect': foo}
ax.set_title(kwargs.pop("title", ""))
......
......@@ -37,6 +37,8 @@ _pow_spaces = [ift.PowerSpace(ift.RGSpace((17, 38), harmonic=True))]
pmp = pytest.mark.parametrize
dtype = list2fixture([np.float64, np.complex128])
np.random.seed(42)
@pmp('sp', _p_RG_spaces)
def testLOSResponse(sp, dtype):
......
......@@ -34,22 +34,10 @@ def _get_rtol(tp):
pmp = pytest.mark.parametrize
dtype = list2fixture([np.float64, np.float32, np.complex64, np.complex128])
op = list2fixture([ift.HartleyOperator, ift.FFTOperator])
fftw = list2fixture([False, True])
def test_switch():
ift.fft.enable_fftw()
assert_(ift.fft._use_fftw is True)
ift.fft.disable_fftw()
assert_(ift.fft._use_fftw is False)
ift.fft.enable_fftw()
assert_(ift.fft._use_fftw is True)
@pmp('d', [0.1, 1, 3.7])
def test_fft1D(d, dtype, op, fftw):
if fftw:
ift.fft.enable_fftw()
def test_fft1D(d, dtype, op):
dim1 = 16
tol = _get_rtol(dtype)
a = ift.RGSpace(dim1, distances=d)
......@@ -69,16 +57,16 @@ def test_fft1D(d, dtype, op, fftw):
domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
ift.fft.disable_fftw()
@pmp('dim1', [12, 15])
@pmp('dim2', [9, 12])
@pmp('d1', [0.1, 1, 3.7])
@pmp('d2', [0.4, 1, 2.7])
def test_fft2D(dim1, dim2, d1, d2, dtype, op, fftw):
if fftw:
ift.fft.enable_fftw()
@pmp('nthreads', [0, 1, 2, 3, 4])
def test_fft2D(dim1, dim2, d1, d2, dtype, op, nthreads):
ift.fft.set_nthreads(nthreads)
assert_(ift.fft.nthreads() == nthreads)
tol = _get_rtol(dtype)
a = ift.RGSpace([dim1, dim2], distances=[d1, d2])
b = ift.RGSpace(
......@@ -97,13 +85,11 @@ def test_fft2D(dim1, dim2, d1, d2, dtype, op, fftw):
domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
ift.fft.disable_fftw()
ift.fft.set_nthreads(1)
@pmp('index', [0, 1, 2])
def test_composed_fft(index, dtype, op, fftw):
if fftw:
ift.fft.enable_fftw()
def test_composed_fft(index, dtype, op):
tol = _get_rtol(dtype)
a = [a1, a2,
a3] = [ift.RGSpace((32,)),
......@@ -115,7 +101,6 @@ def test_composed_fft(index, dtype, op, fftw):
domain=(a1, a2, a3), random_type='normal', std=7, mean=3, dtype=dtype)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
ift.fft.disable_fftw()
@pmp('space', [
......@@ -123,9 +108,7 @@ def test_composed_fft(index, dtype, op, fftw):
ift.RGSpace((15, 27), distances=(.7, .33), harmonic=True),
ift.RGSpace(73, distances=0.5643)
])
def test_normalisation(space, dtype, op, fftw):
if fftw:
ift.fft.enable_fftw()
def test_normalisation(space, dtype, op):
tol = 10*_get_rtol(dtype)
cospace = space.get_default_codomain()
fft = op(space, cospace)
......@@ -138,4 +121,3 @@ def test_normalisation(space, dtype, op, fftw):
assert_allclose(
inp.to_global_data()[zero_idx], out.integrate(), rtol=tol, atol=tol)
assert_allclose(out.local_data, out2.local_data, rtol=tol, atol=tol)
ift.fft.disable_fftw()
......@@ -17,7 +17,7 @@
import numpy as np
import pytest
from numpy.testing import assert_allclose
from numpy.testing import assert_allclose, assert_
import nifty5 as ift
......@@ -26,15 +26,20 @@ np.random.seed(40)
pmp = pytest.mark.parametrize
def _l2error(a, b):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))
@pmp('eps', [1e-2, 1e-4, 1e-7, 1e-10, 1e-11, 1e-12, 2e-13])
@pmp('nu', [12, 128])
@pmp('nv', [4, 12, 128])
@pmp('N', [1, 10, 100])
def test_gridding(nu, nv, N):
def test_gridding(nu, nv, N, eps):
uv = np.random.rand(N, 2) - 0.5
vis = np.random.randn(N) + 1j*np.random.randn(N)
# Nifty
GM = ift.GridderMaker(ift.RGSpace((nu, nv)))
GM = ift.GridderMaker(ift.RGSpace((nu, nv)), eps=eps)
# re-order for performance
idx = GM.getReordering(uv)
uv, vis = uv[idx], vis[idx]
......@@ -48,17 +53,17 @@ def test_gridding(nu, nv, N):
dft = pynu*0.
for i in range(N):
dft += (vis[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))).real
assert_allclose(dft, pynu)
assert_(_l2error(dft, pynu) < eps)
@pmp('eps', [1e-2, 1e-6, 1e-15])
@pmp('eps', [1e-2, 1e-6, 2e-13])
@pmp('nu', [12, 128])
@pmp('nv', [4, 12, 128])
@pmp('N', [1, 10, 100])
def test_build(nu, nv, N, eps):
dom = ift.RGSpace([nu, nv])
uv = np.random.rand(N, 2) - 0.5
GM = ift.GridderMaker(dom)
GM = ift.GridderMaker(dom, eps=eps)
# re-order for performance
idx = GM.getReordering(uv)
uv = uv[idx]
......
......@@ -23,23 +23,23 @@ import nifty5 as ift
def test_simplification():
from nifty5.operators.operator import _ConstantOperator
f1 = ift.Field.full(ift.RGSpace(10),2.)
f1 = ift.Field.full(ift.RGSpace(10), 2.)
op = ift.FFTOperator(f1.domain)
_, op2 = op.simplify_for_constant_input(f1)
assert_equal(isinstance(op2, _ConstantOperator), True)
assert_allclose(op(f1).local_data, op2(f1).local_data)
dom = {"a": ift.RGSpace(10)}
f1 = ift.full(dom,2.)
f1 = ift.full(dom, 2.)
op = ift.FFTOperator(f1.domain["a"]).ducktape("a")
_, op2 = op.simplify_for_constant_input(f1)
assert_equal(isinstance(op2, _ConstantOperator), True)
assert_allclose(op(f1).local_data, op2(f1).local_data)
dom = {"a": ift.RGSpace(10), "b": ift.RGSpace(5)}
f1 = ift.full(dom,2.)
f1 = ift.full(dom, 2.)
pdom = {"a": ift.RGSpace(10)}
f2 = ift.full(pdom,2.)
f2 = ift.full(pdom, 2.)
o1 = ift.FFTOperator(f1.domain["a"])
o2 = ift.FFTOperator(f1.domain["b"])
op = (o1.ducktape("a").ducktape_left("a") +
......
Markdown is supported
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