linear_operator.py 9.03 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
# 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/>.
Theo Steininger's avatar
Theo Steininger committed
13
#
14
# Copyright(C) 2013-2018 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15 16 17
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

19
import abc
Martin Reinecke's avatar
Martin Reinecke committed
20
from ..utilities import NiftyMetaBase
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..field import Field
22
import numpy as np
23 24


Martin Reinecke's avatar
Martin Reinecke committed
25
class LinearOperator(NiftyMetaBase()):
26 27 28 29
    """NIFTY base class for linear operators.

    The base NIFTY operator class is an abstract class from which
    other specific operator subclasses are derived.
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

    Attributes
    ----------
    TIMES : int
        Symbolic constant representing normal operator application
    ADJOINT_TIMES : int
        Symbolic constant representing adjoint operator application
    INVERSE_TIMES : int
        Symbolic constant representing inverse operator application
    ADJOINT_INVERSE_TIMES : int
        Symbolic constant representing adjoint inverse operator application
    INVERSE_ADJOINT_TIMES : int
        same as ADJOINT_INVERSE_TIMES

    Notes
    -----
    The symbolic constants for the operation modes can be combined by the
    "bitwise-or" operator "|", for expressing the capability of the operator
    by means of a single integer number.
49
    """
Theo Steininger's avatar
Theo Steininger committed
50

51 52 53 54 55 56 57 58 59 60 61
    # Field Operator Modes
    TIMES = 1
    ADJOINT_TIMES = 2
    INVERSE_TIMES = 4
    ADJOINT_INVERSE_TIMES = 8
    INVERSE_ADJOINT_TIMES = 8

    # Operator Transform Flags
    ADJOINT_BIT = 1
    INVERSE_BIT = 2

Martin Reinecke's avatar
Martin Reinecke committed
62
    _ilog = (-1, 0, 1, -1, 2, -1, -1, -1, 3)
Martin Reinecke's avatar
Martin Reinecke committed
63
    _validMode = (False, True, True, False, True, False, False, False, True)
Martin Reinecke's avatar
Martin Reinecke committed
64 65 66 67 68 69 70 71
    _modeTable = ((1, 2, 4, 8),
                  (2, 1, 8, 4),
                  (4, 8, 1, 2),
                  (8, 4, 2, 1))
    _capTable = ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15),
                 (0, 2, 1, 3, 8, 10, 9, 11, 4, 6, 5, 7, 12, 14, 13, 15),
                 (0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15),
                 (0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15))
Martin Reinecke's avatar
Martin Reinecke committed
72
    _addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15)
Martin Reinecke's avatar
Martin Reinecke committed
73
    _backwards = 6
Martin Reinecke's avatar
Martin Reinecke committed
74
    _all_ops = 15
75

Martin Reinecke's avatar
Martin Reinecke committed
76 77
    def _dom(self, mode):
        return self.domain if (mode & 9) else self.target
78

Martin Reinecke's avatar
Martin Reinecke committed
79 80
    def _tgt(self, mode):
        return self.domain if (mode & 6) else self.target
81

82 83
    def __init__(self):
        pass
84

85
    @abc.abstractproperty
86
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
87 88 89
        """DomainTuple : the operator's input domain

            The domain on which the Operator's input Field lives."""
90
        raise NotImplementedError
91 92 93

    @abc.abstractproperty
    def target(self):
Martin Reinecke's avatar
Martin Reinecke committed
94 95 96
        """DomainTuple : the operator's output domain

            The domain on which the Operator's output Field lives."""
97 98
        raise NotImplementedError

99
    def _flip_modes(self, trafo):
Martin Reinecke's avatar
Martin Reinecke committed
100
        from .operator_adapter import OperatorAdapter
101
        return self if trafo == 0 else OperatorAdapter(self, trafo)
Martin Reinecke's avatar
Martin Reinecke committed
102

Martin Reinecke's avatar
Martin Reinecke committed
103 104
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
105 106 107 108
        """LinearOperator : the inverse of `self`

        Returns a LinearOperator object which behaves as if it were
        the inverse of this operator."""
Martin Reinecke's avatar
Martin Reinecke committed
109
        return self._flip_modes(self.INVERSE_BIT)
110

Martin Reinecke's avatar
Martin Reinecke committed
111 112
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
113 114 115 116
        """LinearOperator : the adjoint of `self`

        Returns a LinearOperator object which behaves as if it were
        the adjoint of this operator."""
Martin Reinecke's avatar
Martin Reinecke committed
117
        return self._flip_modes(self.ADJOINT_BIT)
118

119 120 121 122 123 124 125 126 127 128
    @staticmethod
    def _toOperator(thing, dom):
        from .diagonal_operator import DiagonalOperator
        from .scaling_operator import ScalingOperator
        if isinstance(thing, LinearOperator):
            return thing
        if isinstance(thing, Field):
            return DiagonalOperator(thing)
        if np.isscalar(thing):
            return ScalingOperator(thing, dom)
Martin Reinecke's avatar
Martin Reinecke committed
129
        return NotImplemented
130

