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

Philipp Arras's avatar
Philipp Arras committed
18
19
from ..field import Field
from ..multi_field import MultiField
Martin Reinecke's avatar
Martin Reinecke committed
20
from .operator import Operator
Philipp Arras's avatar
Fixups    
Philipp Arras committed
21

22

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

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

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

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

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

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

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

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

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

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

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

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

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

121
122
123
    def __radd__(self, other):
        return self.__add__(other)

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

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

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

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

144
145
146
    def force(self, x):
        return self.apply(x.extract(self.domain), self.TIMES)

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

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

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

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

    def __call__(self, x):
Martin Reinecke's avatar
more    
Martin Reinecke committed
172
173
        if isinstance(x, Operator):
            return self @ x
Martin Reinecke's avatar
more    
Martin Reinecke committed
174
        if x.jac is not None:
Martin Reinecke's avatar
more    
Martin Reinecke committed
175
            return x.new(self(x.fld), self).prepend_jac(x.jac)
Martin Reinecke's avatar
more    
Martin Reinecke committed
176
177
        if x.val is not None:
            return self.apply(x, self.TIMES)
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)