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

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

11

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

15
16
17
18
    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.
19
20
21

    Parameters
    ----------
Jakob Knollmueller's avatar
Jakob Knollmueller committed
22
23
24
25
26
27
28
29
30
    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
31
32
        The spectral prior of the inverse gamma distribution. 1.0 corresponds
        to non-informative.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
33
34
        default : 1.0
    q : float
35
36
37
38
39
        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
40
41
42
43
        default : 0.0
    logarithmic : boolean
        Whether smoothness acts on linear or logarithmic scale.
    samples : integer
44
45
        Number of samples used for the estimation of the uncertainty
        corrections.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
46
47
48
49
50
        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
51
52
53
    inverter : ConjugateGradient
        The inversion strategy to invert the curvature and to generate samples.
        default : None
54
55
    """

56
57
    # ---Overwritten properties and methods---

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

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

85
86
    # ---Mandatory properties and methods---

87
    def at(self, position):
88
89
        return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
                              q=self.q, smoothness_prior=self.smoothness_prior,
90
                              logarithmic=self.logarithmic,
91
                              w=self.w, samples=self.samples,
92
                              inverter=self.inverter)
93
94

    @property
95
    @memo
96
    def value(self):
97
        energy = self._theta.sum()
Martin Reinecke's avatar
Martin Reinecke committed
98
        energy += self.position.weight(-1).vdot(self._rho_prime)
99
        energy += 0.5 * self.position.vdot(self._Tt)
100
101
102
        return energy.real

    @property
103
    @memo
104
    def gradient(self):
105
        gradient = -self._theta.weight(-1)
106
        gradient += (self._rho_prime).weight(-1)
107
        gradient += self._Tt
108
        gradient.val = gradient.val.real
109
110
111
        return gradient

    @property
112
    @memo
113
    def curvature(self):
114
115
        return CriticalPowerCurvature(theta=self._theta.weight(-1), T=self.T,
                                      inverter=self.inverter)
116

117
118
119
120
121
122
123
124
125
126
    # ---Added properties and methods---

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

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

127
128
129
    @property
    def w(self):
        if self._w is None:
Theo Steininger's avatar
Theo Steininger committed
130
            self.logger.info("Initializing w")
131
132
133
            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
134
                    self.logger.info("Drawing sample %i" % i)
135
136
137
                    posterior_sample = generate_posterior_sample(
                                                            self.m, self.D)
                    projected_sample = posterior_sample.power_analyze(
Martin Reinecke's avatar
Martin Reinecke committed
138
                     binbounds=self.position.domain[0].binbounds)
139
140
141
142
                    w += (projected_sample) * self.rho
                w /= float(self.samples)
            else:
                w = self.m.power_analyze(
Martin Reinecke's avatar
Martin Reinecke committed
143
                     binbounds=self.position.domain[0].binbounds)
144
145
146
                w *= self.rho
            self._w = w
        return self._w
147

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

    @property
    @memo
    def _rho_prime(self):
        return self.alpha - 1. + self.rho / 2.

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