Merge branch 'laplace_fixes' into 'master'

Laplace fixes

See merge request !189

nifty/operators/laplace_operator/laplace_operator.py

@@ -20,17 +20,17 @@ 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
 
 
 class LaplaceOperator(EndomorphicOperator):
     """A irregular LaplaceOperator with free boundary and excluding monopole.
 
-    This LaplaceOperator implements the second derivative of a Field in PowerSpace
-    on logarithmic or linear scale with vanishing curvature at the boundary, starting
-    at the second entry of the Field. The second derivative of the Field on the irregular grid
-    is calculated using finite differences.
+    This LaplaceOperator implements the second derivative of a Field in
+    PowerSpace  on logarithmic or linear scale with vanishing curvature at the
+    boundary, starting at the second entry of the Field. The second derivative
+    of the Field on the irregular grid is calculated using finite differences.
 
     Parameters
     ----------
@@ -50,15 +50,19 @@ class LaplaceOperator(EndomorphicOperator):
 
         self._logarithmic = bool(logarithmic)
 
+        pos = self.domain[0].kindex.copy()
         if self.logarithmic:
-            self.positions = self.domain[0].kindex.copy()
-            self.positions[1:] = np.log(self.positions[1:])
-            self.positions[0] = -1.
-        else:
-            self.positions = self.domain[0].kindex.copy()
-            self.positions[0] = -1
-
-        self.fwd_dist = self.positions[1:] - self.positions[:-1]
+            pos[1:] = np.log(pos[1:])
+            pos[0] = pos[1]-1.
+
+        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):
@@ -94,24 +98,20 @@ class LaplaceOperator(EndomorphicOperator):
         else:
             axes = x.domain_axes[spaces[0]]
         axis = axes[0]
+        nval = len(self._dposc)
         prefix = (slice(None),) * axis
-        fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape)
-        positions = self.positions.reshape((1,)*axis + self.positions.shape)
+        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()
-        x = x.val
-        ret[prefix + (slice(1, -1),)] = \
-            (-((x[prefix + (slice(1, -1),)] - x[prefix + (slice(0, -2),)]) /
-               fwd_dist[prefix + (slice(0, -1),)]) +
-             ((x[prefix + (slice(2, None),)] - x[prefix + (slice(1, -1),)]) /
-              fwd_dist[prefix + (slice(1, None),)]))
-        ret[prefix + (slice(1, -1),)] /= \
-            (positions[prefix + (slice(2, None),)] -
-             positions[prefix + (slice(None, -2),)])
-        ret *= 2.
-        ret[prefix + (slice(0, 2),)] = 0
-        ret[prefix + (slice(-1, -1),)] = 0
-        ret[prefix + (slice(2, None),)] *= \
-            np.sqrt(fwd_dist)[prefix + (slice(1, None),)]
+        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):
@@ -124,42 +124,19 @@ class LaplaceOperator(EndomorphicOperator):
         else:
             axes = x.domain_axes[spaces[0]]
         axis = axes[0]
+        nval = len(self._dposc)
         prefix = (slice(None),) * axis
-        fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape)
-        positions = self.positions.reshape((1,)*axis + self.positions.shape)
+        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[prefix + (slice(2, None),)] *= \
-            np.sqrt(fwd_dist)[prefix + (slice(1, None),)]
-        y[prefix + (slice(0, 2),)] = 0
-        y[prefix + (slice(-1, -1),)] = 0
-        ret = y.copy_empty()
-        y[prefix + (slice(1, -1),)] /= \
-            (positions[prefix + (slice(2, None),)] -
-             positions[prefix + (slice(None, -2),)])
-        y *= 2
-        ret[prefix + (slice(1, -1),)] = \
-            (-y[prefix + (slice(1, -1),)]/fwd_dist[prefix + (slice(0, -1),)] -
-             y[prefix + (slice(1, -1),)]/fwd_dist[prefix + (slice(1, None),)])
-        ret[prefix + (slice(0, -2),)] += \
-            y[prefix + (slice(1, -1),)] / fwd_dist[prefix + (slice(0, -1),)]
-        ret[prefix + (slice(2, None),)] += \
-            y[prefix + (slice(1, -1),)] / fwd_dist[prefix + (slice(1, None),)]
+        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)
-
-    def _irregular_laplace(self, x):
-        ret = np.zeros_like(x)
-        ret[1:-1] = (-(x[1:-1] - x[0:-2]) / self.fwd_dist[:-1] +
-                      (x[2:] - x[1:-1]) / self.fwd_dist[1:])
-        ret[1:-1] /= self.positions[2:] - self.positions[:-2]
-        ret *= 2.
-        return ret
-
-    def _irregular_adj_laplace(self, x):
-        ret = np.zeros_like(x)
-        y = x.copy()
-        y[1:-1] /= self.positions[2:] - self.positions[:-2]
-        y *= 2
-        ret[1:-1] = -y[1:-1] / self.fwd_dist[:-1] - y[1:-1] / self.fwd_dist[1:]
-        ret[0:-2] += y[1:-1] / self.fwd_dist[:-1]
-        ret[2:] += y[1:-1] / self.fwd_dist[1:]
-        return ret

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_allclose
+from itertools import product
+from test.common import expand
+
+
+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))