log_normal_wiener_filter_energy.py 2.1 KB
Newer Older
1
from ...energies.energy import Energy
2
from ...memoization import memo
3
4
from . import LogNormalWienerFilterCurvature
from ...sugar import create_composed_fft_operator
5

6
7
8
9
10
11
12
13
14
15

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
16
17
18
19
20
21
22
23
24
       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.
25
26
    """

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

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

41
42
    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
43
44
                              S=self.S, fft4exp=self._fft,
                              inverter=self._inverter)
45
46
47
48

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

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

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

    @property
    @memo
67
    def _Sp(self):
68
        return self.S.inverse_times(self.position)