wiener_filter_curvature.py 2.22 KB
Newer Older
1
2
3
from ..operators import EndomorphicOperator, InversionEnabler
from ..field import Field, sqrt
from ..sugar import power_analyze, power_synthesize
4

Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
5

Martin Reinecke's avatar
Martin Reinecke committed
6
class WienerFilterCurvature(EndomorphicOperator):
Jakob Knollmueller's avatar
Jakob Knollmueller committed
7
8
9
10
    """The curvature of the WienerFilterEnergy.

    This operator implements the second derivative of the
    WienerFilterEnergy used in some minimization algorithms or
11
12
    for error estimates of the posterior maps. It is the
    inverse of the propagator operator.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
13
14
15
16

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

24
    def __init__(self, R, N, S, inverter):
Martin Reinecke's avatar
Martin Reinecke committed
25
        super(WienerFilterCurvature, self).__init__()
26
27
28
        self.R = R
        self.N = N
        self.S = S
Martin Reinecke's avatar
Martin Reinecke committed
29
30
        op = R.adjoint*N.inverse*R + S.inverse
        self._op = InversionEnabler(op, inverter, S.times)
31

32
33
    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
34
        return self.S.domain
35
36

    @property
Martin Reinecke's avatar
Martin Reinecke committed
37
38
    def capability(self):
        return self._op.capability
39

Martin Reinecke's avatar
Martin Reinecke committed
40
41
    def apply(self, x, mode):
        return self._op.apply(x, mode)
42

Martin Reinecke's avatar
Martin Reinecke committed
43
    def generate_posterior_sample(self):
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
        """ Generates a posterior sample from a Gaussian distribution with
        given mean and covariance.

        This method generates samples by setting up the observation and
        reconstruction of a mock signal in order to obtain residuals of the
        right correlation which are added to the given mean.

        Returns
        -------
        sample : Field
            Returns the a sample from the Gaussian of given mean and
            covariance.
        """

        power = sqrt(power_analyze(self.S.diagonal()))
        mock_signal = power_synthesize(power, real_signal=True)

        noise = self.N.diagonal().weight(-1)

        mock_noise = Field.from_random(random_type="normal",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
64
                                       domain=self.N.domain, dtype=noise.dtype)
65
66
67
68
69
70
        mock_noise *= sqrt(noise)

        mock_data = self.R(mock_signal) + mock_noise

        mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data))
        mock_m = self.inverse_times(mock_j)
Martin Reinecke's avatar
Martin Reinecke committed
71
        sample = mock_signal - mock_m
72
        return sample