linear_operator.py 8.79 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
Martin Reinecke's avatar
Martin Reinecke committed
20
from ..utilities import NiftyMetaBase
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..field import Field
22
import numpy as np
23
24


Martin Reinecke's avatar
Martin Reinecke committed
25
class LinearOperator(NiftyMetaBase()):
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

Martin Reinecke's avatar
Martin Reinecke committed
51
    _ilog = (-1, 0, 1, -1, 2, -1, -1, -1, 3)
Martin Reinecke's avatar
Martin Reinecke committed
52
    _validMode = (False, True, True, False, True, False, False, False, True)
Martin Reinecke's avatar
Martin Reinecke committed
53
54
55
56
57
58
59
60
    _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
61
    _addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15)
Martin Reinecke's avatar
Martin Reinecke committed
62
    _backwards = 6
Martin Reinecke's avatar
Martin Reinecke committed
63
    _all_ops = 15
Martin Reinecke's avatar
Martin Reinecke committed
64
65
66
67
68
    TIMES = 1
    ADJOINT_TIMES = 2
    INVERSE_TIMES = 4
    ADJOINT_INVERSE_TIMES = 8
    INVERSE_ADJOINT_TIMES = 8
Martin Reinecke's avatar
Martin Reinecke committed
69
70
    ADJOINT_BIT = 1
    INVERSE_BIT = 2
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
79
    def __init__(self):
        pass
80

81
    @abc.abstractproperty
82
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
83
84
85
        """DomainTuple : the operator's input domain

            The domain on which the Operator's input Field lives."""
86
        raise NotImplementedError
87
88
89

    @abc.abstractproperty
    def target(self):
Martin Reinecke's avatar
Martin Reinecke committed
90
91
92
        """DomainTuple : the operator's output domain

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

Martin Reinecke's avatar
Martin Reinecke committed
95
96
97
98
    def _flip_modes(self, mode):
        from .operator_adapter import OperatorAdapter
        return self if mode == 0 else OperatorAdapter(self, mode)

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

Martin Reinecke's avatar
Martin Reinecke committed
107
108
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
109
110
111
112
        """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
113
        return self._flip_modes(self.ADJOINT_BIT)
114

115
116
117
118
119
120
121
122
123
124
    @staticmethod
    def _toOperator(thing, dom):
        from .diagonal_operator import DiagonalOperator
        from .scaling_operator import ScalingOperator
        if isinstance(thing, LinearOperator):
            return thing
        if isinstance(thing, Field):
            return DiagonalOperator(thing)
        if np.isscalar(thing):
            return ScalingOperator(thing, dom)
Martin Reinecke's avatar
Martin Reinecke committed
125
        return NotImplemented
126

Martin Reinecke's avatar
Martin Reinecke committed
127
128
    def __mul__(self, other):
        from .chain_operator import ChainOperator
129
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
130
        return ChainOperator.make([self, other])
131

132
133
134
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
        other = self._toOperator(other, self.target)
Martin Reinecke's avatar
Martin Reinecke committed
135
        return ChainOperator.make([other, self])
136

Martin Reinecke's avatar
Martin Reinecke committed
137
138
    def __add__(self, other):
        from .sum_operator import SumOperator
139
        other = self._toOperator(other, self.domain)
140
        return SumOperator.make([self, other], [False, False])
141

142
143
144
    def __radd__(self, other):
        return self.__add__(other)

Martin Reinecke's avatar
Martin Reinecke committed
145
146
    def __sub__(self, other):
        from .sum_operator import SumOperator
147
        other = self._toOperator(other, self.domain)
148
        return SumOperator.make([self, other], [False, True])
149

150
151
152
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
153
        return SumOperator.make([other, self], [False, True])
154

Martin Reinecke's avatar
Martin Reinecke committed
155
156
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
157
        """int : the supported operation modes
158

159
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
160
161
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
162
        """
Martin Reinecke's avatar
Martin Reinecke committed
163
        raise NotImplementedError
164

Martin Reinecke's avatar
Martin Reinecke committed
165
166
    @abc.abstractmethod
    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
167
        """ Applies the Operator to a given `x`, in a specified `mode`.
168
169
170
171
172
173
174

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

Martin Reinecke's avatar
Martin Reinecke committed
175
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
176
177
178
179
180
            - :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
181
182
183

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
184
        Field
185
186
187
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
188
189
190
        raise NotImplementedError

    def __call__(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
191
        """Same as :meth:`times`"""
Martin Reinecke's avatar
Martin Reinecke committed
192
        return self.apply(x, self.TIMES)
193

Martin Reinecke's avatar
Martin Reinecke committed
194
    def times(self, x):
195
196
197
198
199
200
201
202
203
        """ 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
204
        Field
205
206
            The processed Field living on the Operator's target domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
207
208
209
        return self.apply(x, self.TIMES)

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

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

Martin Reinecke's avatar
Martin Reinecke committed
239
    def adjoint_inverse_times(self, x):
240
241
242
243
244
245
246
247
248
        """ 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
249
        Field
250
251
252
253
254
255
256
            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
257
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
258

Martin Reinecke's avatar
Martin Reinecke committed
259
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
260
        """Same as :meth:`adjoint_inverse_times`"""
Martin Reinecke's avatar
Martin Reinecke committed
261
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
262

Martin Reinecke's avatar
Martin Reinecke committed
263
264
265
266
267
    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")
268

Martin Reinecke's avatar
Martin Reinecke committed
269
    def _check_input(self, x, mode):
270
        if not isinstance(x, Field):
Martin Reinecke's avatar
updates    
Martin Reinecke committed
271
            raise ValueError("supplied object is not a `Field`.")
272

Martin Reinecke's avatar
Martin Reinecke committed
273
274
        self._check_mode(mode)
        if x.domain != self._dom(mode):
Martin Reinecke's avatar
Martin Reinecke committed
275
            raise ValueError("The operator's and field's domains don't match.")