nonlinear_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
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
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
    """

Philipp Arras's avatar
Philipp Arras committed
54
    def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
Philipp Arras's avatar
Philipp Arras committed
55
                 Projection, sigma=0., samples=3, xi_sample_list=None,
Philipp Arras's avatar
Philipp Arras committed
56
                 inverter=None, munit=1., sunit=1.):
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
67
        self.Instrument = Instrument
        self.nonlinearity = nonlinearity
        self.Projection = Projection
68
        self.sigma = sigma
Philipp Arras's avatar
Philipp Arras committed
69
70
        self.munit = munit
        self.sunit = sunit
Philipp Arras's avatar
Philipp Arras committed
71
        if xi_sample_list is None:
Martin Reinecke's avatar
Martin Reinecke committed
72
            if samples is None or samples == 0:
Philipp Arras's avatar
Philipp Arras committed
73
                xi_sample_list = [xi]
74
            else:
Philipp Arras's avatar
Philipp Arras committed
75
                xi_sample_list = [D.generate_posterior_sample() + xi
Philipp Arras's avatar
Philipp Arras committed
76
                                  for _ in range(samples)]
Philipp Arras's avatar
Philipp Arras committed
77
        self.xi_sample_list = xi_sample_list
78
        self.inverter = inverter
79

Philipp Arras's avatar
Philipp Arras committed
80
        A = Projection.adjoint_times(munit * exp(.5 * position))  # unit: munit
Philipp Arras's avatar
Philipp Arras committed
81
        map_s = self.ht(A * xi)
82
83
84
        Tpos = self.T(position)

        self._gradient = None
Philipp Arras's avatar
Philipp Arras committed
85
        for xi_sample in self.xi_sample_list:
Philipp Arras's avatar
Philipp Arras committed
86
            map_s = self.ht(A * xi_sample)
Philipp Arras's avatar
Philipp Arras committed
87
            LinR = LinearizedPowerResponse(
Philipp Arras's avatar
Philipp Arras committed
88
                self.Instrument, self.nonlinearity, self.ht, self.Projection,
Philipp Arras's avatar
Philipp Arras committed
89
                self.position, xi_sample, munit=self.munit, sunit=self.sunit)
Philipp Arras's avatar
Philipp Arras committed
90
91

            residual = self.d - \
Philipp Arras's avatar
Philipp Arras committed
92
                self.Instrument(sunit * self.nonlinearity(map_s))
93
94
95
96
97
98
99
100
101
102
            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
103
        self._value *= 1. / len(self.xi_sample_list)
Philipp Arras's avatar
Philipp Arras committed
104
        self._value += 0.5 * self.position.vdot(Tpos)
Philipp Arras's avatar
Philipp Arras committed
105
        self._gradient *= -1. / len(self.xi_sample_list)
106
        self._gradient += Tpos
107
108

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

Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
118
119
120
121
122
123
124
    @property
    def value(self):
        return self._value

    @property
    def gradient(self):
        return self._gradient
125
126
127
128

    @property
    @memo
    def curvature(self):
129
        return NonlinearPowerCurvature(
130
            self.position, self.ht, self.Instrument, self.nonlinearity,
Philipp Arras's avatar
Philipp Arras committed
131
            self.Projection, self.N, self.T, self.xi_sample_list,
Philipp Arras's avatar
Philipp Arras committed
132
            self.inverter, self.munit, self.sunit)