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

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

5 6
from nifty import DiagonalOperator

7 8
from imagine.likelihoods.likelihood import Likelihood
from imagine.create_ring_profile import create_ring_profile
9

10 11

class EnsembleLikelihood(Likelihood):
12
    def __init__(self, observable_name,  measured_data,
13
                 data_covariance_operator, profile=None):
14
        self.observable_name = observable_name
15
        self.measured_data = measured_data
Theo Steininger's avatar
Theo Steininger committed
16
        self.data_covariance_operator = data_covariance_operator
17 18
        if profile is None:
            profile = create_ring_profile(
19
                            self.measured_data.val.get_full_data())
20
        self.profile = profile
21 22

    def __call__(self, observable):
23 24 25
        field = observable[self.observable_name]
        return self._process_simple_field(field,
                                          self.measured_data,
26 27
                                          self.data_covariance_operator,
                                          self.profile)
28

29
    def _process_simple_field(self, observable, measured_data,
30
                              data_covariance_operator, profile):
Theo Steininger's avatar
Theo Steininger committed
31 32 33 34 35
        # 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
36

37 38
        weight = observable.domain[1].weight(1)

Theo Steininger's avatar
Theo Steininger committed
39
        k = observable.shape[0]
40
        n = observable.shape[1]
Theo Steininger's avatar
Theo Steininger committed
41 42

        obs_val = observable.val.get_full_data()
43
        obs_mean = observable.ensemble_mean().val.get_full_data()
Theo Steininger's avatar
Theo Steininger committed
44

45 46 47 48 49 50
        # divide out profile
        obs_val /= profile
        obs_mean /= profile
        measured_data = measured_data / profile

        u_val = obs_val - obs_mean
51 52

        # compute quantities for OAS estimator
53 54 55
        mu = np.vdot(u_val, u_val)*weight/k
        self.logger.debug("mu: %f" % mu)

Theo Steininger's avatar
Theo Steininger committed
56
        alpha = (np.einsum(u_val, [0, 1], u_val, [2, 1])**2).sum()
57 58 59
        # correct the volume factor: one factor comes from the internal scalar
        # product and one from the trace
        alpha *= weight**2
60

61 62
        numerator = (1 - 2./n)*alpha + mu**2
        denominator = (k + 1 - 2./n) * (alpha - (mu**2)/n)
63 64 65 66

        if denominator == 0:
            rho = 1
        else:
67
            rho = np.min([1, numerator/denominator])
68

Theo Steininger's avatar
Theo Steininger committed
69 70
        self.logger.debug("rho: %f = %f / %f" % (rho, numerator, denominator))

71 72
        # rescale U half/half
        u_val *= np.sqrt(1-rho)
Theo Steininger's avatar
Theo Steininger committed
73 74
        U = observable.copy_empty()
        U.val = u_val
75 76 77 78 79 80 81

        # we assume that data_covariance_operator is a DiagonalOperator
        if not isinstance(data_covariance_operator, DiagonalOperator):
            raise TypeError("data_covariance_operator must be a NIFTY "
                            "DiagonalOperator.")

        A_bare_diagonal = data_covariance_operator.diagonal(bare=True)
82
        A_bare_diagonal.val += rho*mu/n
83 84 85 86 87 88
        A = DiagonalOperator(
                    domain=data_covariance_operator.domain,
                    diagonal=A_bare_diagonal,
                    bare=True, copy=False,
                    default_spaces=data_covariance_operator.default_spaces)

Theo Steininger's avatar
Theo Steininger committed
89 90 91 92 93 94
        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],
95
                            a_u_val, [2, 1])*weight)
Theo Steininger's avatar
Theo Steininger committed
96
        middle = np.linalg.inv(middle)
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
#        result_array = np.zeros(k)
#        for i in xrange(k):
#           c = measured_data - obs_val[i]
        c = measured_data - obs_mean

        # 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()

        # Pure NIFTy is
        # u_a_c = c.dot(a_u, spaces=1)
        # u_a_c_val = u_a_c.val.get_full_data()
        c_weighted_val = c.weight().val.get_full_data()
        u_a_c_val = np.einsum(c_weighted_val, [1], a_u_val, [0, 1])

        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_val *= -1
        second_summand = first_summand.copy_empty()
        second_summand.val = second_summand_val

124 125
        result_1 = -c.vdot(first_summand)
        result_2 = -c.vdot(second_summand)
126
        result = result_1 + result_2
Theo Steininger's avatar
Theo Steininger committed
127 128
        self.logger.info("Calculated (%s): %f + %f = %f" %
                         (self.observable_name, result_1, result_2, result))
129 130
#        result_array[i] = result
#        total_result = result_array.mean()
131 132

        return result