Merge branch 'improve_FuncConvolutionOperator' into 'NIFTy_5'

FuncConvolutionOperator with mean subtraction

See merge request !332

nifty5/operators/convolution_operators.py

@@ -77,4 +77,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.
 

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, assert_equal
+
+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):
+    ones = lambda x: np.ones(x.shape)
+    sig = ift.Field.from_random('normal', domain=space)
+    fco_op = ift.FuncConvolutionOperator(space, ones)
+    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)