linear_operator.py 8.31 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):
Martin Reinecke's avatar
Martin Reinecke committed
59 60 61
        """DomainTuple : the operator's input domain

            The domain on which the Operator's input Field lives."""
62
        raise NotImplementedError
63 64 65

    @abc.abstractproperty
    def target(self):
Martin Reinecke's avatar
Martin Reinecke committed
66 67 68
        """DomainTuple : the operator's output domain

            The domain on which the Operator's output Field lives."""
69 70
        raise NotImplementedError

Martin Reinecke's avatar
Martin Reinecke committed
71 72
    @property
    def inverse(self):
Martin Reinecke's avatar
Martin Reinecke committed
73 74 75 76
        """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
77 78
        from .inverse_operator import InverseOperator
        return InverseOperator(self)
79

Martin Reinecke's avatar
Martin Reinecke committed
80 81
    @property
    def adjoint(self):
Martin Reinecke's avatar
Martin Reinecke committed
82 83 84 85
        """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
86 87
        from .adjoint_operator import AdjointOperator
        return AdjointOperator(self)
88

89 90 91 92 93 94 95 96 97 98
    @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
99
        return NotImplemented
100

Martin Reinecke's avatar
Martin Reinecke committed
101 102
    def __mul__(self, other):
        from .chain_operator import ChainOperator
103
        other = self._toOperator(other, self.domain)
Martin Reinecke's avatar
Martin Reinecke committed
104
        return ChainOperator.make([self, other])
105

106 107 108
    def __rmul__(self, other):
        from .chain_operator import ChainOperator
        other = self._toOperator(other, self.target)
Martin Reinecke's avatar
Martin Reinecke committed
109
        return ChainOperator.make([other, self])
110

Martin Reinecke's avatar
Martin Reinecke committed
111 112
    def __add__(self, other):
        from .sum_operator import SumOperator
113
        other = self._toOperator(other, self.domain)
114
        return SumOperator.make([self, other], [False, False])
115

116 117 118
    def __radd__(self, other):
        return self.__add__(other)

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

Martin Reinecke's avatar
Martin Reinecke committed
124
    # MR FIXME: this might be more complicated ...
125 126 127
    def __rsub__(self, other):
        from .sum_operator import SumOperator
        other = self._toOperator(other, self.domain)
128
        return SumOperator.make(other, self, [False, True])
129

Martin Reinecke's avatar
Martin Reinecke committed
130 131
    @abc.abstractproperty
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
132
        """int : the supported operation modes
133

Martin Reinecke's avatar
Martin Reinecke committed
134 135 136
        Returns the suppoerted subset of :attr:`TIMES`, :attr:`ADJOINT_TIMES`,
        :attr:`INVERSE_TIMES`, and :attr:`ADJOINT_INVERSE_TIMES`,
        joined together by the "|" operator.
137
        """
Martin Reinecke's avatar
Martin Reinecke committed
138
        raise NotImplementedError
139

Martin Reinecke's avatar
Martin Reinecke committed
140 141
    @abc.abstractmethod
    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
142
        """ Applies the Operator to a given `x`, in a specified `mode`.
143 144 145 146 147 148 149

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

Martin Reinecke's avatar
Martin Reinecke committed
150
        mode : int
151 152 153 154 155 156 157 158 159
            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
160
        Field
161 162 163
            The processed Field living on the Operator's target or domain,
            depending on mode.
        """
Martin Reinecke's avatar
Martin Reinecke committed
164 165 166
        raise NotImplementedError

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

Martin Reinecke's avatar
Martin Reinecke committed
170
    def times(self, x):
171 172 173 174 175 176 177 178 179
        """ 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
180
        Field
181 182
            The processed Field living on the Operator's target domain.
        """
Martin Reinecke's avatar
Martin Reinecke committed
183 184 185
        return self.apply(x, self.TIMES)

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

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

Martin Reinecke's avatar
Martin Reinecke committed
215
    def adjoint_inverse_times(self, x):
216 217 218 219 220 221 222 223 224
        """ 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
225
        Field
226 227 228 229 230 231 232
            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
233
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
234

Martin Reinecke's avatar
Martin Reinecke committed
235
    def inverse_adjoint_times(self, x):
Martin Reinecke's avatar
Martin Reinecke committed
236
        """Same as :meth:`adjoint_inverse_times`"""
Martin Reinecke's avatar
Martin Reinecke committed
237
        return self.apply(x, self.ADJOINT_INVERSE_TIMES)
238

Martin Reinecke's avatar
Martin Reinecke committed
239 240 241 242 243
    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")
244

Martin Reinecke's avatar
Martin Reinecke committed
245
    def _check_input(self, x, mode):
246
        if not isinstance(x, Field):
Martin Reinecke's avatar
updates  
Martin Reinecke committed
247
            raise ValueError("supplied object is not a `Field`.")
248

Martin Reinecke's avatar
Martin Reinecke committed
249 250
        self._check_mode(mode)
        if x.domain != self._dom(mode):
Martin Reinecke's avatar
Martin Reinecke committed
251
            raise ValueError("The operator's and field's domains don't match.")
Martin Reinecke's avatar
Martin Reinecke committed
252 253

    def draw_sample(self):
Martin Reinecke's avatar
Martin Reinecke committed
254 255 256
        """Generate a zero-mean sample

        Generates a sample from a Gaussian distribution with zero mean and
Martin Reinecke's avatar
Martin Reinecke committed
257 258 259 260
        covariance given by the operator.

        Returns
        -------
Martin Reinecke's avatar
Martin Reinecke committed
261 262
        Field
            A sample from the Gaussian of given covariance.
Martin Reinecke's avatar
Martin Reinecke committed
263 264
        """
        raise NotImplementedError