propagator_operator.py 5.76 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# NIFTy
# Copyright (C) 2017  Theo Steininger
#
# Author: Theo Steininger
#
# 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/>.
18

19
from nifty.operators import EndomorphicOperator,\
20
21
                            FFTOperator,\
                            InvertibleOperatorMixin
22
23


24
class PropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):
Theo Steininger's avatar
Theo Steininger committed
25
    """ NIFTY Propagator Operator D.
26

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
    The propagator operator D, is known from the Wiener Filter.
    Its inverse functional form might look like:
    D = (S^(-1) + M)^(-1)
    D = (S^(-1) + N^(-1))^(-1)
    D = (S^(-1) + R^(\dagger) N^(-1) R)^(-1)

    Parameters
    ----------
        S : LinearOperator
            Covariance of the signal prior.
        M : LinearOperator
            Likelihood contribution.
        R : LinearOperator
            Response operator translating signal to (noiseless) data.
        N : LinearOperator
            Covariance of the noise prior or the likelihood, respectively.
        inverter : class to invert explicitly defined operators
            (default:ConjugateGradient)
        preconditioner : Field
            numerical preconditioner to speed up convergence
47
48
49
        default_spaces : tuple of ints *optional*
            Defines on which space(s) of a given field the Operator acts by
            default (default: None)
50
51
52

    Attributes
    ----------
53
54
55
56
57
58
59
60
61
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain on which the Operator's input Field lives.
    target : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain in which the outcome of the operator lives. As the Operator
        is endomorphic this is the same as its domain.
    unitary : boolean
        Indicates whether the Operator is unitary or not.
    self_adjoint : boolean
        Indicates whether the operator is self_adjoint or not.
62
63
64

    Raises
    ------
65
66
67
    ValueError
        is raised if
            * neither N nor M is given
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

    Notes
    -----

    Examples
    --------
    >>> x_space = RGSpace(4)
    >>> k_space = RGRGTransformation.get_codomain(x_space)
    >>> f = Field(x_space, val=[2, 4, 6, 8])
    >>> S = create_power_operator(k_space, spec=1)
    >>> N = DiagonalOperaor(f.domain, diag=1)
    >>> D = PropagatorOperator(S=S, N=N) # D^{-1} = S^{-1} + N^{-1}
    >>> D(f).val
    <distributed_data_object>
    array([ 1.,  2.,  3.,  4.]

    See Also
    --------
    Scientific reference
    https://arxiv.org/abs/0806.3474
Theo Steininger's avatar
Theo Steininger committed
88

89
90
    """

91
92
    # ---Overwritten properties and methods---

93
    def __init__(self, S=None, M=None, R=None, N=None, inverter=None,
94
                 preconditioner=None, default_spaces=None):
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
        """
            Sets the standard operator properties and `codomain`, `_A1`, `_A2`,
            and `RN` if required.

            Parameters
            ----------
            S : operator
                Covariance of the signal prior.
            M : operator
                Likelihood contribution.
            R : operator
                Response operator translating signal to (noiseless) data.
            N : operator
                Covariance of the noise prior or the likelihood, respectively.

        """
111
112
113
114
115
116
117
118
119
120
121
        # infer domain, and target
        if M is not None:
            self._domain = M.domain
            self._likelihood_times = M.times

        elif N is None:
            raise ValueError("Either M or N must be given!")

        elif R is not None:
            self._domain = R.domain
            self._likelihood_times = \
122
                lambda z: R.adjoint_times(N.inverse_times(R.times(z)))
123
        else:
124
125
            self._domain = N.domain
            self._likelihood_times = lambda z: N.inverse_times(z)
126

127
128
        self._S = S
        self._fft_S = FFTOperator(self._domain, target=self._S.domain)
129
130
131
132

        if preconditioner is None:
            preconditioner = self._S_times

133
        super(PropagatorOperator, self).__init__(inverter=inverter,
134
135
                                                 preconditioner=preconditioner,
                                                 default_spaces=default_spaces)
136
137
138
139
140

    # ---Mandatory properties and methods---

    @property
    def domain(self):
141
        return self._domain
142
143

    @property
Martin Reinecke's avatar
Martin Reinecke committed
144
    def self_adjoint(self):
145
146
147
148
149
150
151
152
        return True

    @property
    def unitary(self):
        return False

    # ---Added properties and methods---

153
154
155
156
    def _S_times(self, x, spaces=None):
            transformed_x = self._fft_S(x, spaces=spaces)
            y = self._S(transformed_x, spaces=spaces)
            transformed_y = self._fft_S.inverse_times(y, spaces=spaces)
157
158
159
160
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

161
162
163
164
    def _S_inverse_times(self, x, spaces=None):
            transformed_x = self._fft_S(x, spaces=spaces)
            y = self._S.inverse_times(transformed_x, spaces=spaces)
            transformed_y = self._fft_S.inverse_times(y, spaces=spaces)
165
166
167
168
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

169
170
    def _inverse_times(self, x, spaces):
        pre_result = self._S_inverse_times(x, spaces)
171
172
173
        pre_result += self._likelihood_times(x)
        result = x.copy_empty()
        result.set_val(pre_result, copy=False)
174
        return result