critical_power_energy.py 4.46 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
47
    inverter : ConjugateGradient
        The inversion strategy to invert the curvature and to generate samples.
        default : None
48
49
    """

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

68

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

    @property
    def value(self):
78
79
80
        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))
81
82
83
84
        return energy.real

    @property
    def gradient(self):
85
86
87
88
        gradient = - self.theta.weight(-1)
        gradient += (self.alpha - 1. + self.rho / 2.).weight(-1)
        gradient +=  self.T(self.position)
        gradient.val = gradient.val.real
89
90
91
92
        return gradient

    @property
    def curvature(self):
93
        curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T, inverter=self.inverter)
94
95
96
        return curvature

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

112
        return w
113
114
115