nonlinear_power_energy.py 4.08 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
from .. import exp
Philipp Arras's avatar
Philipp Arras committed
2 3
from ..utilities import memo
from ..operators.smoothness_operator import SmoothnessOperator
4 5
from .nonlinear_power_curvature import NonlinearPowerCurvature
from .response_operators import LinearizedPowerResponse
Philipp Arras's avatar
Philipp Arras committed
6
from ..minimization.energy import Energy
7 8 9 10 11


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

Martin Reinecke's avatar
Martin Reinecke committed
12 13 14 15
    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.
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

    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
31 32
        Number of samples used for the estimation of the uncertainty
        corrections.
33 34 35
        default : 3
    """

Martin Reinecke's avatar
Martin Reinecke committed
36 37 38 39
    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)
40 41 42 43
        self.m = m
        self.D = D
        self.d = d
        self.N = N
Martin Reinecke's avatar
Martin Reinecke committed
44
        self.T = SmoothnessOperator(domain=self.position.domain[0],
Martin Reinecke's avatar
Martin Reinecke committed
45
                                    strength=sigma, logarithmic=True)
46 47 48 49 50
        self.FFT = FFT
        self.Instrument = Instrument
        self.nonlinearity = nonlinearity
        self.Projection = Projection

Martin Reinecke's avatar
Martin Reinecke committed
51
        self.power = self.Projection.adjoint_times(exp(0.5*self.position))
52
        if sample_list is None:
Martin Reinecke's avatar
Martin Reinecke committed
53
            if samples is None or samples == 0:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
54
                sample_list = [m]
55
            else:
Martin Reinecke's avatar
Martin Reinecke committed
56
                sample_list = [D.generate_posterior_sample() + m
Martin Reinecke's avatar
Martin Reinecke committed
57
                               for _ in range(samples)]
58 59
        self.sample_list = sample_list
        self.inverter = inverter
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
60
        self._value, self._gradient = self._value_and_grad()
61 62

    def at(self, position):
Martin Reinecke's avatar
Martin Reinecke committed
63 64
        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
65
                              self.Projection, sigma=self.T.strength,
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
66
                              samples=len(self.sample_list),
67
                              sample_list=self.sample_list,
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
68
                              inverter=self.inverter)
69

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
70
    def _value_and_grad(self):
Martin Reinecke's avatar
Martin Reinecke committed
71
        likelihood_gradient = None
72
        for sample in self.sample_list:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
73 74 75
            residual = self.d - \
                self.Instrument(self.nonlinearity(
                    self.FFT.adjoint_times(self.power*sample)))
Martin Reinecke's avatar
Martin Reinecke committed
76
            lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
77 78 79 80
            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
81
            if likelihood_gradient is None:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
82 83
                likelihood = lh
                likelihood_gradient = grad.copy()
Martin Reinecke's avatar
Martin Reinecke committed
84
            else:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
85 86 87 88 89 90 91 92
                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
93

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
94 95 96 97 98 99 100
    @property
    def value(self):
        return self._value

    @property
    def gradient(self):
        return self._gradient
101 102 103 104

    @property
    @memo
    def curvature(self):
105
        return NonlinearPowerCurvature(
Martin Reinecke's avatar
Martin Reinecke committed
106
            self.position, self.FFT, self.Instrument, self.nonlinearity,
107 108
            self.Projection, self.N, self.T, self.sample_list,
            inverter=self.inverter)