diff --git a/nifty/operators/laplace_operator/laplace_operator.py b/nifty/operators/laplace_operator/laplace_operator.py index 00f6eba8fca1464454e08b8cdf8c0ec456a66086..ccefae9a8047176fc1abd19b33a56eb7c845a320 100644 --- a/nifty/operators/laplace_operator/laplace_operator.py +++ b/nifty/operators/laplace_operator/laplace_operator.py @@ -50,15 +50,14 @@ 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 + pos[1:] = np.log(pos[1:]) + pos[0] = pos[1]-1. - self.fwd_dist = self.positions[1:] - self.positions[:-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]) @property def target(self): @@ -94,24 +93,20 @@ class LaplaceOperator(EndomorphicOperator): else: axes = x.domain_axes[spaces[0]] axis = axes[0] + nval = len(self._dist_l) 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_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] 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_c] = (dx_r/dist_r - dx_l/dist_l)/sqrt(dist_c) + ret[prefix + slice(None,2)] = 0. + res[prefix + slice(-1,None)] = 0. return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces) def _adjoint_times(self, x, spaces): @@ -124,42 +119,32 @@ class LaplaceOperator(EndomorphicOperator): else: axes = x.domain_axes[spaces[0]] axis = axes[0] + nval = len(self._dist_l) 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_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] 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 + y[sl_c] *= sqrt(dist_c) + y[prefix + slice(None, 2)] = 0. + y[prefix + slice(-1, None)] = 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[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,)] 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 +Laplace: +L = (dxr/dr - dxl/dl)/ sqrt(dc) + +adjoint Laplace: + +tmp = x/sqrt(dc) +tmp2 = (tmpr-tmp)/dr - (tmp-tmpl)/dl