propagator_operator.py 3.65 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):
25
26
27

    # ---Overwritten properties and methods---

28
29
    def __init__(self, S=None, M=None, R=None, N=None, inverter=None,
                 preconditioner=None):
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
        """
            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.

        """
46
47
48
49
50
51
52
53
54
55
56
        # 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 = \
57
                lambda z: R.adjoint_times(N.inverse_times(R.times(z)))
58
        else:
59
60
            self._domain = N.domain
            self._likelihood_times = lambda z: N.inverse_times(z)
61

62
63
        self._S = S
        self._fft_S = FFTOperator(self._domain, target=self._S.domain)
64
65
66
67

        if preconditioner is None:
            preconditioner = self._S_times

68
69
        super(PropagatorOperator, self).__init__(inverter=inverter,
                                                 preconditioner=preconditioner)
70
71
72
73
74

    # ---Mandatory properties and methods---

    @property
    def domain(self):
75
        return self._domain
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

    @property
    def implemented(self):
        return True

    @property
    def symmetric(self):
        return True

    @property
    def unitary(self):
        return False

    # ---Added properties and methods---

91
92
93
94
    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)
95
96
97
98
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

99
100
101
102
    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)
103
104
105
106
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

107
108
    def _inverse_times(self, x, spaces):
        pre_result = self._S_inverse_times(x, spaces)
109
110
111
        pre_result += self._likelihood_times(x)
        result = x.copy_empty()
        result.set_val(pre_result, copy=False)
112
        return result