ensemble_likelihood.py 2.55 KB
Newer Older
1 2
# -*- coding: utf-8 -*-

Theo Steininger's avatar
Theo Steininger committed
3 4
import numpy as np

5

6 7 8 9
class EnsembleLikelihood(object):
    def __init__(self, observable_name,  measured_data,
                 data_covariance_operator):
        self.observable_name = observable_name
10
        self.measured_data = measured_data
Theo Steininger's avatar
Theo Steininger committed
11
        self.data_covariance_operator = data_covariance_operator
12 13

    def __call__(self, observable):
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's avatar
Theo Steininger committed
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

27 28
        observable = field

Theo Steininger's avatar
Theo Steininger committed
29 30
        k = observable.shape[0]

31
        A = data_covariance_operator
Theo Steininger's avatar
Theo Steininger committed
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):
48
            c = measured_data - obs_val[i]
Theo Steininger's avatar
Theo Steininger committed
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()