laplace_operator.py 5.21 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# 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.
18 19

import numpy as np
20 21 22 23 24
from ...field import Field
from ...spaces.power_space import PowerSpace
from ..endomorphic_operator import EndomorphicOperator
from ... import sqrt
from ... import nifty_utilities as utilities
Jakob Knollmueller's avatar
Jakob Knollmueller committed
25 26


27
class LaplaceOperator(EndomorphicOperator):
Jakob Knollmueller's avatar
Jakob Knollmueller committed
28 29
    """A irregular LaplaceOperator with free boundary and excluding monopole.

Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
30 31 32 33
    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.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
34 35 36 37 38 39 40 41

    Parameters
    ----------
    logarithmic : boolean,
        Whether smoothness is calculated on a logarithmic scale or linear scale
        default : True
    """

42
    def __init__(self, domain, default_spaces=None, logarithmic=True):
43
        super(LaplaceOperator, self).__init__(default_spaces)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
44
        self._domain = self._parse_domain(domain)
Martin Reinecke's avatar
Martin Reinecke committed
45
        if len(self.domain) != 1:
46 47 48 49 50
            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.")

51 52
        self._logarithmic = bool(logarithmic)

Martin Reinecke's avatar
Martin Reinecke committed
53
        pos = self.domain[0].kindex.copy()
54
        if self.logarithmic:
Martin Reinecke's avatar
Martin Reinecke committed
55 56
            pos[1:] = np.log(pos[1:])
            pos[0] = pos[1]-1.
57

58 59 60 61 62
        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
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
63 64 65
        self._dposc[-1] = 0.
        self._dposc[1:] += self._dpos
        self._dposc *= 0.5
Jakob Knollmueller's avatar
Jakob Knollmueller committed
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86

    @property
    def target(self):
        return self._domain

    @property
    def domain(self):
        return self._domain

    @property
    def unitary(self):
        return False

    @property
    def symmetric(self):
        return False

    @property
    def self_adjoint(self):
        return False

87 88 89 90
    @property
    def logarithmic(self):
        return self._logarithmic

91 92 93 94 95 96 97 98 99 100
    def _times(self, x, spaces):
        spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
        if spaces is None:
            # this case means that x lives on only one space, which is
            # identical to the space in the domain of `self`. Otherwise the
            # input check of LinearOperator would have failed.
            axes = x.domain_axes[0]
        else:
            axes = x.domain_axes[spaces[0]]
        axis = axes[0]
101
        nval = len(self._dposc)
102
        prefix = (slice(None),) * axis
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
103 104
        sl_l = prefix + (slice(None, -1),)  # "left" slice
        sl_r = prefix + (slice(1, None),)  # "right" slice
105 106 107
        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
Martin Reinecke's avatar
more  
Martin Reinecke committed
108
        ret = np.empty_like(x.val)
109 110 111 112
        ret[sl_l] = deriv
        ret[prefix + (-1,)] = 0.
        ret[sl_r] -= deriv
        ret /= sqrt(dposc)
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
113
        ret[prefix + (slice(None, 2),)] = 0.
114
        ret[prefix + (-1,)] = 0.
115
        return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces)
116 117 118 119 120 121 122 123 124 125 126

    def _adjoint_times(self, x, spaces):
        spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
        if spaces is None:
            # this case means that x lives on only one space, which is
            # identical to the space in the domain of `self`. Otherwise the
            # input check of LinearOperator would have failed.
            axes = x.domain_axes[0]
        else:
            axes = x.domain_axes[spaces[0]]
        axis = axes[0]
127
        nval = len(self._dposc)
128
        prefix = (slice(None),) * axis
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
129 130
        sl_l = prefix + (slice(None, -1),)  # "left" slice
        sl_r = prefix + (slice(1, None),)  # "right" slice
131 132
        dpos = self._dpos.reshape((1,)*axis + (nval-1,))
        dposc = self._dposc.reshape((1,)*axis + (nval,))
133
        y = x.copy().weight(power=0.5).val
134 135 136 137
        y /= sqrt(dposc)
        y[prefix + (slice(None, 2),)] = 0.
        y[prefix + (-1,)] = 0.
        deriv = (y[sl_r]-y[sl_l])/dpos  # defined between points
Martin Reinecke's avatar
more  
Martin Reinecke committed
138
        ret = np.empty_like(x.val)
139 140 141
        ret[sl_l] = deriv
        ret[prefix + (-1,)] = 0.
        ret[sl_r] -= deriv
142
        return Field(self.domain, val=ret).weight(-1, spaces=spaces)