propagator_operator.py 5.12 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
47
48
49
50
51
52
    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

    Attributes
    ----------

    Raises
    ------
53
54
55
    ValueError
        is raised if
            * neither N nor M is given
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

    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
76

77
78
    """

79
80
    # ---Overwritten properties and methods---

81
    def __init__(self, S=None, M=None, R=None, N=None, inverter=None,
82
                 preconditioner=None, default_spaces=None):
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
        """
            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.

        """
99
100
101
102
103
104
105
106
107
108
109
        # 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 = \
110
                lambda z: R.adjoint_times(N.inverse_times(R.times(z)))
111
        else:
112
113
            self._domain = N.domain
            self._likelihood_times = lambda z: N.inverse_times(z)
114

115
116
        self._S = S
        self._fft_S = FFTOperator(self._domain, target=self._S.domain)
117
118
119
120

        if preconditioner is None:
            preconditioner = self._S_times

121
        super(PropagatorOperator, self).__init__(inverter=inverter,
122
123
                                                 preconditioner=preconditioner,
                                                 default_spaces=default_spaces)
124
125
126
127
128

    # ---Mandatory properties and methods---

    @property
    def domain(self):
129
        return self._domain
130
131

    @property
Martin Reinecke's avatar
Martin Reinecke committed
132
    def self_adjoint(self):
133
134
135
136
137
138
139
140
        return True

    @property
    def unitary(self):
        return False

    # ---Added properties and methods---

141
142
143
144
    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)
145
146
147
148
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

149
150
151
152
    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)
153
154
155
156
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

157
158
    def _inverse_times(self, x, spaces):
        pre_result = self._S_inverse_times(x, spaces)
159
160
161
        pre_result += self._likelihood_times(x)
        result = x.copy_empty()
        result.set_val(pre_result, copy=False)
162
        return result