revise Laplace operator and add some basic testing

nifty/operators/laplace_operator/laplace_operator.py

@@ -20,7 +20,7 @@ import numpy as np
 from nifty.field import Field
 from nifty.spaces.power_space import PowerSpace
 from nifty.operators.endomorphic_operator import EndomorphicOperator
-
+from nifty import sqrt
 import nifty.nifty_utilities as utilities
 
 
@@ -55,9 +55,14 @@ class LaplaceOperator(EndomorphicOperator):
             pos[1:] = np.log(pos[1:])
             pos[0] = pos[1]-1.
 
-        self._dist_l = pos[1:-1]-pos[:-2]
-        self._dist_r = pos[2:]-pos[1:-1]
-        self._dist_c = 0.5*(pos[2:]-pos[:-2])
+        self._dpos = pos[1:]-pos[:-1]  # defined between points
+        # centered distances (also has entries for the first and last point
+        # for convenience, but they will never affect the result)
+        self._dposc = np.empty_like(pos)
+        self._dposc[:-1] = self._dpos
+        self._dposc[-1]=0.
+        self._dposc[1:]+=self._dpos
+        self._dposc*=0.5
 
     @property
     def target(self):
@@ -93,20 +98,20 @@ class LaplaceOperator(EndomorphicOperator):
         else:
             axes = x.domain_axes[spaces[0]]
         axis = axes[0]
-        nval = len(self._dist_l)
+        nval = len(self._dposc)
         prefix = (slice(None),) * axis
-        sl_c = prefix + slice(1,-1)
-        sl_l = prefix + slice(None,-2)
-        sl_r = prefix + slice(2,None)
-        dist_l = self._dist_l.reshape((1,)*axis + (nval,))
-        dist_r = self._dist_r.reshape((1,)*axis + (nval,))
-        dist_c = self._dist_c.reshape((1,)*axis + (nval,))
-        dx_r = x[sl_r] - x[sl_c]
-        dx_l = x[sl_c] - x[sl_l]
+        sl_l = prefix + (slice(None,-1),)  # "left" slice
+        sl_r = prefix + (slice(1,None),)  # "right" slice
+        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
+        dposc = self._dposc.reshape((1,)*axis + (nval,))
+        deriv = (x.val[sl_r]-x.val[sl_l])/dpos  # defined between points
         ret = x.val.copy_empty()
-        ret[sl_c] = (dx_r/dist_r - dx_l/dist_l)/sqrt(dist_c)
-        ret[prefix + slice(None,2)] = 0.
-        res[prefix + slice(-1,None)] = 0.
+        ret[sl_l] = deriv
+        ret[prefix + (-1,)] = 0.
+        ret[sl_r] -= deriv
+        ret /= sqrt(dposc)
+        ret[prefix + (slice(None,2),)] = 0.
+        ret[prefix + (-1,)] = 0.
         return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces)
 
     def _adjoint_times(self, x, spaces):
@@ -119,32 +124,19 @@ class LaplaceOperator(EndomorphicOperator):
         else:
             axes = x.domain_axes[spaces[0]]
         axis = axes[0]
-        nval = len(self._dist_l)
+        nval = len(self._dposc)
         prefix = (slice(None),) * axis
-        sl_c = prefix + slice(1,-1)
-        sl_l = prefix + slice(None,-2)
-        sl_r = prefix + slice(2,None)
-        dist_l = self._dist_l.reshape((1,)*axis + (nval,))
-        dist_r = self._dist_r.reshape((1,)*axis + (nval,))
-        dist_c = self._dist_c.reshape((1,)*axis + (nval,))
-        dx_r = x[sl_r] - x[sl_c]
-        dx_l = x[sl_c] - x[sl_l]
+        sl_l = prefix + (slice(None,-1),)  # "left" slice
+        sl_r = prefix + (slice(1,None),)  # "right" slice
+        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
+        dposc = self._dposc.reshape((1,)*axis + (nval,))
         y = x.copy().weight(power=0.5).val
-        y[sl_c] *= sqrt(dist_c)
-        y[prefix + slice(None, 2)] = 0.
-        y[prefix + slice(-1, None)] = 0.
-        ret = y.copy_empty()
-        y[sl_c] /= dist_c
-        ret[sl_c] = (y[sl_r]-y[sl_c])/dist_r - (y[sl_c]-y[sl_l])/dist_l
-
-        ret[prefix + (0,)] = y[prefix+(1,)] / dist_l[prefix+(0,)]
-        ret[prefix + (-1,)] = y[prefix+(-2,)] / dist_r[prefix + (-1,)]
+        y /= sqrt(dposc)
+        y[prefix + (slice(None, 2),)] = 0.
+        y[prefix + (-1,)] = 0.
+        deriv = (y[sl_r]-y[sl_l])/dpos  # defined between points
+        ret = x.val.copy_empty()
+        ret[sl_l] = deriv
+        ret[prefix + (-1,)] = 0.
+        ret[sl_r] -= deriv
         return Field(self.domain, val=ret).weight(-1, spaces=spaces)
-
-Laplace:
-L = (dxr/dr - dxl/dl)/  sqrt(dc)
-
-adjoint Laplace:
-
-tmp = x/sqrt(dc)
-tmp2 = (tmpr-tmp)/dr - (tmp-tmpl)/dl

test/test_operators/test_laplace_operator.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-2017 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
+# and financially supported by the Studienstiftung des deutschen Volkes.
+
+import unittest
+import numpy as np
+import nifty as ift
+from numpy.testing import assert_equal,\
+    assert_allclose
+from itertools import product
+from test.common import expand
+from nose.plugins.skip import SkipTest
+
+
+class LaplaceOperatorTests(unittest.TestCase):
+    @expand(product([None,False,True], [False,True], [10,100,1000]))
+    def test_Laplace(self, log1, log2, sz):
+        s=ift.RGSpace(sz,harmonic=True)
+        p=ift.PowerSpace(s,logarithmic=log1)
+        L=ift.LaplaceOperator(p,logarithmic=log2)
+        arr=np.random.random(p.shape[0])
+        fp=ift.Field(p,val=arr)
+        assert_allclose (L(fp).vdot(L(fp)),L.adjoint_times(L(fp)).vdot(fp))