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

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
83
    def __init__(self):
        pass
84

85
    @abc.abstractproperty
86
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
87
88
89
        """DomainTuple : the operator's input domain

            The domain on which the Operator's input Field lives."""
90
        raise NotImplementedError
91
92
93

    @abc.abstractproperty
    def target(self):
Martin Reinecke's avatar
Martin Reinecke committed
94
95
96
        """DomainTuple : the operator's output domain

            The domain on which the Operator's output Field lives."""
97
98
        raise NotImplementedError

99
    def _flip_modes(self, op_transform):
Martin Reinecke's avatar
Martin Reinecke committed
100
        from .operator_adapter import OperatorAdapter
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
101
102
        return self if op_transform == 0 \
            else OperatorAdapter(self, op_transform)
Martin Reinecke's avatar
Martin Reinecke committed
103

Martin Reinecke's avatar
Martin Reinecke committed
104
105
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
106
107
108
109
        """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
110
        return self._flip_modes(self.INVERSE_BIT)
111

Martin Reinecke's avatar
Martin Reinecke committed
112
113
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
114
115
116
117
        """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
118
        return self._flip_modes(self.ADJOINT_BIT)
119

120
121
122
123
124
125
126
127
128
129
    @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
130
        return NotImplemented
131

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

137
138
139
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
        other = self._toOperator(other, self.target)
Martin Reinecke's avatar
Martin Reinecke committed
140
        return ChainOperator.make([other, self])
141

Martin Reinecke's avatar
Martin Reinecke committed
142
143
    def __add__(self, other):
        from .sum_operator import SumOperator
144
        other = self._toOperator(other, self.domain)
145
        return SumOperator.make([self, other], [False, False])
146

147
148
149
    def __radd__(self, other):
        return self.__add__(other)

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

155
156
157
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
158
        return SumOperator.make([other, self], [False, True])
159

Martin Reinecke's avatar
Martin Reinecke committed
160
161
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
162
        """int : the supported operation modes
163

164
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
165
166
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
167
        """
Martin Reinecke's avatar
Martin Reinecke committed
168
        raise NotImplementedError
169

Martin Reinecke's avatar
Martin Reinecke committed
170
171
    @abc.abstractmethod
    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
172
        """ Applies the Operator to a given `x`, in a specified `mode`.
173
174
175
176
177
178
179

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

Martin Reinecke's avatar
Martin Reinecke committed
180
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
181
182
183
184
185
            - :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
186
187
188

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
189
        Field
190
191
192
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
193
194
195
        raise NotImplementedError

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
264
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
265
        """Same as :meth:`adjoint_inverse_times`"""
Martin Reinecke's avatar
Martin Reinecke committed
266
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
267

Martin Reinecke's avatar
Martin Reinecke committed
268
269
270
271
272
    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")
273

Martin Reinecke's avatar
Martin Reinecke committed
274
    def _check_input(self, x, mode):
275
        if not isinstance(x, Field):
Martin Reinecke's avatar
updates    
Martin Reinecke committed
276
            raise ValueError("supplied object is not a `Field`.")
277

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