critical_power_energy.py 5.1 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
from ..minimization.energy import Energy
from ..operators.smoothness_operator import SmoothnessOperator
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..operators.power_distributor import PowerDistributor
Martin Reinecke's avatar
Martin Reinecke committed
22
23
24
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,
Philipp Arras's avatar
Philipp Arras committed
40
41
        The map whose power spectrum is inferred. Needs to live in harmonic
        signal space.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
42
43
44
45
46
    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
47
48
        The spectral prior of the inverse gamma distribution. 1.0 corresponds
        to non-informative.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
49
50
        default : 1.0
    q : float
51
52
53
54
55
        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
56
57
58
59
        default : 0.0
    logarithmic : boolean
        Whether smoothness acts on linear or logarithmic scale.
    samples : integer
60
61
        Number of samples used for the estimation of the uncertainty
        corrections.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
62
63
64
65
66
        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
67
68
69
    inverter : ConjugateGradient
        The inversion strategy to invert the curvature and to generate samples.
        default : None
70
71
    """

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

Martin Reinecke's avatar
Martin Reinecke committed
88
        if w is None:
Martin Reinecke's avatar
Martin Reinecke committed
89
90
            Dist = PowerDistributor(target=self.m.domain,
                                    power_space=self.position.domain[0])
Martin Reinecke's avatar
Martin Reinecke committed
91
92
93
            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
94
                    sample = self.D.draw_sample() + self.m
Martin Reinecke's avatar
Martin Reinecke committed
95
                    w += Dist.adjoint_times(abs(sample)**2)
Martin Reinecke's avatar
Martin Reinecke committed
96

97
                w *= 1./self.samples
Martin Reinecke's avatar
Martin Reinecke committed
98
            else:
Martin Reinecke's avatar
Martin Reinecke committed
99
                w = Dist.adjoint_times(abs(self.m)**2)
Martin Reinecke's avatar
Martin Reinecke committed
100
101
        self._w = w

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

Philipp Arras's avatar
Philipp Arras committed
105
        energy = self._theta.sum()
106
107
        energy += self.position.sum() * (self.alpha-0.5)
        energy += 0.5*self.position.vdot(Tt)
Martin Reinecke's avatar
Martin Reinecke committed
108
109
        self._value = energy.real

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

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

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

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

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