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

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

24

Martin Reinecke's avatar
Martin Reinecke committed
25
class LinearOperator(Operator):
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
    def _flip_modes(self, trafo):
Martin Reinecke's avatar
Martin Reinecke committed
83
        from .operator_adapter import OperatorAdapter
84
        return self if trafo == 0 else OperatorAdapter(self, trafo)
Martin Reinecke's avatar
Martin Reinecke committed
85

Martin Reinecke's avatar
Martin Reinecke committed
86
87
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
88
89
90
91
        """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
92
        return self._flip_modes(self.INVERSE_BIT)
93

Martin Reinecke's avatar
Martin Reinecke committed
94
95
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
96
97
98
99
        """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
100
        return self._flip_modes(self.ADJOINT_BIT)
101

102
    def __matmul__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
103
104
105
106
        if isinstance(other, LinearOperator):
            from .chain_operator import ChainOperator
            return ChainOperator.make([self, other])
        return Operator.__matmul__(self, other)
107
108

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

Martin Reinecke's avatar
Martin Reinecke committed
114
    def _myadd(self, other, oneg):
Martin Reinecke's avatar
Martin Reinecke committed
115
116
        from .sum_operator import SumOperator
        return SumOperator.make((self, other), (False, oneg))
Martin Reinecke's avatar
Martin Reinecke committed
117

Martin Reinecke's avatar
Martin Reinecke committed
118
    def __add__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
119
        if isinstance(other, LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
120
            return self._myadd(other, False)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
121
        return Operator.__add__(self, other)
122

123
124
125
    def __radd__(self, other):
        return self.__add__(other)

Martin Reinecke's avatar
Martin Reinecke committed
126
    def __sub__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
127
        if isinstance(other, LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
128
            return self._myadd(other, True)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
129
        return Operator.__sub__(self, other)
130

131
    def __rsub__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
132
        if isinstance(other, LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
133
            return other._myadd(self, True)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
134
        return Operator.__rsub__(self, other)
135

Martin Reinecke's avatar
Martin Reinecke committed
136
    @property
Martin Reinecke's avatar
Martin Reinecke committed
137
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
138
        """int : the supported operation modes
139

140
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
141
142
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
143
        """
Martin Reinecke's avatar
Martin Reinecke committed
144
        return self._capability
145

146
147
148
    def force(self, x):
        return self.apply(x.extract(self.domain), self.TIMES)

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

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

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

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

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

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

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
260
261
    def _check_input(self, x, mode):
        self._check_mode(mode)
Martin Reinecke's avatar
Martin Reinecke committed
262
        self._check_domain_equality(self._dom(mode), x.domain)