critical_power_energy.py 5.12 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

Martin Reinecke's avatar
Martin Reinecke committed
19
20
21
22
23
24
from ..minimization.energy import Energy
from ..operators.smoothness_operator import SmoothnessOperator
from ..operators.power_projection_operator import PowerProjectionOperator
from .critical_power_curvature import CriticalPowerCurvature
from ..utilities import memo
from .. import Field, exp
25

26

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

30
31
32
33
    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.
34
35
36

    Parameters
    ----------
Jakob Knollmueller's avatar
Jakob Knollmueller committed
37
    position : Field,
38
        The current position of this energy. (Logarithm of power spectrum)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
39
    m : Field,
Martin Reinecke's avatar
Martin Reinecke committed
40
        The map whose power spectrum has to be inferred
Jakob Knollmueller's avatar
Jakob Knollmueller committed
41
42
43
44
45
    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
46
47
        The spectral prior of the inverse gamma distribution. 1.0 corresponds
        to non-informative.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
48
49
        default : 1.0
    q : float
50
51
52
53
54
        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
55
56
57
58
        default : 0.0
    logarithmic : boolean
        Whether smoothness acts on linear or logarithmic scale.
    samples : integer
59
60
        Number of samples used for the estimation of the uncertainty
        corrections.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
61
62
63
64
65
        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
66
67
68
    inverter : ConjugateGradient
        The inversion strategy to invert the curvature and to generate samples.
        default : None
69
70
    """

71
    def __init__(self, position, m, D=None, alpha=1.0, q=0.,
72
73
                 smoothness_prior=0., logarithmic=True, samples=3, w=None,
                 inverter=None):
74
        super(CriticalPowerEnergy, self).__init__(position=position)
75
        self.m = m
76
77
        self.D = D
        self.samples = samples
Martin Reinecke's avatar
Martin Reinecke committed
78
79
        self.alpha = float(alpha)
        self.q = float(q)
Martin Reinecke's avatar
Martin Reinecke committed
80
81
        self._smoothness_prior = smoothness_prior
        self._logarithmic = logarithmic
82
        self.T = SmoothnessOperator(domain=self.position.domain[0],
83
                                    strength=smoothness_prior,
Jakob Knollmueller's avatar
Jakob Knollmueller committed
84
                                    logarithmic=logarithmic)
85
        self._inverter = inverter
86

Martin Reinecke's avatar
Martin Reinecke committed
87
88
89
90
91
92
        if w is None:
            P = PowerProjectionOperator(domain=self.m.domain,
                                        power_space=self.position.domain[0])
            if self.D is not None:
                w = Field.zeros(self.position.domain, dtype=self.m.dtype)
                for i in range(self.samples):
Martin Reinecke's avatar
Martin Reinecke committed
93
                    sample = self.D.generate_posterior_sample() + self.m
Martin Reinecke's avatar
Martin Reinecke committed
94
95
96
97
98
99
100
                    w += P(abs(sample)**2)

                w *= 1./self.samples
            else:
                w = P(abs(self.m)**2)
        self._w = w

101
        self._theta = exp(-self.position) * (self.q + self._w*0.5)
Martin Reinecke's avatar
Martin Reinecke committed
102
103
        Tt = self.T(self.position)

104
105
        energy = self._theta.vdot(Field.ones_like(self._theta))
        energy += self.position.vdot(Field.ones_like(self.position)) *(self.alpha-0.5)
Martin Reinecke's avatar
Martin Reinecke committed
106
107
108
        energy += 0.5*self.position.vdot(Tt)
        self._value = energy.real

109
        gradient = -self._theta
Martin Reinecke's avatar
Martin Reinecke committed
110
111
        gradient += self.alpha-0.5
        gradient += Tt
112
        self._gradient = gradient
Martin Reinecke's avatar
Martin Reinecke committed
113

114
    def at(self, position):
115
        return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
Martin Reinecke's avatar
Martin Reinecke committed
116
117
118
                              q=self.q,
                              smoothness_prior=self._smoothness_prior,
                              logarithmic=self._logarithmic,
119
                              samples=self.samples, w=self._w,
120
                              inverter=self._inverter)
121
122
123

    @property
    def value(self):
Martin Reinecke's avatar
Martin Reinecke committed
124
        return self._value
125
126
127

    @property
    def gradient(self):
Martin Reinecke's avatar
Martin Reinecke committed
128
        return self._gradient
129
130

    @property
Martin Reinecke's avatar
Martin Reinecke committed
131
    @memo
132
    def curvature(self):
133
134
        return CriticalPowerCurvature(theta=self._theta, T=self.T,
                                      inverter=self._inverter)