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

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

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

    Parameters
    ----------
Jakob Knollmueller's avatar
Jakob Knollmueller committed
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
45
    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
46
47
48
    inverter : ConjugateGradient
        The inversion strategy to invert the curvature and to generate samples.
        default : None
49
50
    """

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

69

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

    @property
    def value(self):
Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
        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))
82
83
84
85
        return energy.real

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

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

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

113
        return w
114
115
116