Merge remote-tracking branch 'origin/NIFTy_5' into improving_consistency_checks

.gitignore

@@ -2,6 +2,7 @@
 git_version.py
 
 # custom
+*.txt
 setup.cfg
 .idea
 .DS_Store

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-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/*
 

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,35 +73,19 @@ 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:
 
-    pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
+    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
 
 MPI support is added via:
 

demos/bench_gridder.py

@@ -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))

demos/getting_started_3.py

@@ -103,13 +103,15 @@ if __name__ == '__main__':
     data = signal_response(mock_position) + N.draw_sample()
 
     # Minimization parameters
-    ic_sampling = ift.GradientNormController(iteration_limit=100)
-    ic_newton = ift.GradInfNormController(
-        name='Newton', tol=1e-7, iteration_limit=35)
+    ic_sampling = ift.AbsDeltaEnergyController(
+        name='Sampling', deltaE=0.05, iteration_limit=100)
+    ic_newton = ift.AbsDeltaEnergyController(
+        name='Newton', deltaE=0.5, iteration_limit=35)
     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.)

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}},
@@ -11,7 +17,7 @@ NIFTy-related publications
       date = {2018-04-05},
     }
 
-    @ARTICLE{2013A&A...554A..26S,
+    @article{2013A&A...554A..26S,
       author = {{Selig}, M. and {Bell}, M.~R. and {Junklewitz}, H. and {Oppermann}, N. and {Reinecke}, M. and {Greiner}, M. and {Pachajoa}, C. and {En{\ss}lin}, T.~A.},
       title = "{NIFTY - Numerical Information Field Theory. A versatile PYTHON library for signal inference}",
       journal = {\aap},
@@ -29,7 +35,7 @@ NIFTy-related publications
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
     }
 
-    @ARTICLE{2017arXiv170801073S,
+    @article{2017arXiv170801073S,
       author = {{Steininger}, T. and {Dixit}, J. and {Frank}, P. and {Greiner}, M. and {Hutschenreuter}, S. and {Knollm{\"u}ller}, J. and {Leike}, R. and {Porqueres}, N. and {Pumpe}, D. and {Reinecke}, M. and {{\v S}raml}, M. and {Varady}, C. and {En{\ss}lin}, T.},
       title = "{NIFTy 3 - Numerical Information Field Theory - A Python framework for multicomponent signal inference on HPC clusters}",
       journal = {ArXiv e-prints},

docs/source/installation.rst

@@ -9,35 +9,19 @@ 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
+    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
 
 MPI support is added via::
 

nifty5/__init__.py

@@ -45,20 +45,22 @@ 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, PartialExtractor)
 from .operators.value_inserter import ValueInserter
 from .operators.energy_operators import (
     EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
-    BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
+    BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator,
+    Squared2NormOperator)
 from .operators.convolution_operators import FuncConvolutionOperator
 
 from .probing import probe_with_posterior_samples, probe_diagonal, \
-    StatCalculator
+    StatCalculator, approximation2endo
 
 from .minimization.line_search import LineSearch
 from .minimization.iteration_controllers import (
     IterationController, GradientNormController, DeltaEnergyController,
-    GradInfNormController)
+    GradInfNormController, AbsDeltaEnergyController)
 from .minimization.minimizer import Minimizer
 from .minimization.conjugate_gradient import ConjugateGradient
 from .minimization.nonlinear_cg import NonlinearCG
@@ -86,7 +88,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
@@ -97,6 +98,8 @@ from .logger import logger
 
 from .linearization import Linearization
 
+from .operator_spectrum import operator_spectrum
+
 from . import internal_config
 _scheme = internal_config.parallelization_scheme()
 if _scheme == "Samples":

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

nifty5/domains/log_rg_space.py

@@ -50,6 +50,8 @@ class LogRGSpace(StructuredDomain):
 
         self._bindistances = tuple(bindistances)
         self._t_0 = tuple(t_0)
+        if min(self._bindistances) <= 0:
+            raise ValueError('Non-positive bindistances encountered')
 
         self._dim = int(reduce(lambda x, y: x*y, self._shape))
         self._dvol = float(reduce(lambda x, y: x*y, self._bindistances))
@@ -80,8 +82,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
@@ -91,10 +93,10 @@ class LogRGSpace(StructuredDomain):
         Returns
         -------
         LogRGSpace
-            The parter domain
+            The partner domain
         """
         codomain_bindistances = 1./(self.bindistances*self.shape)
-        return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
+        return LogRGSpace(self.shape, codomain_bindistances, self._t_0, not self.harmonic)
 
     def get_k_length_array(self):
         """Generates array of distances to origin of the space.

nifty5/domains/power_space.py

@@ -165,6 +165,8 @@ class PowerSpace(StructuredDomain):
 
         if binbounds is not None:
             binbounds = tuple(binbounds)
+            if min(binbounds) < 0:
+                raise ValueError('Negative binbounds encountered')
 
         key = (harmonic_partner, binbounds)
         if self._powerIndexCache.get(key) is None:

nifty5/domains/rg_space.py

@@ -54,6 +54,8 @@ class RGSpace(StructuredDomain):
         if np.isscalar(shape):
             shape = (shape,)
         self._shape = tuple(int(i) for i in shape)
+        if min(self._shape) < 0:
+            raise ValueError('Negative number of pixels encountered')
 
         if distances is None:
             if self.harmonic:
@@ -66,6 +68,8 @@ class RGSpace(StructuredDomain):
             temp = np.empty(len(self.shape), dtype=np.float64)
             temp[:] = distances
             self._distances = tuple(temp)
+        if min(self._distances) <= 0:
+            raise ValueError('Non-positive distances encountered')
 
         self._dvol = float(reduce(lambda x, y: x*y, self._distances))
         self._size = int(reduce(lambda x, y: x*y, self._shape))
@@ -177,7 +181,7 @@ class RGSpace(StructuredDomain):
         Returns
         -------
         RGSpace
-            The parter domain
+            The partner domain
         """
         distances = 1. / (np.array(self.shape)*np.array(self.distances))
         return RGSpace(self.shape, distances, not self.harmonic)

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
 
@@ -68,14 +70,25 @@ def _full_implementation(op, domain_dtype, target_dtype, atol, rtol,
 
 
 def _check_linearity(op, domain_dtype, atol, rtol):
+    needed_cap = op.TIMES
+    if (op.capability & needed_cap) != needed_cap:
+        return
     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,7 +122,11 @@ 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)
+    _check_linearity(op.adjoint, target_dtype, atol, rtol)
+    _check_linearity(op.inverse, target_dtype, atol, rtol)
+    _check_linearity(op.adjoint.inverse, domain_dtype, atol, rtol)
     _full_implementation(op, domain_dtype, target_dtype, atol, rtol,
                          only_r_linear)
     _full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol,
@@ -180,6 +197,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)

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 = int(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)

nifty5/library/gridder.py

@@ -20,109 +20,92 @@ 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
 
 
 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