critical_power_energy.py 4.23 KB
Newer Older
1
from nifty.energies.energy import Energy
2
from nifty.library.operator_library import CriticalPowerCurvature,\
3
                                            SmoothnessOperator
Jakob Knollmueller's avatar
Jakob Knollmueller committed
4

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

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
11
12
    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
13
    inverse gamma prior. It is used to infer the correlation structure of a correlated signal.
14
15
16

    Parameters
    ----------
Jakob Knollmueller's avatar
Jakob Knollmueller committed
17
18
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
    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
45
46
    """

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

    def at(self, position):
65
66
67
68
69
        return self.__class__(position, self.m, D=self.D,
                              alpha =self.alpha,
                              q=self.q,
                              sigma=self.sigma, w=self.w,
                              samples=self.samples)
70
71
72

    @property
    def value(self):
73
74
75
        energy = exp(-self.position).dot(self.q + self.w / 2., bare= True)
        energy += self.position.dot(self.alpha - 1. + self.rho / 2., bare=True)
        energy += 0.5 * self.position.dot(self.T(self.position))
76
77
78
79
        return energy.real

    @property
    def gradient(self):
80
81
82
83
        gradient = - self.theta.weight(-1)
        gradient += (self.alpha - 1. + self.rho / 2.).weight(-1)
        gradient +=  self.T(self.position)
        gradient.val = gradient.val.real
84
85
86
87
        return gradient

    @property
    def curvature(self):
88
        curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T)
89
90
91
        return curvature

    def _calculate_w(self, m, D, samples):
92
        w = Field(domain=self.position.domain, val=0. , dtype=m.dtype)
93
94
95
        if D is not None:
            for i in range(samples):
                posterior_sample = generate_posterior_sample(m, D)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
96
                projected_sample = posterior_sample.power_analyze(
97
                    logarithmic=self.position.domain[0].config["logarithmic"],
Jakob Knollmueller's avatar
Jakob Knollmueller committed
98
99
                    nbin= self.position.domain[0].config["nbin"])
                w += (projected_sample) * self.rho
100
101
            w /= float(samples)
        else:
Jakob Knollmueller's avatar
Jakob Knollmueller committed
102
            w = m.power_analyze(
103
                    logarithmic=self.position.domain[0].config["logarithmic"],
Jakob Knollmueller's avatar
Jakob Knollmueller committed
104
                        nbin=self.position.domain[0].config["nbin"])
Jakob Knollmueller's avatar
Jakob Knollmueller committed
105
            w *= self.rho
106

107
        return w
108
109
110