ensemble_likelihood.py 2.55 KB
 Theo Steininger committed Feb 15, 2017 1 2 ``````# -*- coding: utf-8 -*- `````` Theo Steininger committed Feb 21, 2017 3 4 ``````import numpy as np `````` Theo Steininger committed Feb 15, 2017 5 `````` `````` Theo Steininger committed Feb 21, 2017 6 7 8 9 ``````class EnsembleLikelihood(object): def __init__(self, observable_name, measured_data, data_covariance_operator): self.observable_name = observable_name `````` Theo Steininger committed Feb 15, 2017 10 `````` self.measured_data = measured_data `````` Theo Steininger committed Feb 21, 2017 11 `````` self.data_covariance_operator = data_covariance_operator `````` Theo Steininger committed Feb 15, 2017 12 13 `````` def __call__(self, observable): `````` Theo Steininger committed Feb 21, 2017 14 15 16 17 18 19 20 `````` field = observable[self.observable_name] return self._process_simple_field(field, self.measured_data, self.data_covariance_operator) def _process_simple_field(self, field, measured_data, data_covariance_operator): `````` Theo Steininger committed Feb 21, 2017 21 22 23 24 25 26 `````` # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization # B = A^{-1} + U U^dagger # A = data_covariance # B^{-1} c = (A_inv - # A_inv U (I_k + U^dagger A_inv U)^{-1} U^dagger A_inv) c `````` Theo Steininger committed Feb 21, 2017 27 28 `````` observable = field `````` Theo Steininger committed Feb 21, 2017 29 30 `````` k = observable.shape[0] `````` Theo Steininger committed Feb 21, 2017 31 `````` A = data_covariance_operator `````` Theo Steininger committed Feb 21, 2017 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 `````` obs_val = observable.val.get_full_data() obs_mean = observable.mean(spaces=0).val.get_full_data() u_val = observable.val.get_full_data() - obs_mean U = observable.copy_empty() U.val = u_val a_u = A.inverse_times(U, spaces=1) # build middle-matrix (kxk) a_u_val = a_u.val.get_full_data() middle = (np.eye(k) + np.einsum(u_val.conjugate(), [0, 1], a_u_val, [2, 1])) middle = np.linalg.inv(middle) result_array = np.zeros(k) for i in xrange(k): `````` Theo Steininger committed Feb 21, 2017 48 `````` c = measured_data - obs_val[i] `````` Theo Steininger committed Feb 21, 2017 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 `````` # assuming that A == A^dagger, this can be shortend # a_c = A.inverse_times(c) # u_a_c = a_c.dot(U, spaces=1) # u_a_c = u_a_c.conjugate() # and: double conjugate shouldn't make a difference # u_a_c = c.conjugate().dot(a_u, spaces=1).conjugate() u_a_c = c.dot(a_u, spaces=1) u_a_c_val = u_a_c.val.get_full_data() first_summand = A.inverse_times(c) second_summand_val = np.einsum(middle, [0, 1], u_a_c_val, [1]) second_summand_val = np.einsum(a_u_val, [0, 1], second_summand_val, [0]) second_summand = first_summand.copy_empty() second_summand.val = second_summand_val result = c.dot(first_summand - second_summand) result_array[i] = result return -result_array.mean()``````