Commit 05eeeabb authored by Theo Steininger's avatar Theo Steininger

LaplaceOperator -> pep8. Moved to nifty.operators

parent 16164871
Pipeline #14528 failed with stage
in 4 minutes and 46 seconds
......@@ -54,8 +54,8 @@ from operators import *
from probing import *
from library import *
from sugar import *
import library
import plotting
from wiener_filter_curvature import WienerFilterCurvature
from critical_power_curvature import CriticalPowerCurvature
from laplace_operator import LaplaceOperator
from smoothness_operator import SmoothnessOperator
\ No newline at end of file
......@@ -7,8 +7,8 @@ class WienerFilterCurvature(InvertibleOperatorMixin, EndomorphicOperator):
This operator implements the second derivative of the
WienerFilterEnergy used in some minimization algorithms or
for error estimates of the posterior maps. It corresponds to the
inverse propagator.
for error estimates of the posterior maps. It is the
inverse of the propagator operator.
Parameters
......
......@@ -32,10 +32,8 @@ from invertible_operator_mixin import InvertibleOperatorMixin
from projection_operator import ProjectionOperator
from propagator_operator import PropagatorOperator
from propagator_operator import HarmonicPropagatorOperator
from composed_operator import ComposedOperator
from response_operator import ResponseOperator
from laplace_operator import LaplaceOperator
# -*- coding: utf-8 -*-
from laplace_operator import LaplaceOperator
# 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 numpy as np
from nifty import Field,\
......@@ -5,21 +22,6 @@ from nifty import Field,\
PowerSpace
import nifty.nifty_utilities as utilities
# def _irregular_nabla(x,k):
# #applies forward differences and does nonesense at the edge. Thus needs cutting
# y = -x
# y[:-1] += x[1:]
# y[1:-1] /= - k[1:-1] + k[2:]
# return y
#
# def _irregular_adj_nabla(z, k):
# #applies backwards differences*(-1) and does nonesense at the edge. Thus needs cutting
# x = z.copy()
# x[1:-1] /= - k[1:-1] + k[2:]
# y = -x
# y[1:] += x[:-1]
# return y
class LaplaceOperator(EndomorphicOperator):
"""A irregular LaplaceOperator with free boundary and excluding monopole.
......@@ -36,18 +38,20 @@ class LaplaceOperator(EndomorphicOperator):
default : True
"""
def __init__(self, domain,
default_spaces=None, logarithmic = True):
def __init__(self, domain, default_spaces=None, logarithmic=True):
super(LaplaceOperator, self).__init__(default_spaces)
if (domain is not None):
if (not isinstance(domain, PowerSpace)):
raise TypeError("The domain has to live over a PowerSpace")
self._domain = self._parse_domain(domain)
if logarithmic :
if len(self.domain) != 0:
raise ValueError("The domain must contain exactly one PowerSpace.")
if not isinstance(self.domain[0], PowerSpace):
raise TypeError("The domain must contain exactly one PowerSpace.")
if logarithmic:
self.positions = self.domain[0].kindex.copy()
self.positions[1:] = np.log(self.positions[1:])
self.positions[0] = -1.
else :
else:
self.positions = self.domain[0].kindex.copy()
self.positions[0] = -1
......@@ -73,7 +77,6 @@ class LaplaceOperator(EndomorphicOperator):
def self_adjoint(self):
return False
def _times(self, x, spaces):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
if spaces is None:
......@@ -89,15 +92,20 @@ class LaplaceOperator(EndomorphicOperator):
positions = self.positions.reshape((1,)*axis + self.positions.shape)
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[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),)]
return Field(self.domain, val=ret).weight(power=-0.5,spaces=spaces)
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),)]
return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces)
def _adjoint_times(self, x, spaces):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
......@@ -113,23 +121,28 @@ class LaplaceOperator(EndomorphicOperator):
fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape)
positions = self.positions.reshape((1,)*axis + self.positions.shape)
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[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[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),)]
return Field(self.domain, val=ret).weight(-1,spaces=spaces)
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),)]
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] = (-(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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment