critical_power_energy.py 4.6 KB
Newer Older
1 2
import numpy as np

3
from nifty.energies.energy import Energy
4 5
from nifty.operators.smoothness_operator import SmoothnessOperator
from nifty.library.critical_filter import CriticalPowerCurvature
6
from nifty.energies.memoization import memo
Jakob Knollmueller's avatar
Jakob Knollmueller committed
7

8
from nifty.sugar import generate_posterior_sample
9
from nifty import Field, exp
10

11

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

15 16 17 18
    It describes the energy of the logarithmic amplitudes of the power spectrum
    of a Gaussian random field under reconstruction uncertainty with smoothness
    and inverse gamma prior. It is used to infer the correlation structure of a
    correlated signal.
19 20 21

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

56 57 58
    def __init__(self, position, m, D=None, alpha=1.0, q=0.,
                 smoothness_prior=0., logarithmic=True, samples=3, w=None):
        super(CriticalPowerEnergy, self).__init__(position=position)
59
        self.m = m
60 61
        self.D = D
        self.samples = samples
62
        self.smoothness_prior = np.float(smoothness_prior)
63 64
        self.alpha = Field(self.position.domain, val=alpha)
        self.q = Field(self.position.domain, val=q)
65 66
        self.T = SmoothnessOperator(domain=self.position.domain[0],
                                    strength=self.smoothness_prior,
Jakob Knollmueller's avatar
Jakob Knollmueller committed
67
                                    logarithmic=logarithmic)
68
        self.rho = self.position.domain[0].rho
69
        self._w = w if w is not None else None
70

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

    @property
    def value(self):
78
        energy = self._theta.sum()
79 80
        energy += self.position.vdot(self._rho_prime, bare=True)
        energy += 0.5 * self.position.vdot(self._Tt)
81 82 83 84
        return energy.real

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

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

97 98 99 100 101 102 103 104 105
    @property
    def w(self):
        if self._w is None:
            w = Field(domain=self.position.domain, val=0., dtype=self.m.dtype)
            if self.D is not None:
                for i in range(self.samples):
                    posterior_sample = generate_posterior_sample(
                                                            self.m, self.D)
                    projected_sample = posterior_sample.power_analyze(
Martin Reinecke's avatar
Martin Reinecke committed
106
                     binbounds=self.position.domain[0].binbounds)
107 108 109 110
                    w += (projected_sample) * self.rho
                w /= float(self.samples)
            else:
                w = self.m.power_analyze(
Martin Reinecke's avatar
Martin Reinecke committed
111
                     binbounds=self.position.domain[0].binbounds)
112 113 114
                w *= self.rho
            self._w = w
        return self._w
115

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

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

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