critical_power_energy.py 4.47 KB
Newer Older
1
from nifty.energies.energy import Energy
2
3
4
from nifty.operators.smoothness_operator import SmoothnessOperator
from nifty.library.critical_filter import CriticalPowerCurvature

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):
Martin Reinecke's avatar
Martin Reinecke committed
80
81
82
        energy = exp(-self.position).vdot(self.q + self.w / 2., bare= True)
        energy += self.position.vdot(self.alpha - 1. + self.rho / 2., bare=True)
        energy += 0.5 * self.position.vdot(self.T(self.position))
83
84
85
86
        return energy.real

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

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

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

114
        return w
115
116
117