nonlinear_power_energy.py 4.91 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
from .. import exp
Philipp Arras's avatar
Philipp Arras committed
20 21
from ..utilities import memo
from ..operators.smoothness_operator import SmoothnessOperator
22 23
from .nonlinear_power_curvature import NonlinearPowerCurvature
from .response_operators import LinearizedPowerResponse
Philipp Arras's avatar
Philipp Arras committed
24
from ..minimization.energy import Energy
25 26 27 28 29


class NonlinearPowerEnergy(Energy):
    """The Energy of the power spectrum according to the critical filter.

Martin Reinecke's avatar
Martin Reinecke committed
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 37 38 39 40 41 42 43 44 45 46 47 48

    Parameters
    ----------
    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
    sigma : float
        The parameter of the smoothness prior.
        default : ??? None? ???????
    samples : integer
Martin Reinecke's avatar
Martin Reinecke committed
49 50
        Number of samples used for the estimation of the uncertainty
        corrections.
51 52 53
        default : 3
    """

Martin Reinecke's avatar
Martin Reinecke committed
54 55 56 57
    def __init__(self, position, d, N, m, D, FFT, Instrument, nonlinearity,
                 Projection, sigma=0., samples=3, sample_list=None,
                 inverter=None):
        super(NonlinearPowerEnergy, self).__init__(position)
58 59 60 61
        self.m = m
        self.D = D
        self.d = d
        self.N = N
Martin Reinecke's avatar
Martin Reinecke committed
62
        self.T = SmoothnessOperator(domain=self.position.domain[0],
Martin Reinecke's avatar
Martin Reinecke committed
63
                                    strength=sigma, logarithmic=True)
64 65 66 67
        self.FFT = FFT
        self.Instrument = Instrument
        self.nonlinearity = nonlinearity
        self.Projection = Projection
Martin Reinecke's avatar
Martin Reinecke committed
68
        self._sigma = sigma
69

Martin Reinecke's avatar
Martin Reinecke committed
70
        self.power = self.Projection.adjoint_times(exp(0.5*self.position))
71
        if sample_list is None:
Martin Reinecke's avatar
Martin Reinecke committed
72
            if samples is None or samples == 0:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
73
                sample_list = [m]
74
            else:
Martin Reinecke's avatar
Martin Reinecke committed
75
                sample_list = [D.generate_posterior_sample() + m
Martin Reinecke's avatar
Martin Reinecke committed
76
                               for _ in range(samples)]
77 78
        self.sample_list = sample_list
        self.inverter = inverter
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
79
        self._value, self._gradient = self._value_and_grad()
80 81

    def at(self, position):
Martin Reinecke's avatar
Martin Reinecke committed
82 83
        return self.__class__(position, self.d, self.N, self.m, self.D,
                              self.FFT, self.Instrument, self.nonlinearity,
Martin Reinecke's avatar
Martin Reinecke committed
84
                              self.Projection, sigma=self._sigma,
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
85
                              samples=len(self.sample_list),
86
                              sample_list=self.sample_list,
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
87
                              inverter=self.inverter)
88

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
89
    def _value_and_grad(self):
Martin Reinecke's avatar
Martin Reinecke committed
90
        likelihood_gradient = None
91
        for sample in self.sample_list:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
92 93 94
            residual = self.d - \
                self.Instrument(self.nonlinearity(
                    self.FFT.adjoint_times(self.power*sample)))
Martin Reinecke's avatar
Martin Reinecke committed
95
            lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
96 97 98 99
            LinR = LinearizedPowerResponse(
                self.Instrument, self.nonlinearity, self.FFT, self.Projection,
                self.position, sample)
            grad = LinR.adjoint_times(self.N.inverse_times(residual))
Martin Reinecke's avatar
Martin Reinecke committed
100
            if likelihood_gradient is None:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
101 102
                likelihood = lh
                likelihood_gradient = grad.copy()
Martin Reinecke's avatar
Martin Reinecke committed
103
            else:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
104 105 106 107 108 109 110 111
                likelihood += lh
                likelihood_gradient += grad
        Tpos = self.T(self.position)
        likelihood *= 1./len(self.sample_list)
        likelihood += 0.5*self.position.vdot(Tpos)
        likelihood_gradient *= -1./len(self.sample_list)
        likelihood_gradient += Tpos
        return likelihood, likelihood_gradient
112

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
113 114 115 116 117 118 119
    @property
    def value(self):
        return self._value

    @property
    def gradient(self):
        return self._gradient
120 121 122 123

    @property
    @memo
    def curvature(self):
124
        return NonlinearPowerCurvature(
Martin Reinecke's avatar
Martin Reinecke committed
125
            self.position, self.FFT, self.Instrument, self.nonlinearity,
126 127
            self.Projection, self.N, self.T, self.sample_list,
            inverter=self.inverter)