noise_energy.py 2.72 KB
Newer Older
Philipp Arras's avatar
Philipp Arras committed
1
2
3
from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..minimization.energy import Energy
4
5
6


class NoiseEnergy(Energy):
Martin Reinecke's avatar
Martin Reinecke committed
7
8
9
    def __init__(self, position, d, m, D, t, FFT, Instrument, nonlinearity,
                 alpha, q, Projection, samples=3, sample_list=None,
                 inverter=None):
Martin Reinecke's avatar
Martin Reinecke committed
10
        super(NoiseEnergy, self).__init__(position=position)
11
12
13
        self.m = m
        self.D = D
        self.d = d
Martin Reinecke's avatar
Martin Reinecke committed
14
        self.N = DiagonalOperator(diagonal=exp(self.position))
15
16
17
18
19
20
21
22
23
24
25
26
        self.t = t
        self.samples = samples
        self.FFT = FFT
        self.Instrument = Instrument
        self.nonlinearity = nonlinearity

        self.alpha = alpha
        self.q = q
        self.Projection = Projection
        self.power = self.Projection.adjoint_times(exp(0.5 * self.t))
        self.one = Field(self.position.domain, val=1.)
        if sample_list is None:
Martin Reinecke's avatar
Martin Reinecke committed
27
28
            if samples is None or samples == 0:
                sample_list = [m]
29
            else:
Martin Reinecke's avatar
Martin Reinecke committed
30
                sample_list = [D.generate_posterior_sample() + m
Martin Reinecke's avatar
Martin Reinecke committed
31
                               for _ in range(samples)]
32
33
        self.sample_list = sample_list
        self.inverter = inverter
Martin Reinecke's avatar
Martin Reinecke committed
34
        self._value, self._gradient = self._value_and_grad()
35
36

    def at(self, position):
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
40
41
        return self.__class__(
            position, self.d, self.m, self.D, self.t, self.FFT,
            self.Instrument, self.nonlinearity, self.alpha, self.q,
            self.Projection, sample_list=self.sample_list,
            samples=self.samples, inverter=self.inverter)
42

Martin Reinecke's avatar
Martin Reinecke committed
43
44
    def _value_and_grad(self):
        likelihood_gradient = None
45
        for sample in self.sample_list:
Martin Reinecke's avatar
Martin Reinecke committed
46
47
48
49
50
51
52
53
54
55
56
            residual = self.d - \
                self.Instrument(self.nonlinearity(
                    self.FFT.adjoint_times(self.power*sample)))
            lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
            grad = -0.5 * self.N.inverse_times(residual.conjugate() * residual)
            if likelihood_gradient is None:
                likelihood = lh
                likelihood_gradient = grad.copy()
            else:
                likelihood += lh
                likelihood_gradient += grad
57

Martin Reinecke's avatar
Martin Reinecke committed
58
59
60
61
62
63
64
65
66
        likelihood = ((likelihood / float(len(self.sample_list))) +
                      0.5 * self.position.integrate() +
                      (self.alpha - 1.).vdot(self.position) +
                      self.q.vdot(exp(-self.position)))
        likelihood_gradient = (
            likelihood_gradient / float(len(self.sample_list)) +
            (self.alpha-0.5) -
            self.q * (exp(-self.position)))
        return likelihood, likelihood_gradient
67
68

    @property
Martin Reinecke's avatar
Martin Reinecke committed
69
70
    def value(self):
        return self._value
71
72

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