noise_energy.py 3.04 KB
 Philipp Arras committed Nov 21, 2017 1 2 3 4 ``````from .. import Field, exp from ..operators.diagonal_operator import DiagonalOperator from ..minimization.energy import Energy from ..utilities import memo `````` Philipp Arras committed Nov 21, 2017 5 6 7 `````` class NoiseEnergy(Energy): `````` Martin Reinecke committed Nov 29, 2017 8 9 10 `````` def __init__(self, position, d, m, D, t, FFT, Instrument, nonlinearity, alpha, q, Projection, samples=3, sample_list=None, inverter=None): `````` Philipp Arras committed Nov 21, 2017 11 12 13 14 `````` super(NoiseEnergy, self).__init__(position=position.copy()) self.m = m self.D = D self.d = d `````` Martin Reinecke committed Nov 22, 2017 15 `````` self.N = DiagonalOperator(diagonal=exp(self.position)) `````` Philipp Arras committed Nov 21, 2017 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) `````` Philipp Arras committed Nov 21, 2017 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 committed Nov 29, 2017 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) `````` Philipp Arras committed Nov 21, 2017 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 committed Nov 29, 2017 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))) `````` Philipp Arras committed Nov 21, 2017 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 committed Nov 29, 2017 71 72 73 `````` return (likelihood_gradient / float(len(self.sample_list)) + 0.5 * self.one + (self.alpha - self.one) - self.q * (exp(-self.position))) `````` Philipp Arras committed Nov 21, 2017 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``````