Merge branch 'NIFTy_5' into documentation_feedback

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 \
+    libfftw3-dev python3-mpi4py python3-matplotlib python3-pynfft \
   # more optional NIFTy dependencies
   && pip3 install pyfftw \
   && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \

nifty5/__init__.py

@@ -74,7 +74,8 @@ from .minimization.metric_gaussian_kl import MetricGaussianKL
 from .sugar import *
 from .plot import Plot
 
-from .library.smooth_linear_amplitude import SLAmplitude, CepstrumOperator
+from .library.smooth_linear_amplitude import (
+    SLAmplitude, LinearSLAmplitude, CepstrumOperator)
 from .library.inverse_gamma_operator import InverseGammaOperator
 from .library.los_response import LOSResponse
 from .library.dynamic_operator import (dynamic_operator,
@@ -85,6 +86,7 @@ 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 . import extra
 

nifty5/data_objects/numpy_do.py

@@ -149,7 +149,7 @@ def ensure_default_distributed(arr):
 
 
 def absmax(arr):
-    return np.linalg.norm(arr.rehape(-1), ord=np.inf)
+    return np.linalg.norm(arr.reshape(-1), ord=np.inf)
 
 
 def norm(arr, ord=2):

nifty5/library/nfft.py

+# 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)

nifty5/library/smooth_linear_amplitude.py

@@ -169,6 +169,18 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
         which returns on its target a power spectrum which consists out of a
         smooth and a linear part.
     '''
+    return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0, sm=sm,
+                             sv=sv, im=im, iv=iv, keys=keys).exp()
+
+
+def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
+                      keys=['tau', 'phi']):
+    '''LinearOperator for parametrizing smooth log-amplitudes (square roots of
+    power spectra).
+
+    Logarithm of SLAmplitude
+    See documentation of SLAmplitude for more details
+    '''
     if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)):
         raise TypeError
 
@@ -196,4 +208,4 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
     loglog_ampl = 0.5*(smooth + linear)
 
     # Go from loglog-space to linear-linear-space
-    return et @ loglog_ampl.exp()
+    return et @ loglog_ampl

nifty5/operators/energy_operators.py

@@ -20,6 +20,7 @@ import numpy as np
 from .. import utilities
 from ..domain_tuple import DomainTuple
 from ..field import Field
+from ..multi_field import MultiField
 from ..linearization import Linearization
 from ..sugar import makeDomain, makeOp
 from .linear_operator import LinearOperator
@@ -121,7 +122,7 @@ class GaussianEnergy(EnergyOperator):
     """
 
     def __init__(self, mean=None, covariance=None, domain=None):
-        if mean is not None and not isinstance(mean, Field):
+        if mean is not None and not isinstance(mean, (Field, MultiField)):
             raise TypeError
         if covariance is not None and not isinstance(covariance,
                                                      LinearOperator):
@@ -307,7 +308,6 @@ class StandardHamiltonian(EnergyOperator):
         Tells an internal :class:`SamplingEnabler` which convergence criterion
         to use to draw Gaussian samples.
 
-
     See also
     --------
     `Encoding prior knowledge in the structure of the likelihood`,

nifty5/operators/operator.py

@@ -384,7 +384,6 @@ class _OpSum(Operator):
         v = x._val if lin else x
         v1 = v.extract(self._op1.domain)
         v2 = v.extract(self._op2.domain)
-        res = None
         if not lin:
             return self._op1(v1).unite(self._op2(v2))
         wm = x.want_metric
@@ -393,7 +392,10 @@ class _OpSum(Operator):
         op = lin1._jac._myadd(lin2._jac, False)
         res = lin1.new(lin1._val.unite(lin2._val), op(x.jac))
         if lin1._metric is not None and lin2._metric is not None:
-            res = res.add_metric(lin1._metric + lin2._metric)
+            from .sandwich_operator import SandwichOperator
+            met = lin1._metric._myadd(lin2._metric, False)
+            met = SandwichOperator.make(x.jac, met)
+            res = res.add_metric(met)
         return res
 
     def _simplify_for_constant_input_nontrivial(self, c_inp):

nifty5/operators/scaling_operator.py

@@ -66,9 +66,12 @@ class ScalingOperator(EndomorphicOperator):
             return x
         if fct == 0.:
             return full(self.domain, 0.)
-        if (mode & 10) != 0:
+
+        MODES_WITH_ADJOINT = self.ADJOINT_TIMES | self.ADJOINT_INVERSE_TIMES
+        MODES_WITH_INVERSE = self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES
+        if (mode & MODES_WITH_ADJOINT) != 0:
             fct = np.conj(fct)
-        if (mode & 12) != 0:
+        if (mode & MODES_WITH_INVERSE) != 0:
             fct = 1./fct
         return x*fct
 

nifty5/sugar.py

@@ -29,6 +29,7 @@ from .multi_field import MultiField
 from .operators.block_diagonal_operator import BlockDiagonalOperator
 from .operators.diagonal_operator import DiagonalOperator
 from .operators.distributors import PowerDistributor
+from .plot import Plot
 
 __all__ = ['PS_field', 'power_analyze', 'create_power_operator',
            'create_harmonic_smoothing_operator', 'from_random',
@@ -37,7 +38,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator',
            'sin', 'cos', 'tan', 'sinh', 'cosh',
            'absolute', 'one_over', 'clip', 'sinc',
            'conjugate', 'get_signal_variance', 'makeOp', 'domain_union',
-           'get_default_codomain']
+           'get_default_codomain', 'single_plot']
 
 
 def PS_field(pspace, func):
@@ -434,3 +435,14 @@ def get_default_codomain(domainoid, space=None):
     ret = [dom for dom in domainoid]
     ret[space] = domainoid[space].get_default_codomain()
     return DomainTuple.make(ret)
+
+
+def single_plot(field, **kwargs):
+    """Creates a single plot using `Plot`.
+    Keyword arguments are passed to both `Plot.add` and `Plot.output`.
+    """
+    p = Plot()
+    p.add(field, **kwargs)
+    if 'title' in kwargs:
+        del(kwargs['title'])
+    p.output(**kwargs)

test/test_operators/test_adjoint.py

@@ -279,3 +279,10 @@ 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)