linear_operator.py 8.94 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
    """NIFTY base class for linear operators.

    The base NIFTY operator class is an abstract class from which
    other specific operator subclasses are derived.
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

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

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

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

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

73
74
    def __init__(self):
        pass
75

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

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

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

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

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

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

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

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

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

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

135
136
137
    def __radd__(self, other):
        return self.__add__(other)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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