linear_operator.py 8.26 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-2019 Max-Planck-Society
Theo Steininger's avatar
Theo Steininger committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17

Martin Reinecke's avatar
Martin Reinecke committed
18
from .operator import Operator
Philipp Arras's avatar
Fixups    
Philipp Arras committed
19

20

Martin Reinecke's avatar
Martin Reinecke committed
21
class LinearOperator(Operator):
22
23
24
25
    """NIFTY base class for linear operators.

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
72
73
    def _dom(self, mode):
        return self.domain if (mode & 9) else self.target
74

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

78
    def _flip_modes(self, trafo):
Martin Reinecke's avatar
Martin Reinecke committed
79
        from .operator_adapter import OperatorAdapter
80
        return self if trafo == 0 else OperatorAdapter(self, trafo)
Martin Reinecke's avatar
Martin Reinecke committed
81

Martin Reinecke's avatar
Martin Reinecke committed
82
83
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
84
85
86
87
        """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
88
        return self._flip_modes(self.INVERSE_BIT)
89

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

98
    def __matmul__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
99
100
101
102
        if isinstance(other, LinearOperator):
            from .chain_operator import ChainOperator
            return ChainOperator.make([self, other])
        return Operator.__matmul__(self, other)
103
104

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

Martin Reinecke's avatar
Martin Reinecke committed
110
    def _myadd(self, other, oneg):
Martin Reinecke's avatar
Martin Reinecke committed
111
112
        from .sum_operator import SumOperator
        return SumOperator.make((self, other), (False, oneg))
Martin Reinecke's avatar
Martin Reinecke committed
113

Martin Reinecke's avatar
Martin Reinecke committed
114
    def __add__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
115
        if isinstance(other, LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
116
            return self._myadd(other, False)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
117
        return Operator.__add__(self, other)
118

119
120
121
    def __radd__(self, other):
        return self.__add__(other)

Martin Reinecke's avatar
Martin Reinecke committed
122
    def __sub__(self, other):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
123
        if isinstance(other, LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
124
            return self._myadd(other, True)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
125
        return Operator.__sub__(self, other)
126

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

Martin Reinecke's avatar
Martin Reinecke committed
132
    @property
Martin Reinecke's avatar
Martin Reinecke committed
133
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
134
        """int : the supported operation modes
135

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

142
143
144
    def force(self, x):
        return self.apply(x.extract(self.domain), self.TIMES)

Martin Reinecke's avatar
Martin Reinecke committed
145
    def apply(self, x, mode):
146
        """Applies the Operator to a given `x`, in a specified `mode`.
147
148
149
150

        Parameters
        ----------
        x : Field
Philipp Arras's avatar
Philipp Arras committed
151
            The input Field, defined on the Operator's domain or target,
152
153
            depending on mode.

Martin Reinecke's avatar
Martin Reinecke committed
154
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
155
156
157
158
159
            - :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
160
161
162

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
163
        Field
Philipp Arras's avatar
Philipp Arras committed
164
            The processed Field defined on the Operator's target or domain,
165
166
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
167
168
169
        raise NotImplementedError

    def __call__(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
170
        """Same as :meth:`times`"""
Martin Reinecke's avatar
Martin Reinecke committed
171
        from ..field import Field
Martin Reinecke's avatar
Martin Reinecke committed
172
        from ..multi_field import MultiField
Martin Reinecke's avatar
Martin Reinecke committed
173
174
        if isinstance(x, (Field, MultiField)):
            return self.apply(x, self.TIMES)
Martin Reinecke's avatar
Martin Reinecke committed
175
        from ..linearization import Linearization
Martin Reinecke's avatar
Martin Reinecke committed
176
        if isinstance(x, Linearization):
Philipp Arras's avatar
Philipp Arras committed
177
            return x.new(self(x._val), self).prepend_jac(x.jac)
178
        return self@x
179

Martin Reinecke's avatar
Martin Reinecke committed
180
    def times(self, x):
181
        """Applies the Operator to a given Field.
182
183
184
185

        Parameters
        ----------
        x : Field
Philipp Arras's avatar
Philipp Arras committed
186
            The input Field, defined on the Operator's domain.
187
188
189

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
190
        Field
Philipp Arras's avatar
Philipp Arras committed
191
            The processed Field defined on the Operator's target domain.
192
        """
Martin Reinecke's avatar
Martin Reinecke committed
193
194
195
        return self.apply(x, self.TIMES)

    def inverse_times(self, x):
196
197
198
199
200
        """Applies the inverse Operator to a given Field.

        Parameters
        ----------
        x : Field
Philipp Arras's avatar
Philipp Arras committed
201
            The input Field, defined on the Operator's target domain
202
203
204

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
205
        Field
Philipp Arras's avatar
Philipp Arras committed
206
            The processed Field defined on the Operator's domain.
207
        """
Martin Reinecke's avatar
Martin Reinecke committed
208
209
210
        return self.apply(x, self.INVERSE_TIMES)

    def adjoint_times(self, x):
211
212
213
214
215
        """Applies the adjoint-Operator to a given Field.

        Parameters
        ----------
        x : Field
Philipp Arras's avatar
Philipp Arras committed
216
            The input Field, defined on the Operator's target domain
217
218
219

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
220
        Field
Philipp Arras's avatar
Philipp Arras committed
221
            The processed Field defined on the Operator's domain.
222
        """
Martin Reinecke's avatar
Martin Reinecke committed
223
        return self.apply(x, self.ADJOINT_TIMES)
224

Martin Reinecke's avatar
Martin Reinecke committed
225
    def adjoint_inverse_times(self, x):
226
        """Applies the adjoint-inverse Operator to a given Field.
227
228
229
230

        Parameters
        ----------
        x : Field
Philipp Arras's avatar
Philipp Arras committed
231
            The input Field, defined on the Operator's domain.
232
233
234

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
235
        Field
Philipp Arras's avatar
Philipp Arras committed
236
            The processed Field defined on the Operator's target domain.
237
238
239
240
241
242

        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
243
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
244

Martin Reinecke's avatar
Martin Reinecke committed
245
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
246
        """Same as :meth:`adjoint_inverse_times`"""
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
250
    def _check_mode(self, mode):
        if not self._validMode[mode]:
Martin Reinecke's avatar
Martin Reinecke committed
251
            raise NotImplementedError("invalid operator mode specified")
Martin Reinecke's avatar
Martin Reinecke committed
252
        if mode & self.capability == 0:
Martin Reinecke's avatar
Martin Reinecke committed
253
254
            raise NotImplementedError(
                "requested operator mode is not supported")
255

Martin Reinecke's avatar
Martin Reinecke committed
256
257
    def _check_input(self, x, mode):
        self._check_mode(mode)
Martin Reinecke's avatar
Martin Reinecke committed
258
        self._check_domain_equality(self._dom(mode), x.domain)