linear_operator.py 8.55 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 NiftyMeta
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..field import Field
Martin Reinecke's avatar
Martin Reinecke committed
22
from future.utils import with_metaclass
23
import numpy as np
24
25


26
class LinearOperator(with_metaclass(
Martin Reinecke's avatar
Martin Reinecke committed
27
        NiftyMeta, type('NewBase', (object,), {}))):
28
29
30
31
32
    """NIFTY base class for linear operators.

    The base NIFTY operator class is an abstract class from which
    other specific operator subclasses are derived.
    """
Theo Steininger's avatar
Theo Steininger committed
33

Martin Reinecke's avatar
Martin Reinecke committed
34
35
36
37
38
39
    _validMode = (False, True, True, False, True, False, False, False, True)
    _inverseMode = (0, 4, 8, 0, 1, 0, 0, 0, 2)
    _inverseCapability = (0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15)
    _adjointMode = (0, 2, 1, 0, 8, 0, 0, 0, 4)
    _adjointCapability = (0, 2, 1, 3, 8, 10, 9, 11, 4, 6, 5, 7, 12, 14, 13, 15)
    _addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15)
Martin Reinecke's avatar
Martin Reinecke committed
40
    _backwards = 6
Martin Reinecke's avatar
Martin Reinecke committed
41
    _all_ops = 15
Martin Reinecke's avatar
Martin Reinecke committed
42
43
44
45
46
    TIMES = 1
    ADJOINT_TIMES = 2
    INVERSE_TIMES = 4
    ADJOINT_INVERSE_TIMES = 8
    INVERSE_ADJOINT_TIMES = 8
47

Martin Reinecke's avatar
Martin Reinecke committed
48
49
    def _dom(self, mode):
        return self.domain if (mode & 9) else self.target
50

Martin Reinecke's avatar
Martin Reinecke committed
51
52
    def _tgt(self, mode):
        return self.domain if (mode & 6) else self.target
53

54
55
    def __init__(self):
        pass
56

57
    @abc.abstractproperty
58
    def domain(self):
59
        """
Martin Reinecke's avatar
Martin Reinecke committed
60
        DomainTuple
Theo Steininger's avatar
Theo Steininger committed
61
62
            The domain on which the Operator's input Field lives.
            Every Operator which inherits from the abstract LinearOperator
63
64
            base class must have this attribute.
        """
65
        raise NotImplementedError
66
67
68

    @abc.abstractproperty
    def target(self):
69
        """
Martin Reinecke's avatar
Martin Reinecke committed
70
        DomainTuple
Theo Steininger's avatar
Theo Steininger committed
71
72
            The domain on which the Operator's output Field lives.
            Every Operator which inherits from the abstract LinearOperator
73
74
            base class must have this attribute.
        """
75
76
        raise NotImplementedError

Martin Reinecke's avatar
Martin Reinecke committed
77
78
    @property
    def inverse(self):
79
        """
Martin Reinecke's avatar
Martin Reinecke committed
80
        LinearOperator
81
82
83
            Returns a LinearOperator object which behaves as if it were
            the inverse of this operator.
        """
Martin Reinecke's avatar
Martin Reinecke committed
84
85
        from .inverse_operator import InverseOperator
        return InverseOperator(self)
86

Martin Reinecke's avatar
Martin Reinecke committed
87
88
    @property
    def adjoint(self):
89
        """
Martin Reinecke's avatar
Martin Reinecke committed
90
        LinearOperator
91
92
93
            Returns a LinearOperator object which behaves as if it were
            the adjoint of this operator.
        """
Martin Reinecke's avatar
Martin Reinecke committed
94
95
        from .adjoint_operator import AdjointOperator
        return AdjointOperator(self)
96

97
98
99
100
101
102
103
104
105
106
    @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
107
        return NotImplemented
108

Martin Reinecke's avatar
Martin Reinecke committed
109
110
    def __mul__(self, other):
        from .chain_operator import ChainOperator
111
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
112
        return ChainOperator.make([self, other])
113

114
115
116
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
        other = self._toOperator(other, self.target)
Martin Reinecke's avatar
Martin Reinecke committed
117
        return ChainOperator.make([other, self])
118

Martin Reinecke's avatar
Martin Reinecke committed
119
120
    def __add__(self, other):
        from .sum_operator import SumOperator
121
        other = self._toOperator(other, self.domain)
122
        return SumOperator.make([self, other], [False, False])
123

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

Martin Reinecke's avatar
Martin Reinecke committed
127
128
    def __sub__(self, other):
        from .sum_operator import SumOperator
129
        other = self._toOperator(other, self.domain)
130
        return SumOperator.make([self, other], [False, True])
131

Martin Reinecke's avatar
Martin Reinecke committed
132
    # MR FIXME: this might be more complicated ...
133
134
135
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
136
        return SumOperator.make(other, self, [False, True])
137

Martin Reinecke's avatar
Martin Reinecke committed
138
139
    @abc.abstractproperty
    def capability(self):
140
141
142
143
        """ Specifies the application modes supported by this operator

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
144
        int
145
146
147
148
            This is any subset of LinearOperator.{TIMES, ADJOINT_TIMES,
            INVERSE_TIMES, ADJOINT_INVERSE_TIMES, INVERSE_ADJOINT_TIMES},
            joined together by the "|" operator.
        """
Martin Reinecke's avatar
Martin Reinecke committed
149
        raise NotImplementedError
150

Martin Reinecke's avatar
Martin Reinecke committed
151
152
    @abc.abstractmethod
    def apply(self, x, mode):
153
154
155
156
157
158
159
160
        """ Applies the Operator to a given Field, in a specified mode.

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

Martin Reinecke's avatar
Martin Reinecke committed
161
        mode : int
162
163
164
165
166
167
168
169
170
            LinearOperator.TIMES: normal application
            LinearOperator.ADJOINT_TIMES: adjoint application
            LinearOperator.INVERSE_TIMES: inverse application
            LinearOperator.ADJOINT_INVERSE_TIMES or
            LinearOperator.INVERSE_ADJOINT_TIMES:
            adjoint inverse application

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
171
        Field
172
173
174
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
175
176
177
        raise NotImplementedError

    def __call__(self, x):
178
        """Same as times()"""
Martin Reinecke's avatar
Martin Reinecke committed
179
        return self.apply(x, self.TIMES)
180

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
246
    def inverse_adjoint_times(self, x):
247
        """Same as adjoint_inverse_times()"""
Martin Reinecke's avatar
Martin Reinecke committed
248
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
249

Martin Reinecke's avatar
Martin Reinecke committed
250
251
252
253
254
    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")
255

Martin Reinecke's avatar
Martin Reinecke committed
256
    def _check_input(self, x, mode):
257
        if not isinstance(x, Field):
Martin Reinecke's avatar
updates    
Martin Reinecke committed
258
            raise ValueError("supplied object is not a `Field`.")
259

Martin Reinecke's avatar
Martin Reinecke committed
260
261
        self._check_mode(mode)
        if x.domain != self._dom(mode):
Martin Reinecke's avatar
Martin Reinecke committed
262
            raise ValueError("The operator's and field's domains don't match.")
Martin Reinecke's avatar
Martin Reinecke committed
263
264

    def draw_sample(self):
Martin Reinecke's avatar
Martin Reinecke committed
265
        """ Generates a sample from a Gaussian distribution with zero mean and
Martin Reinecke's avatar
Martin Reinecke committed
266
267
268
269
        covariance given by the operator.

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
270
271
        Field
            A sample from the Gaussian of given covariance.
Martin Reinecke's avatar
Martin Reinecke committed
272
273
        """
        raise NotImplementedError