log_normal_wiener_filter_energy.py 2.54 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
from ..minimization.energy import Energy
from ..utilities import memo
from .log_normal_wiener_filter_curvature import LogNormalWienerFilterCurvature
from ..sugar import create_composed_fft_operator
Martin Reinecke's avatar
Martin Reinecke committed
5
from ..field import exp
Martin Reinecke's avatar
Martin Reinecke committed
6

7
8
9
10
11
12
13
14
15
16

class LogNormalWienerFilterEnergy(Energy):
    """The Energy for the log-normal Wiener filter.

    It covers the case of linear measurement with
    Gaussian noise and Gaussain signal prior with known covariance.

    Parameters
    ----------
    position: Field,
Martin Reinecke's avatar
Martin Reinecke committed
17
18
19
20
21
22
23
24
25
       The current position.
    d: Field,
       the data.
    R: Operator,
       The response operator, describtion of the measurement process.
    N: EndomorphicOperator,
       The noise covariance in data space.
    S: EndomorphicOperator,
       The prior signal covariance in harmonic space.
26
27
    """

Martin Reinecke's avatar
adjust    
Martin Reinecke committed
28
    def __init__(self, position, d, R, N, S, inverter, fft4exp=None):
29
        super(LogNormalWienerFilterEnergy, self).__init__(position=position)
30
31
32
33
        self.d = d
        self.R = R
        self.N = N
        self.S = S
Martin Reinecke's avatar
adjust    
Martin Reinecke committed
34
        self._inverter = inverter
35

36
37
38
39
40
41
        if fft4exp is None:
            self._fft = create_composed_fft_operator(self.S.domain,
                                                     all_to='position')
        else:
            self._fft = fft4exp

42
43
    def at(self, position):
        return self.__class__(position=position, d=self.d, R=self.R, N=self.N,
Martin Reinecke's avatar
adjust    
Martin Reinecke committed
44
45
                              S=self.S, fft4exp=self._fft,
                              inverter=self._inverter)
46
47
48
49

    @property
    @memo
    def value(self):
Theo Steininger's avatar
Theo Steininger committed
50
        return 0.5*(self.position.vdot(self._Sp) +
Martin Reinecke's avatar
Martin Reinecke committed
51
                    self._Rexppd.vdot(self._NRexppd))
52
53
54
55

    @property
    @memo
    def gradient(self):
Martin Reinecke's avatar
Martin Reinecke committed
56
        return self._Sp + self._exppRNRexppd
57
58

    @property
Martin Reinecke's avatar
Martin Reinecke committed
59
    @memo
60
    def curvature(self):
61
        return LogNormalWienerFilterCurvature(
Martin Reinecke's avatar
Martin Reinecke committed
62
63
            R=self.R, N=self.N, S=self.S, fft=self._fft,
            expp_sspace = self._expp_sspace, inverter=self._inverter)
64
65
66

    @property
    @memo
67
    def _Sp(self):
68
        return self.S.inverse_times(self.position)
Martin Reinecke's avatar
Martin Reinecke committed
69

Martin Reinecke's avatar
Martin Reinecke committed
70
71
72
73
74
    @property
    @memo
    def _expp_sspace(self):
        return exp(self._fft(self.position))

Martin Reinecke's avatar
Martin Reinecke committed
75
76
77
    @property
    @memo
    def _Rexppd(self):
Martin Reinecke's avatar
Martin Reinecke committed
78
        expp = self._fft.adjoint_times(self._expp_sspace)
Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
82
83
84
85
86
87
88
89
        return self.R(expp) - self.d

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

    @property
    @memo
    def _exppRNRexppd(self):
        return self._fft.adjoint_times(
Martin Reinecke's avatar
Martin Reinecke committed
90
                    self._expp_sspace *
Martin Reinecke's avatar
Martin Reinecke committed
91
                    self._fft(self.R.adjoint_times(self._NRexppd)))