linear_operator.py 8.53 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

Martin Reinecke's avatar
Martin Reinecke committed
19
from __future__ import absolute_import, division, print_function
20

21
import numpy as np
22
23

from ..compat import *
Martin Reinecke's avatar
Martin Reinecke committed
24
from .operator import Operator
Philipp Arras's avatar
Fixups    
Philipp Arras committed
25

26

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

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

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

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

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

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

84
85
    def __init__(self):
        pass
86

87
    def _flip_modes(self, trafo):
Martin Reinecke's avatar
Martin Reinecke committed
88
        from .operator_adapter import OperatorAdapter
89
        return self if trafo == 0 else OperatorAdapter(self, trafo)
Martin Reinecke's avatar
Martin Reinecke committed
90

Martin Reinecke's avatar
Martin Reinecke committed
91
92
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
93
94
95
96
        """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
97
        return self._flip_modes(self.INVERSE_BIT)
98

Martin Reinecke's avatar
Martin Reinecke committed
99
100
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
101
102
103
104
        """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
105
        return self._flip_modes(self.ADJOINT_BIT)
106

107
    def __matmul__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
108
109
110
111
        if isinstance(other, LinearOperator):
            from .chain_operator import ChainOperator
            return ChainOperator.make([self, other])
        return Operator.__matmul__(self, other)
112
113

    def __rmatmul__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
114
        if isinstance(other, LinearOperator):
115
            from .chain_operator import ChainOperator
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
116
117
            return ChainOperator.make([other, self])
        return Operator.__rmatmul__(self, other)
118

Martin Reinecke's avatar
Martin Reinecke committed
119
    def __add__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
120
121
122
123
        if isinstance(other, LinearOperator):
            from .sum_operator import SumOperator
            return SumOperator.make([self, other], [False, False])
        return Operator.__add__(self, other)
124

125
126
127
    def __radd__(self, other):
        return self.__add__(other)

Martin Reinecke's avatar
Martin Reinecke committed
128
    def __sub__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
129
130
131
132
        if isinstance(other, LinearOperator):
            from .sum_operator import SumOperator
            return SumOperator.make([self, other], [False, True])
        return Operator.__sub__(self, other)
133

134
    def __rsub__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
135
136
137
138
        if isinstance(other, LinearOperator):
            from .sum_operator import SumOperator
            return SumOperator.make([other, self], [False, True])
        return Operator.__rsub__(self, other)
139

Martin Reinecke's avatar
Martin Reinecke committed
140
    @property
Martin Reinecke's avatar
Martin Reinecke committed
141
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
142
        """int : the supported operation modes
143

144
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
145
146
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
147
        """
Martin Reinecke's avatar
Martin Reinecke committed
148
        raise NotImplementedError
149

Martin Reinecke's avatar
Martin Reinecke committed
150
    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
151
        """ Applies the Operator to a given `x`, in a specified `mode`.
152
153
154
155
156
157
158

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

Martin Reinecke's avatar
Martin Reinecke committed
159
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
160
161
162
163
164
            - :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
165
166
167

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
168
        Field
169
170
171
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
172
173
174
        raise NotImplementedError

    def __call__(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
175
        """Same as :meth:`times`"""
Martin Reinecke's avatar
Martin Reinecke committed
176
        from ..field import Field
Martin Reinecke's avatar
Martin Reinecke committed
177
        from ..multi_field import MultiField
Martin Reinecke's avatar
Martin Reinecke committed
178
179
        if isinstance(x, (Field, MultiField)):
            return self.apply(x, self.TIMES)
Martin Reinecke's avatar
Martin Reinecke committed
180
        from ..linearization import Linearization
Martin Reinecke's avatar
Martin Reinecke committed
181
        if isinstance(x, Linearization):
Martin Reinecke's avatar
Martin Reinecke committed
182
183
            return Linearization(self(x._val), self(x._jac))
        return self.__matmul__(x)
184

Martin Reinecke's avatar
Martin Reinecke committed
185
    def times(self, x):
186
187
188
189
190
191
192
193
194
        """ 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
195
        Field
196
197
            The processed Field living on the Operator's target domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
198
199
200
        return self.apply(x, self.TIMES)

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

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

Martin Reinecke's avatar
Martin Reinecke committed
230
    def adjoint_inverse_times(self, x):
231
232
233
234
235
236
237
238
239
        """ 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
240
        Field
241
242
243
244
245
246
247
            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
248
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
249

Martin Reinecke's avatar
Martin Reinecke committed
250
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
251
        """Same as :meth:`adjoint_inverse_times`"""
Martin Reinecke's avatar
Martin Reinecke committed
252
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
253

Martin Reinecke's avatar
Martin Reinecke committed
254
255
    def _check_mode(self, mode):
        if not self._validMode[mode]:
Martin Reinecke's avatar
Martin Reinecke committed
256
            raise NotImplementedError("invalid operator mode specified")
Martin Reinecke's avatar
Martin Reinecke committed
257
        if mode & self.capability == 0:
Martin Reinecke's avatar
Martin Reinecke committed
258
259
            raise NotImplementedError(
                "requested operator mode is not supported")
260

Martin Reinecke's avatar
Martin Reinecke committed
261
262
    def _check_input(self, x, mode):
        self._check_mode(mode)
Martin Reinecke's avatar
Martin Reinecke committed
263
        if self._dom(mode) is not x.domain:
Martin Reinecke's avatar
Martin Reinecke committed
264
            raise ValueError("The operator's and field's domains don't match.")