Martin Reinecke's avatar
Martin Reinecke committed
131 132
    def __mul__(self, other):
        from .chain_operator import ChainOperator
133 134
        if np.isscalar(other) and other == 1.:
            return self
135
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
136
        return ChainOperator.make([self, other])
137

138 139
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
140 141
        if np.isscalar(other) and other == 1.:
            return self
142
        other = self._toOperator(other, self.target)
Martin Reinecke's avatar
Martin Reinecke committed
143
        return ChainOperator.make([other, self])
144

Martin Reinecke's avatar
Martin Reinecke committed
145 146
    def __add__(self, other):
        from .sum_operator import SumOperator
147 148
        if np.isscalar(other) and other == 0.:
            return self
149
        other = self._toOperator(other, self.domain)
150
        return SumOperator.make([self, other], [False, False])
151

152 153 154
    def __radd__(self, other):
        return self.__add__(other)

Martin Reinecke's avatar
Martin Reinecke committed
155 156
    def __sub__(self, other):
        from .sum_operator import SumOperator
157 158
        if np.isscalar(other) and other == 0.:
            return self
159
        other = self._toOperator(other, self.domain)
160
        return SumOperator.make([self, other], [False, True])
161

162 163 164
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
165
        return SumOperator.make([other, self], [False, True])
166

Martin Reinecke's avatar
Martin Reinecke committed
167 168
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
169
        """int : the supported operation modes
170

171
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
172 173
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
174
        """
Martin Reinecke's avatar
Martin Reinecke committed
175
        raise NotImplementedError
176

Martin Reinecke's avatar
Martin Reinecke committed
177 178
    @abc.abstractmethod
    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
179
        """ Applies the Operator to a given `x`, in a specified `mode`.
180 181 182 183 184 185 186

        Parameters
        ----------
        x : Field
            The input Field, living on the Operator's domain or target,
            depending on mode.

Martin Reinecke's avatar
Martin Reinecke committed
187
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
188 189 190 191 192
            - :attr:`TIMES`: normal application
            - :attr:`ADJOINT_TIMES`: adjoint application
            - :attr:`INVERSE_TIMES`: inverse application
            - :attr:`ADJOINT_INVERSE_TIMES` or
              :attr:`INVERSE_ADJOINT_TIMES`: adjoint inverse application
193 194 195

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
196
        Field
197 198 199
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
200 201 202
        raise NotImplementedError

    def __call__(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
203
        """Same as :meth:`times`"""
Martin Reinecke's avatar
Martin Reinecke committed
204
        return self.apply(x, self.TIMES)
205

Martin Reinecke's avatar
Martin Reinecke committed
206
    def times(self, x):
207 208 209 210 211 212 213 214 215
        """ Applies the Operator to a given Field.

        Parameters
        ----------
        x : Field
            The input Field, living on the Operator's domain.

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
216
        Field
217 218
            The processed Field living on the Operator's target domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
219 220 221
        return self.apply(x, self.TIMES)

    def inverse_times(self, x):
222 223 224 225 226 227 228 229 230
        """Applies the inverse Operator to a given Field.

        Parameters
        ----------
        x : Field
            The input Field, living on the Operator's target domain

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
231
        Field
232 233
            The processed Field living on the Operator's domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
234 235 236
        return self.apply(x, self.INVERSE_TIMES)

    def adjoint_times(self, x):
237 238 239 240 241 242 243 244 245
        """Applies the adjoint-Operator to a given Field.

        Parameters
        ----------
        x : Field
            The input Field, living on the Operator's target domain

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
246
        Field
247 248
            The processed Field living on the Operator's domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
249
        return self.apply(x, self.ADJOINT_TIMES)
250

Martin Reinecke's avatar
Martin Reinecke committed
251
    def adjoint_inverse_times(self, x):
252 253 254 255 256 257 258 259 260
        """ Applies the adjoint-inverse Operator to a given Field.

        Parameters
        ----------
        x : Field
            The input Field, living on the Operator's domain.

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
261
        Field
262 263 264 265 266 267 268
            The processed Field living on the Operator's target domain.

        Notes
        -----
        If the operator has an `inverse` then the inverse adjoint is identical
        to the adjoint inverse. We provide both names for convenience.
        """
Martin Reinecke's avatar
Martin Reinecke committed
269
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
270

Martin Reinecke's avatar
Martin Reinecke committed
271
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
272
        """Same as :meth:`adjoint_inverse_times`"""
Martin Reinecke's avatar
Martin Reinecke committed
273
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
274

Martin Reinecke's avatar
Martin Reinecke committed
275 276 277 278 279
    def _check_mode(self, mode):
        if not self._validMode[mode]:
            raise ValueError("invalid operator mode specified")
        if mode & self.capability == 0:
            raise ValueError("requested operator mode is not supported")
280

Martin Reinecke's avatar
Martin Reinecke committed
281 282 283
    def _check_input(self, x, mode):
        self._check_mode(mode)
        if x.domain != self._dom(mode):
Martin Reinecke's avatar
Martin Reinecke committed
284
            raise ValueError("The operator's and field's domains don't match.")