critical_power_energy.py 4.97 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
from ..minimization.energy import Energy
from ..operators.smoothness_operator import SmoothnessOperator
from ..operators.power_projection_operator import PowerProjectionOperator
from ..operators.inversion_enabler import InversionEnabler
from .critical_power_curvature import CriticalPowerCurvature
from ..utilities import memo
from .. import Field, exp
from ..sugar import generate_posterior_sample
9

10

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

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

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

55
56
    # ---Overwritten properties and methods---

57
    def __init__(self, position, m, D=None, alpha=1.0, q=0.,
58
59
                 smoothness_prior=0., logarithmic=True, samples=3, w=None,
                 inverter=None):
60
        super(CriticalPowerEnergy, self).__init__(position=position)
61
        self.m = m
62
63
        self.D = D
        self.samples = samples
64
65
        self.alpha = Field(self.position.domain, val=alpha)
        self.q = Field(self.position.domain, val=q)
66
        self.T = SmoothnessOperator(domain=self.position.domain[0],
67
                                    strength=smoothness_prior,
Jakob Knollmueller's avatar
Jakob Knollmueller committed
68
                                    logarithmic=logarithmic)
Martin Reinecke's avatar
Martin Reinecke committed
69
70
        self.P = PowerProjectionOperator(domain=self.m.domain,
                                         power_space=self.position.domain[0])
71
72
        self._w = w
        self._inverter = inverter
73

74
75
    # ---Mandatory properties and methods---

76
    def at(self, position):
77
78
        return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
                              q=self.q, smoothness_prior=self.smoothness_prior,
79
                              logarithmic=self.logarithmic,
80
81
                              w=self.w, samples=self.samples,
                              inverter=self._inverter)
82
83

    @property
Martin Reinecke's avatar
Martin Reinecke committed
84
    @memo
85
    def value(self):
Martin Reinecke's avatar
Martin Reinecke committed
86
87
        energy = Field.ones_like(self.position).vdot(self._theta)
        energy += self.position.vdot(self.alpha-0.5)
88
        energy += 0.5 * self.position.vdot(self._Tt)
89
90
91
        return energy.real

    @property
Martin Reinecke's avatar
Martin Reinecke committed
92
    @memo
93
    def gradient(self):
Martin Reinecke's avatar
Martin Reinecke committed
94
95
        gradient = -self._theta
        gradient += self.alpha-0.5
96
        gradient += self._Tt
Martin Reinecke's avatar
Martin Reinecke committed
97
        return gradient.real
98
99

    @property
Martin Reinecke's avatar
Martin Reinecke committed
100
    @memo
101
    def curvature(self):
Martin Reinecke's avatar
Martin Reinecke committed
102
103
104
        curv = CriticalPowerCurvature(theta=self._theta, T=self.T)
        return InversionEnabler(curv, inverter=self._inverter,
                                preconditioner=curv.preconditioner)
105

106
107
108
109
110
111
112
113
114
115
    # ---Added properties and methods---

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

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

116
117
118
    @property
    def w(self):
        if self._w is None:
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
119
            # self.logger.info("Initializing w")
120
121
122
            w = Field(domain=self.position.domain, val=0., dtype=self.m.dtype)
            if self.D is not None:
                for i in range(self.samples):
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
123
                    # self.logger.info("Drawing sample %i" % i)
124
125
                    posterior_sample = generate_posterior_sample(
                                                            self.m, self.D)
Martin Reinecke's avatar
Martin Reinecke committed
126
127
                    w += self.P(abs(posterior_sample) ** 2)

128
129
                w /= float(self.samples)
            else:
Martin Reinecke's avatar
Martin Reinecke committed
130
                w = self.P(abs(self.m)**2)
131
132
            self._w = w
        return self._w
133

134
135
136
    @property
    @memo
    def _theta(self):
137
        return exp(-self.position) * (self.q + self.w / 2.)
138
139
140
141
142

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