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


class NoiseEnergy(Energy):
Martin Reinecke's avatar
Martin Reinecke committed
8
9
10
    def __init__(self, position, d, m, D, t, FFT, Instrument, nonlinearity,
                 alpha, q, Projection, samples=3, sample_list=None,
                 inverter=None):
11
12
13
14
        super(NoiseEnergy, self).__init__(position=position.copy())
        self.m = m
        self.D = D
        self.d = d
Martin Reinecke's avatar
Martin Reinecke committed
15
        self.N = DiagonalOperator(diagonal=exp(self.position))
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
        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:
            sample_list = []
            if samples is None:
                sample_list.append(self.m)
            else:
                for i in range(samples):
33
                    sample = D.generate_posterior_sample(m)
34
35
36
37
38
39
40
                    sample = FFT(Field(FFT.domain, val=(
                        FFT.adjoint_times(sample).val)))
                    sample_list.append(sample)
        self.sample_list = sample_list
        self.inverter = inverter

    def at(self, position):
Martin Reinecke's avatar
Martin Reinecke committed
41
42
43
44
45
        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)
46
47
48
49
50
51
52

    @property
    @memo
    def value(self):
        likelihood = 0.
        for sample in self.sample_list:
            likelihood += self._likelihood(sample)
Martin Reinecke's avatar
Martin Reinecke committed
53
54
55
56
        return ((likelihood / float(len(self.sample_list))) +
                0.5 * self.one.vdot(self.position) +
                (self.alpha - self.one).vdot(self.position) +
                self.q.vdot(exp(-self.position)))
57
58
59
60
61
62
63
64
65
66
67
68
69
70

    def _likelihood(self, m):
        residual = self.d - \
            self.Instrument(self.nonlinearity(
                self.FFT.adjoint_times(self.power * m)))
        energy = 0.5 * residual.vdot(self.N.inverse_times(residual))
        return energy.real

    @property
    @memo
    def gradient(self):
        likelihood_gradient = Field(self.position.domain, val=0.)
        for sample in self.sample_list:
            likelihood_gradient += self._likelihood_gradient(sample)
Martin Reinecke's avatar
Martin Reinecke committed
71
72
73
        return (likelihood_gradient / float(len(self.sample_list)) +
                0.5 * self.one + (self.alpha - self.one) -
                self.q * (exp(-self.position)))
74
75
76
77
78
79
80
81
82
83
84
85
86

    def _likelihood_gradient(self, m):
        residual = self.d - \
            self.Instrument(self.nonlinearity(
                self.FFT.adjoint_times(self.power * m)))
        gradient = -  0.5 * \
            self.N.inverse_times(residual.conjugate() * residual)
        return gradient

    @property
    @memo
    def curvature(self):
        pass