diff --git a/Dockerfile b/Dockerfile index b1521076f3a7770f71550d92a6c58602c57c95eb..afe594fb517bb4681951cd84086b999949f16dfc 100644 --- a/Dockerfile +++ b/Dockerfile @@ -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/* diff --git a/README.md b/README.md index 7a7c50d9902db9f8d5f9f76504a971f7bb26ebc7..269dda4cfb7bbb5fd36ff399dcab9865bdc169ba 100644 --- a/README.md +++ b/README.md @@ -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 diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py index 02f5f8f70d1671717ada7d15f7102509851a7109..48ae4f3b1b8a3ff9390581695fcdf31e53c25d81 100644 --- a/demos/bench_gridder.py +++ b/demos/bench_gridder.py @@ -5,12 +5,11 @@ import numpy as np import nifty5 as ift -ift.fft.enable_fftw() np.random.seed(40) N0s, a0s, b0s, c0s = [], [], [], [] -for ii in range(10, 23): +for ii in range(10, 26): nu = 1024 nv = 1024 N = int(2**ii) @@ -28,17 +27,15 @@ 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) diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py index a4f6458ff3e6713d3bc228198b3e43ca7fedc5c5..00875a34c4227304927c2b1b9a9e5b34e0a0c595 100644 --- a/demos/getting_started_3.py +++ b/demos/getting_started_3.py @@ -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.) diff --git a/docs/source/citations.rst b/docs/source/citations.rst index 98115a7ba28d562117c2b34101518536e01f31be..3d7315907b5decaefd99785b327cf8fbb44d5b30 100644 --- a/docs/source/citations.rst +++ b/docs/source/citations.rst @@ -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}}, diff --git a/docs/source/installation.rst b/docs/source/installation.rst index cc5d3d5ba5b0e6ac7cf82489145dbb21964b3494..73bf035232f63845762e74845ab259c4b50fe725 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -9,33 +9,17 @@ 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 -Support for the radio interferometry gridder is added via: +Support for the radio interferometry gridder is added via:: pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git diff --git a/nifty5/domains/log_rg_space.py b/nifty5/domains/log_rg_space.py index 65cfb46a50aa3a9e16b26a32b6d6d5990ae798d1..4cd66c03009bd46784ff8bab3220341b6591e573 100644 --- a/nifty5/domains/log_rg_space.py +++ b/nifty5/domains/log_rg_space.py @@ -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 diff --git a/nifty5/extra.py b/nifty5/extra.py index d48ea3cf9369b306a1f02db6d1c224e2f6a96770..8d8760f26cde3925b03ed55a5253c8b13c70f0e4 100644 --- a/nifty5/extra.py +++ b/nifty5/extra.py @@ -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) diff --git a/nifty5/fft.py b/nifty5/fft.py index 0fc2381409bc92ef66b28f930de95115fcd4885c..714e502b4274a45c73540eedfbbdc193f6d8db35 100644 --- a/nifty5/fft.py +++ b/nifty5/fft.py @@ -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) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index ed51ee6f1e340462a0e8f0365fea79027489c8d0..c9d33ed9283bf18bd7b573b819e91a3594f387e3 100644 --- a/nifty5/library/gridder.py +++ b/nifty5/library/gridder.py @@ -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