Merge remote-tracking branch 'origin/HEAD' into operator_spectra

README.md

@@ -85,7 +85,7 @@ Support for spherical harmonic transforms 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:
 

docs/source/installation.rst

@@ -21,7 +21,7 @@ Support for spherical harmonic transforms 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/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))

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

nifty5/linearization.py

@@ -20,6 +20,7 @@ import numpy as np
 from .field import Field
 from .multi_field import MultiField
 from .sugar import makeOp
+from . import utilities
 
 
 class Linearization(object):
@@ -108,7 +109,6 @@ class Linearization(object):
         return self._metric
 
     def __getitem__(self, name):
-        from .operators.simple_linear_operators import ducktape
         return self.new(self._val[name], self._jac.ducktape_left(name))
 
     def __neg__(self):
@@ -298,6 +298,10 @@ class Linearization(object):
             tmp2 = makeOp(1. - (tmp == min) - (tmp == max))
         return self.new(tmp, tmp2(self._jac))
 
+    def sqrt(self):
+        tmp = self._val.sqrt()
+        return self.new(tmp, makeOp(0.5/tmp)(self._jac))
+
     def sin(self):
         tmp = self._val.sin()
         tmp2 = self._val.cos()
@@ -315,7 +319,11 @@ class Linearization(object):
 
     def sinc(self):
         tmp = self._val.sinc()
-        tmp2 = (self._val.cos()-tmp)/self._val
+        tmp2 = ((np.pi*self._val).cos()-tmp)/self._val
+        ind = self._val.local_data == 0
+        loc = tmp2.local_data.copy()
+        loc[ind] = 0
+        tmp2 = Field.from_local_data(tmp.domain, loc)
         return self.new(tmp, makeOp(tmp2)(self._jac))
 
     def log(self):
@@ -342,8 +350,16 @@ class Linearization(object):
         return self.new(tmp2, makeOp(0.5*(1.-tmp**2))(self._jac))
 
     def absolute(self):
+        if utilities.iscomplextype(self._val.dtype):
+            raise TypeError("Argument must not be complex")
         tmp = self._val.absolute()
         tmp2 = self._val.sign()
+
+        ind = self._val.local_data == 0
+        loc = tmp2.local_data.copy().astype(float)
+        loc[ind] = np.nan
+        tmp2 = Field.from_local_data(tmp.domain, loc)
+
         return self.new(tmp, makeOp(tmp2)(self._jac))
 
     def one_over(self):

nifty5/operators/convolution_operators.py

@@ -15,19 +15,16 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import numpy as np
 
-from ..domains.rg_space import RGSpace
-from ..domains.lm_space import LMSpace
-from ..domains.hp_space import HPSpace
+from .. import utilities
+from ..domain_tuple import DomainTuple
 from ..domains.gl_space import GLSpace
+from ..domains.hp_space import HPSpace
+from ..domains.rg_space import RGSpace
+from .diagonal_operator import DiagonalOperator
 from .endomorphic_operator import EndomorphicOperator
 from .harmonic_operators import HarmonicTransformOperator
-from .diagonal_operator import DiagonalOperator
 from .simple_linear_operators import WeightApplier
-from ..domain_tuple import DomainTuple
-from ..field import Field
-from .. import utilities
 
 
 def FuncConvolutionOperator(domain, func, space=None):
@@ -77,4 +74,25 @@ def _ConvolutionOperator(domain, kernel, space=None):
     HT = HarmonicTransformOperator(lm, domain[space], space)
     diag = DiagonalOperator(kernel*domain[space].total_volume, lm, (space,))
     wgt = WeightApplier(domain, space, 1)
