linear_operator.py 8.86 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
52
53
54
55
56
    _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
57
    _backwards = 6
Martin Reinecke's avatar
Martin Reinecke committed
58
    _all_ops = 15
Martin Reinecke's avatar
Martin Reinecke committed
59
60
61
62
63
    TIMES = 1
    ADJOINT_TIMES = 2
    INVERSE_TIMES = 4
    ADJOINT_INVERSE_TIMES = 8
    INVERSE_ADJOINT_TIMES = 8
64

Martin Reinecke's avatar
Martin Reinecke committed
65
66
    def _dom(self, mode):
        return self.domain if (mode & 9) else self.target
67

Martin Reinecke's avatar
Martin Reinecke committed
68
69
    def _tgt(self, mode):
        return self.domain if (mode & 6) else self.target
70

71
72
    def __init__(self):
        pass
73

74
    @abc.abstractproperty
75
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
76
77
78
        """DomainTuple : the operator's input domain

            The domain on which the Operator's input Field lives."""
79
        raise NotImplementedError
80
81
82

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

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

Martin Reinecke's avatar
Martin Reinecke committed
88
89
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
90
91
92
93
        """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
94
95
        from .inverse_operator import InverseOperator
        return InverseOperator(self)
96

Martin Reinecke's avatar
Martin Reinecke committed
97
98
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
99
100
101
102
        """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
103
104
        from .adjoint_operator import AdjointOperator
        return AdjointOperator(self)
105

106
107
108
109
110
111
112
113
114
115
    @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
116
        return NotImplemented
117

Martin Reinecke's avatar
Martin Reinecke committed
118
119
    def __mul__(self, other):
        from .chain_operator import ChainOperator
120
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
121
        return ChainOperator.make([self, other])
122

123
124
125
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
        other = self._toOperator(other, self.target)
Martin Reinecke's avatar
Martin Reinecke committed
126
        return ChainOperator.make([other, self])
127

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

133
134
135
    def __radd__(self, other):
        return self.__add__(other)

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

Martin Reinecke's avatar
Martin Reinecke committed
141
    # MR FIXME: this might be more complicated ...
142
143
144
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
145
        return SumOperator.make(other, self, [False, True])
146

Martin Reinecke's avatar
Martin Reinecke committed
147
148
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
149
        """int : the supported operation modes
150

151
        Returns the supported subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
Martin Reinecke's avatar
Martin Reinecke committed
152
153
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
154
        """
Martin Reinecke's avatar
Martin Reinecke committed
155
        raise NotImplementedError
156

Martin Reinecke's avatar
Martin Reinecke committed
157
158
    @abc.abstractmethod
    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
159
        """ Applies the Operator to a given `x`, in a specified `mode`.
160
161
162
163
164
165
166

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

Martin Reinecke's avatar
Martin Reinecke committed
167
        mode : int
Martin Reinecke's avatar
Martin Reinecke committed
168
169
170
171
172
            - :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
173
174
175

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
176
        Field
177
178
179
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
180
181
182
        raise NotImplementedError

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

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

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
251
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
252
        """Same as :meth:`adjoint_inverse_times`"""
Martin Reinecke's avatar
Martin Reinecke committed
253
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
254

Martin Reinecke's avatar
Martin Reinecke committed
255
256
257
258
259
    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")
260

Martin Reinecke's avatar
Martin Reinecke committed
261
    def _check_input(self, x, mode):
262
        if not isinstance(x, Field):
Martin Reinecke's avatar
updates    
Martin Reinecke committed
263
            raise ValueError("supplied object is not a `Field`.")
264

Martin Reinecke's avatar
Martin Reinecke committed
265
266
        self._check_mode(mode)
        if x.domain != self._dom(mode):
Martin Reinecke's avatar
Martin Reinecke committed
267
            raise ValueError("The operator's and field's domains don't match.")
Martin Reinecke's avatar
Martin Reinecke committed
268
269

    def draw_sample(self):
Martin Reinecke's avatar
Martin Reinecke committed
270
271
272
        """Generate a zero-mean sample

        Generates a sample from a Gaussian distribution with zero mean and
Martin Reinecke's avatar
Martin Reinecke committed
273
274
275
276
        covariance given by the operator.

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
277
278
        Field
            A sample from the Gaussian of given covariance.
Martin Reinecke's avatar
Martin Reinecke committed
279
280
        """
        raise NotImplementedError