critical_power_energy.py 5.04 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 38 39
    position : Field,
        The current position of this energy.
    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
        energy = self._theta.integrate()
Martin Reinecke's avatar
Martin Reinecke committed
105 106 107 108
        energy += self.position.integrate()*(self.alpha-0.5)
        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)