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
Philipp Arras's avatar
Fixups    
Philipp Arras committed
20

21
import numpy as np
22

Philipp Arras's avatar
Fixups    
Philipp Arras committed
23
24
from ..utilities import NiftyMetaBase

25

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

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

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

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

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

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

83
84
    def __init__(self):
        pass
85

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

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

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

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

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

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

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

120
121
122
123
124
125
126
127
128
    @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
129
130
    def __mul__(self, other):
        from .chain_operator import ChainOperator
131
132
133
        if np.isscalar(other) and other == 1.:
            return self
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
134
        return ChainOperator.make([self, other])
135

136
137
138
139
140
141
142
    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
143
144
    def __add__(self, other):
        from .sum_operator import SumOperator
145
146
147
        if np.isscalar(other) and other == 0.:
            return self
        other = self._toOperator(other, self.domain)
148
        return SumOperator.make([self, other], [False, False])
149

150
151
152
    def __radd__(self, other):
        return self.__add__(other)

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

160
161
162
163
164
    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
165
166
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
167
        """int : the supported operation modes
168

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

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

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

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

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

    def __call__(self, x):
Philipp Arras's avatar
Fixups    
Philipp Arras committed
201
        from ..nonlinear_operators import LinearModel, NonlinearOperator
Martin Reinecke's avatar
Martin Reinecke committed
202
        """Same as :meth:`times`"""
Philipp Arras's avatar
Fixups    
Philipp Arras committed
203
204
        if isinstance(x, NonlinearOperator):
            return LinearModel(x, self)
Martin Reinecke's avatar
Martin Reinecke committed
205
        return self.apply(x, self.TIMES)
206

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
276
277
278
279
280
    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")
281

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