noise_energy.py 3.19 KB
Newer Older
Philipp Arras's avatar
Philipp Arras committed
1
2
3
4
5
from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..sugar import generate_posterior_sample
from ..minimization.energy import Energy
from ..utilities import memo
6
7
8
9
10
11
12
13
14
15


class NoiseEnergy(Energy):
    def __init__(self, position, d, m, D, t, FFT, Instrument, nonlinearity, alpha, q, Projection,
                 samples=3, sample_list=None, inverter=None):
        super(NoiseEnergy, self).__init__(position=position.copy())
        dummy = self.position.norm()
        self.m = m
        self.D = D
        self.d = d
Martin Reinecke's avatar
Martin Reinecke committed
16
        self.N = DiagonalOperator(diagonal=exp(self.position))
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
        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):
                    sample = generate_posterior_sample(m, D)
                    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):
        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)

    @property
    @memo
    def value(self):
        likelihood = 0.
        for sample in self.sample_list:
            likelihood += self._likelihood(sample)
        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)))

    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)
        return (likelihood_gradient / float(len(self.sample_list))
                + 0.5 * self.one + (self.alpha - self.one) - self.q * (exp(-self.position)))

    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