critical_power_energy.py 5.61 KB
Newer Older
1
2
from ...energies.energy import Energy
from ...operators.smoothness_operator import SmoothnessOperator
3
from ...operators.diagonal_operator import DiagonalOperator
4
from ...operators.power_projection_operator import PowerProjection
5
from . import CriticalPowerCurvature
6
from ...memoization import memo
7
from ...minimization import ConjugateGradient
8

9
10
from ...sugar import generate_posterior_sample
from ... import Field, exp
11

12

13
class CriticalPowerEnergy(Energy):
Jakob Knollmueller's avatar
Jakob Knollmueller committed
14
    """The Energy of the power spectrum according to the critical filter.
15

16
17
18
19
    It describes the energy of the logarithmic amplitudes of the power spectrum
    of a Gaussian random field under reconstruction uncertainty with smoothness
    and inverse gamma prior. It is used to infer the correlation structure of a
    correlated signal.
20
21
22

    Parameters
    ----------
Jakob Knollmueller's avatar
Jakob Knollmueller committed
23
24
25
26
27
28
29
30
31
    position : Field,
        The current position of this energy.
    m : Field,
        The map whichs power spectrum has to be inferred
    D : EndomorphicOperator,
        The curvature of the Gaussian encoding the posterior covariance.
        If not specified, the map is assumed to be no reconstruction.
        default : None
    alpha : float
32
33
        The spectral prior of the inverse gamma distribution. 1.0 corresponds
        to non-informative.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
34
35
        default : 1.0
    q : float
36
37
38
39
40
        The cutoff parameter of the inverse gamma distribution. 0.0 corresponds
        to non-informative.
        default : 0.0
    smoothness_prior : float
        Controls the strength of the smoothness prior
Jakob Knollmueller's avatar
Jakob Knollmueller committed
41
42
43
44
        default : 0.0
    logarithmic : boolean
        Whether smoothness acts on linear or logarithmic scale.
    samples : integer
45
46
        Number of samples used for the estimation of the uncertainty
        corrections.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
47
48
49
50
51
        default : 3
    w : Field
        The contribution from the map with or without uncertainty. It is used
        to pass on the result of the costly sampling during the minimization.
        default : None
52
53
54
    inverter : ConjugateGradient
        The inversion strategy to invert the curvature and to generate samples.
        default : None
55
56
    """

57
58
    # ---Overwritten properties and methods---

59
    def __init__(self, position, m, D=None, alpha=1.0, q=0.,
60
                 smoothness_prior=0., logarithmic=True, samples=3, w=None,
61
62
63
64
                 inverter=None, gradient=None, curvature=None):
        super(CriticalPowerEnergy, self).__init__(position=position,
                                                  gradient=gradient,
                                                  curvature=curvature)
65
        self.m = m
66
67
        self.D = D
        self.samples = samples
68
69
        self.alpha = Field(self.position.domain, val=alpha)
        self.q = Field(self.position.domain, val=q)
70
        self.T = SmoothnessOperator(domain=self.position.domain[0],
71
                                    strength=smoothness_prior,
Jakob Knollmueller's avatar
Jakob Knollmueller committed
72
                                    logarithmic=logarithmic)
73
        self.rho = self.position.domain[0].rho
Martin Reinecke's avatar
Martin Reinecke committed
74
75
        self.P = PowerProjection(domain=self.m.domain,
                                 target=self.position.domain)
76
        self._w = w if w is not None else None
77
78
        if inverter is None:
            preconditioner = DiagonalOperator(self._theta.domain,
79
                                              diagonal=self._theta,
80
81
82
                                              copy=False)
            inverter = ConjugateGradient(preconditioner=preconditioner)
        self._inverter = inverter
Martin Reinecke's avatar
Martin Reinecke committed
83
        self.one = Field(self.position.domain, val=1.)
84
        self.constants = self.one/2. + self.alpha - 1
85
86
87
88

    @property
    def inverter(self):
        return self._inverter
89

90
91
    # ---Mandatory properties and methods---

92
    def at(self, position):
93
94
        return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
                              q=self.q, smoothness_prior=self.smoothness_prior,
95
                              logarithmic=self.logarithmic,
96
                              w=self.w, samples=self.samples,
97
                              inverter=self.inverter)
98
99

    @property
100
    @memo
101
    def value(self):
102
        energy = self.one.vdot(self._theta)
103
        energy += self.position.vdot(self.constants)
104
        energy += 0.5 * self.position.vdot(self._Tt)
105
106
107
        return energy.real

    @property
108
    @memo
109
    def gradient(self):
110
        gradient = -self._theta
111
        gradient += (self.constants)
112
        gradient += self._Tt
113
        gradient.val = gradient.val.real
114
115
116
        return gradient

    @property
117
    @memo
118
    def curvature(self):
119
        return CriticalPowerCurvature(theta=self._theta, T=self.T,
120
                                      inverter=self.inverter)
121

122
123
124
125
126
127
128
129
130
131
    # ---Added properties and methods---

    @property
    def logarithmic(self):
        return self.T.logarithmic

    @property
    def smoothness_prior(self):
        return self.T.strength

132
133
134
    @property
    def w(self):
        if self._w is None:
Theo Steininger's avatar
Theo Steininger committed
135
            self.logger.info("Initializing w")
136
137
138
            w = Field(domain=self.position.domain, val=0., dtype=self.m.dtype)
            if self.D is not None:
                for i in range(self.samples):
Theo Steininger's avatar
Theo Steininger committed
139
                    self.logger.info("Drawing sample %i" % i)
140
141
                    posterior_sample = generate_posterior_sample(
                                                            self.m, self.D)
142
143
                    w += self.P(abs(posterior_sample) ** 2)

144
145
                w /= float(self.samples)
            else:
146
147
                w = self.P(abs(self.m)**2)
            self._w = w
148
        return self._w
149

150
151
152
    @property
    @memo
    def _theta(self):
153
        return exp(-self.position) * (self.q + self.w / 2.)
154
155
156
157
158

    @property
    @memo
    def _Tt(self):
        return self.T(self.position)