noise_energy.py 3.56 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.

Philipp Arras's avatar
Philipp Arras committed
19
20
21
from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..minimization.energy import Energy
22
23
24


class NoiseEnergy(Energy):
Martin Reinecke's avatar
Martin Reinecke committed
25
    def __init__(self, position, d, m, D, t, FFT, Instrument, nonlinearity,
Philipp Arras's avatar
Philipp Arras committed
26
                 alpha, q, Projection, munit=1., sunit=1., dunit=1., samples=3, sample_list=None,
Martin Reinecke's avatar
Martin Reinecke committed
27
                 inverter=None):
Martin Reinecke's avatar
Martin Reinecke committed
28
        super(NoiseEnergy, self).__init__(position=position)
29
30
31
        self.m = m
        self.D = D
        self.d = d
Philipp Arras's avatar
Philipp Arras committed
32
        self.N = DiagonalOperator(diagonal=dunit**2 * exp(self.position))
33
34
35
36
37
        self.t = t
        self.samples = samples
        self.FFT = FFT
        self.Instrument = Instrument
        self.nonlinearity = nonlinearity
Philipp Arras's avatar
Philipp Arras committed
38
39
40
        self.munit = munit
        self.sunit = sunit
        self.dunit = dunit
41
42
43
44

        self.alpha = alpha
        self.q = q
        self.Projection = Projection
Philipp Arras's avatar
Philipp Arras committed
45
        self.power = self.Projection.adjoint_times(munit * exp(0.5 * self.t))
46
47
        self.one = Field(self.position.domain, val=1.)
        if sample_list is None:
Martin Reinecke's avatar
Martin Reinecke committed
48
49
            if samples is None or samples == 0:
                sample_list = [m]
50
            else:
Martin Reinecke's avatar
Martin Reinecke committed
51
                sample_list = [D.generate_posterior_sample() + m
Martin Reinecke's avatar
Martin Reinecke committed
52
                               for _ in range(samples)]
53
54
        self.sample_list = sample_list
        self.inverter = inverter
Philipp Arras's avatar
Philipp Arras committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

        A = Projection.adjoint_times(munit * exp(.5*self.t)) # unit: munit
        map_s = FFT.inverse_times(A * m)

        self._gradient = None
        for sample in self.sample_list:
            map_s = FFT.inverse_times(A * sample)

            residual = self.d - self.Instrument(sunit * self.nonlinearity(map_s))
            lh = .5 * residual.vdot(self.N.inverse_times(residual))
            grad = -.5 * self.N.inverse_times(residual.conjugate() * residual)

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

        self._value *= 1./len(self.sample_list)
        self._value += .5 * self.position.integrate() + (self.alpha - 1.).vdot(self.position) + self.q.vdot(exp(-self.position))

        self._gradient *= 1./len(self.sample_list)
        self._gradient += (self.alpha-0.5) - self.q * (exp(-self.position))
79
80

    def at(self, position):
Martin Reinecke's avatar
Martin Reinecke committed
81
82
83
        return self.__class__(
            position, self.d, self.m, self.D, self.t, self.FFT,
            self.Instrument, self.nonlinearity, self.alpha, self.q,
Philipp Arras's avatar
Philipp Arras committed
84
85
            self.Projection, munit=self.munit, sunit=self.sunit,
            dunit=self.dunit, sample_list=self.sample_list,
Martin Reinecke's avatar
Martin Reinecke committed
86
            samples=self.samples, inverter=self.inverter)
87
88

    @property
Martin Reinecke's avatar
Martin Reinecke committed
89
90
    def value(self):
        return self._value
91
92

    @property
Martin Reinecke's avatar
Martin Reinecke committed
93
94
    def gradient(self):
        return self._gradient