propagator_operator.py 4 KB
Newer Older
1
2
# -*- coding: utf-8 -*-
from nifty.minimization import ConjugateGradient
3
from nifty.field import Field
4
5
from nifty.operators import EndomorphicOperator,\
                            FFTOperator
6
7


8
class PropagatorOperator(EndomorphicOperator):
9
10
11

    # ---Overwritten properties and methods---

12
13
    def __init__(self, S=None, M=None, R=None, N=None, inverter=None,
                 preconditioner=None):
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
        """
            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.

        """
30
31
32
33
34
35
36
37
38
39
40
        # 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 = \
41
                lambda z: R.adjoint_times(N.inverse_times(R.times(z)))
42
        else:
43
44
            self._domain = N.domain
            self._likelihood_times = lambda z: N.inverse_times(z)
45

46
47
        self._S = S
        self._fft_S = FFTOperator(self._domain, target=self._S.domain)
48
49
50
51
52

        if preconditioner is None:
            preconditioner = self._S_times

        self.preconditioner = preconditioner
53
54
55
56

        if inverter is not None:
            self.inverter = inverter
        else:
57
58
            self.inverter = ConjugateGradient(
                                preconditioner=self.preconditioner)
59
60
61
62
63

    # ---Mandatory properties and methods---

    @property
    def domain(self):
64
        return self._domain
65
66
67

    @property
    def field_type(self):
68
        return ()
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

    @property
    def implemented(self):
        return True

    @property
    def symmetric(self):
        return True

    @property
    def unitary(self):
        return False

    # ---Added properties and methods---

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
    def _S_times(self, x, spaces=None, types=None):
            transformed_x = self._fft_S(x,
                                        spaces=spaces,
                                        types=types)
            y = self._S(transformed_x, spaces=spaces, types=types)
            transformed_y = self._fft_S.inverse_times(y,
                                                      spaces=spaces,
                                                      types=types)
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

    def _S_inverse_times(self, x, spaces=None, types=None):
            transformed_x = self._fft_S(x,
                                        spaces=spaces,
                                        types=types)
            y = self._S.inverse_times(transformed_x,
                                      spaces=spaces,
                                      types=types)
            transformed_y = self._fft_S.inverse_times(y,
                                                      spaces=spaces,
                                                      types=types)
            result = x.copy_empty()
            result.set_val(transformed_y, copy=False)
            return result

    def _times(self, x, spaces, types, x0=None):
        if x0 is None:
112
            x0 = Field(self.target, val=0., dtype=x.dtype)
113

114
115
116
        (result, convergence) = self.inverter(A=self.inverse_times,
                                              b=x,
                                              x0=x0)
117
118
        return result

119
    def _inverse_times(self, x, spaces, types):
120
121
122
123
        pre_result = self._S_inverse_times(x, spaces, types)
        pre_result += self._likelihood_times(x)
        result = x.copy_empty()
        result.set_val(pre_result, copy=False)
124
        return result