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


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

    Parameters
    ----------
Martin Reinecke's avatar
Martin Reinecke committed
37
    position : Field
38
        The current position of this energy.
Martin Reinecke's avatar
Martin Reinecke committed
39
    m : Field
40
        The map whichs power spectrum has to be inferred
Martin Reinecke's avatar
Martin Reinecke committed
41
    D : EndomorphicOperator
42
43
44
45
46
47
        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? ???????
Martin Reinecke's avatar
Martin Reinecke committed
48
    samples : int
Martin Reinecke's avatar
Martin Reinecke committed
49
50
        Number of samples used for the estimation of the uncertainty
        corrections.
51
52
        default : 3
    """
Martin Reinecke's avatar
Martin Reinecke committed
53
    # MR FIXME: docstring incomplete and outdated
Philipp Arras's avatar
Philipp Arras committed
54
    def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
Martin Reinecke's avatar
Martin Reinecke committed
55
                 Distributor, sigma=0., samples=3, xi_sample_list=None,
Philipp Arras's avatar
Philipp Arras committed
56
                 inverter=None):
Martin Reinecke's avatar
Martin Reinecke committed
57
        super(NonlinearPowerEnergy, self).__init__(position)
Philipp Arras's avatar
Philipp Arras committed
58
        self.xi = xi
59
60
61
        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
        self.ht = ht
65
66
        self.Instrument = Instrument
        self.nonlinearity = nonlinearity
Martin Reinecke's avatar
Martin Reinecke committed
67
        self.Distributor = Distributor
68
        self.sigma = sigma
Philipp Arras's avatar
Philipp Arras committed
69
        if xi_sample_list is None:
Martin Reinecke's avatar
Martin Reinecke committed
70
            if samples is None or samples == 0:
Philipp Arras's avatar
Philipp Arras committed
71
                xi_sample_list = [xi]
72
            else:
Martin Reinecke's avatar
Martin Reinecke committed
73
                xi_sample_list = [D.draw_sample() + xi
Philipp Arras's avatar
Philipp Arras committed
74
                                  for _ in range(samples)]
Philipp Arras's avatar
Philipp Arras committed
75
        self.xi_sample_list = xi_sample_list
76
        self.inverter = inverter
77

Martin Reinecke's avatar
Martin Reinecke committed
78
        A = Distributor(exp(.5 * position))
Philipp Arras's avatar
Philipp Arras committed
79
        map_s = self.ht(A * xi)
80
81
82
        Tpos = self.T(position)

        self._gradient = None
Philipp Arras's avatar
Philipp Arras committed
83
        for xi_sample in self.xi_sample_list:
Philipp Arras's avatar
Philipp Arras committed
84
            map_s = self.ht(A * xi_sample)
Philipp Arras's avatar
Philipp Arras committed
85
            LinR = LinearizedPowerResponse(
Martin Reinecke's avatar
Martin Reinecke committed
86
                self.Instrument, self.nonlinearity, self.ht, self.Distributor,
Philipp Arras's avatar
Philipp Arras committed
87
                self.position, xi_sample)
Philipp Arras's avatar
Philipp Arras committed
88
89

            residual = self.d - \
Philipp Arras's avatar
Philipp Arras committed
90
                self.Instrument(self.nonlinearity(map_s))
91
92
93
94
95
96
97
98
99
100
            lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
            grad = LinR.adjoint_times(self.N.inverse_times(residual))

            if self._gradient is None:
                self._value = lh
                self._gradient = grad.copy()
            else:
                self._value += lh
                self._gradient += grad

Philipp Arras's avatar
Philipp Arras committed
101
        self._value *= 1. / len(self.xi_sample_list)
Philipp Arras's avatar
Philipp Arras committed
102
        self._value += 0.5 * self.position.vdot(Tpos)
Philipp Arras's avatar
Philipp Arras committed
103
        self._gradient *= -1. / len(self.xi_sample_list)
104
        self._gradient += Tpos
Martin Reinecke's avatar
Martin Reinecke committed
105
        self._gradient.lock()
106
107

    def at(self, position):
Philipp Arras's avatar
Philipp Arras committed
108
        return self.__class__(position, self.d, self.N, self.xi, self.D,
109
                              self.ht, self.Instrument, self.nonlinearity,
Martin Reinecke's avatar
Martin Reinecke committed
110
                              self.Distributor, sigma=self.sigma,
Philipp Arras's avatar
Philipp Arras committed
111
112
                              samples=len(self.xi_sample_list),
                              xi_sample_list=self.xi_sample_list,
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
113
                              inverter=self.inverter)
114

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

    @property
    def gradient(self):
        return self._gradient
122
123
124
125

    @property
    @memo
    def curvature(self):
126
        return NonlinearPowerCurvature(
127
            self.position, self.ht, self.Instrument, self.nonlinearity,
Martin Reinecke's avatar
Martin Reinecke committed
128
            self.Distributor, self.N, self.T, self.xi_sample_list,
Philipp Arras's avatar
Philipp Arras committed
129
            self.inverter)