-    return HT(diag(HT.adjoint(wgt)))
+    op = HT(diag(HT.adjoint(wgt)))
+    return _ApplicationWithoutMeanOperator(op)
+
+
+class _ApplicationWithoutMeanOperator(EndomorphicOperator):
+    def __init__(self, op):
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+        if op.domain != op.target:
+            raise TypeError("Operator needs to be endomorphic")
+        self._domain = op.domain
+        self._op = op
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+        mean = x.mean()
+        return mean + self._op.apply(x - mean, mode)
+
+    def __repr__(self):
+        from ..utilities import indent
+        return "\n".join((
+            "_ApplicationWithoutMeanOperator:",
+            indent(self._op.__repr__())))

nifty5/operators/distributors.py

@@ -30,7 +30,7 @@ class DOFDistributor(LinearOperator):
     """Operator which distributes actual degrees of freedom (dof) according to
     some distribution scheme into a higher dimensional space. This distribution
     scheme is defined by the dofdex, a degree of freedom index, which
-    associates the entries within the operators domain to locations in its
+    associates the entries within the operator's domain to locations in its
     target. This operator's domain is a DOFSpace, which is defined according to
     the distribution scheme.
 

nifty5/operators/sum_operator.py

@@ -23,7 +23,6 @@ from ..sugar import domain_union
 from ..utilities import indent
 from .block_diagonal_operator import BlockDiagonalOperator
 from .linear_operator import LinearOperator
-from .simple_linear_operators import NullOperator
 
 
 class SumOperator(LinearOperator):

nifty5/plot.py

@@ -151,17 +151,17 @@ def _rgb_data(spectral_cube):
         res /= tmp
         return res
 
+    shp = spectral_cube.shape[:-1]+(3,)
     spectral_cube = spectral_cube.reshape((-1, spectral_cube.shape[-1]))
     xyz = getxyz(spectral_cube.shape[-1])
     xyz_data = np.tensordot(spectral_cube, xyz, axes=[-1, -1])
     xyz_data /= xyz_data.max()
     xyz_data = to_logscale(xyz_data, max(1e-3, xyz_data.min()), 1.)
     rgb_data = xyz_data.copy()
-    it = np.nditer(xyz_data[:, 0], flags=['multi_index'])
     for x in range(xyz_data.shape[0]):
         rgb_data[x] = _gammacorr(np.matmul(MATRIX_SRGB_D65, xyz_data[x]))
     rgb_data = rgb_data.clip(0., 1.)
-    return rgb_data.reshape(spectral_cube.shape[:-1]+(-1,))
+    return rgb_data.reshape(shp)
 
 
 def _find_closest(A, target):

nifty5/utilities.py

@@ -240,6 +240,8 @@ def special_add_at(a, axis, index, b):
 
 
 _iscomplex_tpl = (np.complex64, np.complex128)
+
+
 def iscomplextype(dtype):
     return dtype.type in _iscomplex_tpl
 

test/test_field.py

@@ -57,8 +57,6 @@ def test_power_synthesize_analyze(space1, space2):
     fp1 = ift.PS_field(p1, _spec1)
     p2 = ift.PowerSpace(space2)
     fp2 = ift.PS_field(p2, _spec2)
-    outer = np.outer(fp1.to_global_data(), fp2.to_global_data())
-    fp = ift.Field.from_global_data((p1, p2), outer)
 
     op1 = ift.create_power_operator((space1, space2), _spec1, 0)
     op2 = ift.create_power_operator((space1, space2), _spec2, 1)
@@ -345,11 +343,11 @@ def test_funcs(num, dom, func):
 @pmp('dtype', [np.float64, np.complex128])
 def test_from_random(rtype, dtype):
     sp = ift.RGSpace(3)
-    f = ift.Field.from_random(rtype, sp, dtype=dtype)
+    ift.Field.from_random(rtype, sp, dtype=dtype)
 
 
 def test_field_of_objects():
     arr = np.array(['x', 'y', 'z'])
     sp = ift.RGSpace(3)
     with assert_raises(TypeError):
-        f = ift.Field.from_global_data(sp, arr)
+        ift.Field.from_global_data(sp, arr)

test/test_linearization.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) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+import pytest
+from numpy.testing import assert_, assert_allclose
+
+import nifty5 as ift
+
+pmp = pytest.mark.parametrize
+
+
+def _lin2grad(lin):
+    return lin.jac(ift.full(lin.domain, 1.)).to_global_data()
+
+
+def jt(lin, check):
+    assert_allclose(_lin2grad(lin), check)
+
+
+def test_special_gradients():
+    dom = ift.UnstructuredDomain((1,))
+    f = ift.full(dom, 2.4)
+    var = ift.Linearization.make_var(f)
+    s = f.to_global_data()
+
+    jt(var.clip(0, 10), np.ones_like(s))
+    jt(var.clip(-1, 0), np.zeros_like(s))
+
+    assert_allclose(
+        _lin2grad(ift.Linearization.make_var(0*f).sinc()), np.zeros(s.shape))
+    assert_(np.isnan(_lin2grad(ift.Linearization.make_var(0*f).absolute())))
+    assert_allclose(
+        _lin2grad(ift.Linearization.make_var(0*f + 10).absolute()),
+        np.ones(s.shape))
+    assert_allclose(
+        _lin2grad(ift.Linearization.make_var(0*f - 10).absolute()),
+        -np.ones(s.shape))
+
+
+@pmp('f', [
+    'log', 'exp', 'sqrt', 'sin', 'cos', 'tan', 'sinc', 'sinh', 'cosh', 'tanh',
+    'absolute', 'one_over', 'sigmoid'
+])
+def test_actual_gradients(f):
+    dom = ift.UnstructuredDomain((1,))
+    fld = ift.full(dom, 2.4)
+    eps = 1e-8
+    var0 = ift.Linearization.make_var(fld)
+    var1 = ift.Linearization.make_var(fld + eps)
+    f0 = getattr(var0, f)().val.to_global_data()
+    f1 = getattr(var1, f)().val.to_global_data()
+    df0 = (f1 - f0)/eps
+    df1 = _lin2grad(getattr(var0, f)())
+    assert_allclose(df0, df1, rtol=100*eps)

test/test_operators/test_convolution_operators.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) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from numpy.testing import assert_allclose
+
+import nifty5 as ift
+import numpy as np
+
+from ..common import list2fixture
+
+space = list2fixture([
+    ift.RGSpace(4),
+    ift.HPSpace(4),
+    ift.GLSpace(4)
+])
+
+
+def test_const_func(space):
+    sig = ift.Field.from_random('normal', domain=space)
+    fco_op = ift.FuncConvolutionOperator(space, lambda x: np.ones(x.shape))
+    vals = fco_op(sig).to_global_data()
+    vals = np.round(vals, decimals=5)
+    assert len(np.unique(vals)) == 1
+
+
+def gauss(x, sigma):
+    normalization = np.sqrt(2. * np.pi) * sigma
+    return np.exp(-0.5 * x * x / sigma**2) / normalization
+
+
+def test_gaussian_smoothing():
+    N = 128
+    sigma = N / 10**4
+    dom = ift.RGSpace(N)
+    sig = ift.exp(ift.Field.from_random('normal', dom))
+    fco_op = ift.FuncConvolutionOperator(dom, lambda x: gauss(x, sigma))
+    sm_op = ift.HarmonicSmoothingOperator(dom, sigma)
+    assert_allclose(fco_op(sig).to_global_data(),
+                    sm_op(sig).to_global_data(),
+                    rtol=1e-05)
+    assert_allclose(fco_op.adjoint_times(sig).to_global_data(),
+                    sm_op.adjoint_times(sig).to_global_data(),
+                    rtol=1e-05)

test/test_operators/test_representation.py

@@ -211,7 +211,6 @@ def testExpTransform(args, dtype):
       ((ift.LogRGSpace(10, [2.], [1.]), ift.UnstructuredDomain(13)), 0),
       ((ift.UnstructuredDomain(13), ift.LogRGSpace(17, [3.], [.7])), 1)])
 def testQHTOperator(args):
-    dtype = np.float64
     tgt = ift.DomainTuple.make(args[0])
     _check_repr(ift.QHTOperator(tgt, args[1]))
 
@@ -225,14 +224,11 @@ def testRegridding(args):
     _check_repr(ift.RegriddingOperator(*args))
 
 
-@pmp(
-    'fdomain',
-    [
-        ift.DomainTuple.make((ift.RGSpace(
-            (3, 5, 4)), ift.RGSpace((16,), distances=(7.,))),),
-        ift.DomainTuple.make(ift.HPSpace(12),)
-    ],
-)
+@pmp('fdomain', [
+    ift.DomainTuple.make((ift.RGSpace(
+        (3, 5, 4)), ift.RGSpace((16,), distances=(7.,))),),
+    ift.DomainTuple.make(ift.HPSpace(12),)
+])
 @pmp('domain', [
     ift.DomainTuple.make((ift.RGSpace((2,)), ift.GLSpace(10)),),
     ift.DomainTuple.make(ift.RGSpace((10, 12), distances=(0.1, 1.)),)
@@ -251,5 +247,5 @@ def testValueInserter(sp, seed):
         if ss == 1:
             ind.append(0)
         else:
-            ind.append(np.random.randint(0, ss-1))
+            ind.append(np.random.randint(0, ss - 1))
     _check_repr(ift.ValueInserter(sp, ind))

test/test_operators/test_simplification.py

@@ -15,7 +15,6 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import pytest
 from numpy.testing import assert_allclose, assert_equal
 
 import nifty5 as ift

test/test_spaces/test_lm_space.py

@@ -66,8 +66,7 @@ def get_k_length_array_configs():
 
 @pmp('attribute', ['lmax', 'mmax', 'size'])
 def test_property_ret_type(attribute):
-    l = ift.LMSpace(7, 5)
-    assert_(isinstance(getattr(l, attribute), int))
+    assert_(isinstance(getattr(ift.LMSpace(7, 5), attribute), int))
 
 
 @pmp('lmax, mmax, expected', CONSTRUCTOR_CONFIGS)
@@ -76,9 +75,8 @@ def test_constructor(lmax, mmax, expected):
         with assert_raises(expected['error']):
             ift.LMSpace(lmax, mmax)
     else:
-        l = ift.LMSpace(lmax, mmax)
         for key, value in expected.items():
-            assert_equal(getattr(l, key), value)
+            assert_equal(getattr(ift.LMSpace(lmax, mmax), key), value)
 
 
 def test_dvol():
@@ -87,5 +85,5 @@ def test_dvol():
 
 @pmp('lmax, expected', get_k_length_array_configs())
 def test_k_length_array(lmax, expected):
-    l = ift.LMSpace(lmax)
-    assert_allclose(l.get_k_length_array().to_global_data(), expected)
+    assert_allclose(ift.LMSpace(lmax).get_k_length_array().to_global_data(),
+                    expected)

test/test_spaces/test_power_space.py

@@ -17,11 +17,12 @@
 
 from itertools import chain, product
 
-import pytest
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
 
+import nifty5 as ift
+
 pmp = pytest.mark.parametrize
 
 HARMONIC_SPACES = [

test/test_sugar.py

@@ -16,8 +16,7 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
-import pytest
-from numpy.testing import assert_allclose, assert_equal, assert_raises
+from numpy.testing import assert_equal
 
 import nifty5 as ift
 
@@ -25,8 +24,8 @@ import nifty5 as ift
 def test_get_signal_variance():
     space = ift.RGSpace(3)
     hspace = space.get_default_codomain()
-    spec1 = lambda x: np.ones_like(x)
-    assert_equal(ift.get_signal_variance(spec1, hspace), 3.)
+    sv = ift.get_signal_variance(lambda x: np.ones_like(x), hspace)
+    assert_equal(sv, 3.)
 
     space = ift.RGSpace(3, distances=1.)
     hspace = space.get_default_codomain()