log_normal_wiener_filter_curvature.py 2.55 KB
Newer Older
1
from ...operators import EndomorphicOperator,\
2
                            InvertibleOperatorMixin
3
from ...memoization import memo
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
4
from ...basic_arithmetics import exp
5
from ...sugar import create_composed_fft_operator
6

7
8
9

class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
                                     EndomorphicOperator):
10
11
12
13
14
15
16
17
18
19
    """The curvature of the LogNormalWienerFilterEnergy.

    This operator implements the second derivative of the
    LogNormalWienerFilterEnergy used in some minimization algorithms or for
    error estimates of the posterior maps. It is the inverse of the propagator
    operator.

    Parameters
    ----------
    R: LinearOperator,
Martin Reinecke's avatar
Martin Reinecke committed
20
21
22
       The response operator of the Wiener filter measurement.
    N: EndomorphicOperator
       The noise covariance.
23
    S: DiagonalOperator,
Martin Reinecke's avatar
Martin Reinecke committed
24
       The prior signal covariance
25
26
    """

Martin Reinecke's avatar
adjust  
Martin Reinecke committed
27
    def __init__(self, R, N, S, d, position, inverter, fft4exp=None, **kwargs):
28
29
30
        self.R = R
        self.N = N
        self.S = S
31
32
        self.d = d
        self.position = position
33
34
35
36
37
38
39

        if fft4exp is None:
            self._fft = create_composed_fft_operator(self.domain,
                                                     all_to='position')
        else:
            self._fft = fft4exp

40
        super(LogNormalWienerFilterCurvature, self).__init__(
41
42
43
44
45
                                                 inverter=inverter,
                                                 **kwargs)

    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
46
        return self.S.domain
47
48
49
50
51
52
53
54
55

    @property
    def self_adjoint(self):
        return True

    @property
    def unitary(self):
        return False

56
    def _times(self, x):
57
        part1 = self.S.inverse_times(x)
Theo Steininger's avatar
Theo Steininger committed
58
        # part2 = self._exppRNRexppd * x
59
60
61
62
63
        part3 = self._fft.adjoint_times(self._expp_sspace * self._fft(x))
        part3 = self._fft.adjoint_times(
                    self._expp_sspace *
                    self._fft(self.R.adjoint_times(
                                self.N.inverse_times(self.R(part3)))))
Theo Steininger's avatar
Theo Steininger committed
64
        return part1 + part3  # + part2
65

66
67
68
    @property
    @memo
    def _expp_sspace(self):
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
69
        return exp(self._fft(self.position))
70

71
72
73
    @property
    @memo
    def _Rexppd(self):
Martin Reinecke's avatar
Martin Reinecke committed
74
75
        expp = self._fft.adjoint_times(self._expp_sspace)
        return self.R(expp) - self.d
76
77
78
79
80
81
82
83
84

    @property
    @memo
    def _NRexppd(self):
        return self.N.inverse_times(self._Rexppd)

    @property
    @memo
    def _exppRNRexppd(self):
85
86
87
        return self._fft.adjoint_times(
                    self._expp_sspace *
                    self._fft(self.R.adjoint_times(self._NRexppd)))