linear_operator.py 8.87 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
21
import numpy as np
22
23


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

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

    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.
48
    """
Theo Steininger's avatar
Theo Steininger committed
49

50
51
52
53
54
55
56
57
58
59
60
    # 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
61
    _ilog = (-1, 0, 1, -1, 2, -1, -1, -1, 3)
Martin Reinecke's avatar
Martin Reinecke committed
62
    _validMode = (False, True, True, False, True, False, False, False, True)
Martin Reinecke's avatar
Martin Reinecke committed
63
64
65
66
67
68
69
70
    _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
71
    _addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15)
Martin Reinecke's avatar
Martin Reinecke committed
72
    _backwards = 6
Martin Reinecke's avatar
Martin Reinecke committed
73
    _all_ops = 15
74

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

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

81
82
    def __init__(self):
        pass
83

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
102
103
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
104
105
106
107
        """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
108
        return self._flip_modes(self.INVERSE_BIT)
109

Martin Reinecke's avatar
Martin Reinecke committed
110
111
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
112
113
114
115
        """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
116
        return self._flip_modes(self.ADJOINT_BIT)
117

118
119
120
121
122
123
124
125
126
    @staticmethod
    def _toOperator(thing, dom):
        from .scaling_operator import ScalingOperator
        if isinstance(thing, LinearOperator):
            return thing
        if np.isscalar(thing):
            return ScalingOperator(thing, dom)
        return NotImplemented

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

134
135
136
137
138
139
140
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
        if np.isscalar(other) and other == 1.:
            return self
        other = self._toOperator(other, self.target)
        return ChainOperator.make([other, self])

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

148
149
150
    def __radd__(self, other):
        return self.__add__(other)

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

158
159
160
161
162
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
        return SumOperator.make([other, self], [False, True])

Martin Reinecke's avatar
Martin Reinecke committed
163
164
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
165
        """int : the supported operation modes
166

167
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
168
169
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
170
        """
Martin Reinecke's avatar
Martin Reinecke committed
171
        raise NotImplementedError
172

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

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

Martin Reinecke's avatar
Martin Reinecke committed
183
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
184
185
186
187
188
            - :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
189
190
191

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
192
        Field
193
194
195
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
196
197
198
        raise NotImplementedError

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

Martin Reinecke's avatar
Martin Reinecke committed
202
    def times(self, x):
203
204
205
206
207
208
209
210
211
        """ 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
212
        Field
213
214
            The processed Field living on the Operator's target domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
215
216
217
        return self.apply(x, self.TIMES)

    def inverse_times(self, x):
218
219
220
221
222
223
224
225
226
        """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
227
        Field
228
229
            The processed Field living on the Operator's domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
230
231
232
        return self.apply(x, self.INVERSE_TIMES)

    def adjoint_times(self, x):
233
234
235
236
237
238
239
240
241
        """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
242
        Field
243
244
            The processed Field living on the Operator's domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
245
        return self.apply(x, self.ADJOINT_TIMES)
246

Martin Reinecke's avatar
Martin Reinecke committed
247
    def adjoint_inverse_times(self, x):
248
249
250
251
252
253
254
255
256
        """ 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
257
        Field
258
259
260
261
262
263
264
            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
265
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
266

Martin Reinecke's avatar
Martin Reinecke committed
267
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
268
        """Same as :meth:`adjoint_inverse_times`"""
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
272
273
274
275
    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")
276

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