critical_power_energy.py 4.78 KB
Newer Older
1
from nifty.energies.energy import Energy
2
3
from nifty.operators.smoothness_operator import SmoothnessOperator
from nifty.library.critical_filter import CriticalPowerCurvature
4
from nifty.energies.memoization import memo
Jakob Knollmueller's avatar
Jakob Knollmueller committed
5

6
from nifty.sugar import generate_posterior_sample
7
from nifty import Field, exp
8

9

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
13
14
    It describes the energy of the logarithmic amplitudes of the power spectrum of
    a Gaussian random field under reconstruction uncertainty with smoothness and
Jakob Knollmueller's avatar
Jakob Knollmueller committed
15
    inverse gamma prior. It is used to infer the correlation structure of a correlated signal.
16
17
18

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
52
    def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0.,
53
                 logarithmic = True, samples=3, w=None, inverter=None):
54
        super(CriticalPowerEnergy, self).__init__(position = position)
55
        self.m = m
56
57
58
        self.D = D
        self.samples = samples
        self.sigma = sigma
59
60
        self.alpha = Field(self.position.domain, val=alpha)
        self.q = Field(self.position.domain, val=q)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
61
62
        self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma,
                                    logarithmic=logarithmic)
63
        self.rho = self.position.domain[0].rho
64
        self.inverter = inverter
65
66
        self.w = w
        if self.w is None:
67
            self.w = self._calculate_w(self.m, self.D, self.samples)
68
        self.theta = (exp(-self.position) * (self.q + self.w / 2.))
69

70

71
    def at(self, position):
72
73
74
75
76
        return self.__class__(position, self.m, D=self.D,
                              alpha =self.alpha,
                              q=self.q,
                              sigma=self.sigma, w=self.w,
                              samples=self.samples)
77
78
79

    @property
    def value(self):
80
81
82
        energy = self._theta.vdot(Field(self.position.domain,val=1.), bare= True)
        energy += self.position.vdot(self._rho_prime, bare=True)
        energy += 0.5 * self.position.vdot(self._Tt)
83
84
85
86
        return energy.real

    @property
    def gradient(self):
87
88
89
        gradient = - self._theta.weight(-1)
        gradient += (self._rho_prime).weight(-1)
        gradient +=  self._Tt
90
        gradient.val = gradient.val.real
91
92
93
94
        return gradient

    @property
    def curvature(self):
95
96
        curvature = CriticalPowerCurvature(theta=self._theta.weight(-1), T=self.T,
                                           inverter=self.inverter)
97
98
99
        return curvature

    def _calculate_w(self, m, D, samples):
100
        w = Field(domain=self.position.domain, val=0. , dtype=m.dtype)
101
102
        if D is not None:
            for i in range(samples):
103
                posterior_sample = generate_posterior_sample(m, D, inverter = self.inverter)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
104
                projected_sample = posterior_sample.power_analyze(
105
                    logarithmic=self.position.domain[0].config["logarithmic"],
Jakob Knollmueller's avatar
Jakob Knollmueller committed
106
107
                    nbin= self.position.domain[0].config["nbin"])
                w += (projected_sample) * self.rho
108
109
            w /= float(samples)
        else:
Jakob Knollmueller's avatar
Jakob Knollmueller committed
110
            w = m.power_analyze(
111
                    logarithmic=self.position.domain[0].config["logarithmic"],
Jakob Knollmueller's avatar
Jakob Knollmueller committed
112
                        nbin=self.position.domain[0].config["nbin"])
Jakob Knollmueller's avatar
Jakob Knollmueller committed
113
            w *= self.rho
114

115
        return w
116
117


118
119
120
121
122
123
124
125
126
127
128
129
130
131
    @property
    @memo
    def _theta(self):
        return (exp(-self.position) * (self.q + self.w / 2.))

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

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