ensemble_likelihood.py 4.93 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# IMAGINE is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
18

Theo Steininger's avatar
Theo Steininger committed
19 20
import numpy as np

21
from nifty import Field
22

23
from imagine.likelihoods.likelihood import Likelihood
24

25 26

class EnsembleLikelihood(Likelihood):
27
    def __init__(self, observable_name,  measured_data,
28
                 data_covariance, profile=None):
29
        self.observable_name = observable_name
30
        self.measured_data = self._strip_data(measured_data)
31 32 33
        if isinstance(data_covariance, Field):
            data_covariance = data_covariance.val.get_full_data()
        self.data_covariance = data_covariance
34
        self.use_determinant = True
35 36

    def __call__(self, observable):
37 38 39
        field = observable[self.observable_name]
        return self._process_simple_field(field,
                                          self.measured_data,
40
                                          self.data_covariance)
41

42
    def _process_simple_field(self, observable, measured_data,
43
                              data_covariance):
Theo Steininger's avatar
Theo Steininger committed
44 45 46 47 48
        # 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's avatar
bug fix  
Theo Steininger committed
49
        data_covariance = data_covariance.copy()
Theo Steininger's avatar
Theo Steininger committed
50
        k = observable.shape[0]
51
        n = observable.shape[1]
Theo Steininger's avatar
Theo Steininger committed
52 53

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

56
        U = obs_val - obs_mean
57 58
        U *= np.sqrt(n)

59
        # compute quantities for OAS estimator
60 61
        mu = np.vdot(U, U)/k/n
        alpha = (np.einsum(U, [0, 1], U, [2, 1])**2).sum()
62
        alpha /= k**2
63

64 65
        numerator = (1 - 2./n)*alpha + (mu*n)**2
        denominator = (k + 1 - 2./n) * (alpha - ((mu*n)**2)/n)
66 67 68 69

        if denominator == 0:
            rho = 1
        else:
70
            rho = np.min([1, numerator/denominator])
Theo Steininger's avatar
Theo Steininger committed
71 72
        self.logger.debug("rho: %f = %f / %f" % (rho, numerator, denominator))

73
        # rescale U half/half
74
        V = U * np.sqrt(1-rho) / np.sqrt(k)
75

76
        self.logger.info(('data_cov', np.mean(data_covariance),
Theo Steininger's avatar
Theo Steininger committed
77
                          'rho*mu', rho*mu,
78 79
                          'rho', rho,
                          'mu', mu,
80
                          'alpha', alpha))
81
        B = data_covariance + rho*mu
82

83
        V_B = V/B
Theo Steininger's avatar
Theo Steininger committed
84 85 86

        # build middle-matrix (kxk)
        middle = (np.eye(k) +
87 88
                  np.einsum(V.conjugate(), [0, 1],
                            V_B, [2, 1]))
Theo Steininger's avatar
Theo Steininger committed
89
        middle = np.linalg.inv(middle)
90 91
        c = measured_data - obs_mean

92 93
        # If the data was incomplete, i.e. contains np.NANs, set those values
        # to zero.
94
        c = np.nan_to_num(c)
95 96 97 98 99 100 101 102 103 104 105
        # 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()
106
        V_B_c = np.einsum(c, [1], V_B, [0, 1])
107

108 109 110
        first_summand_val = c/B
        second_summand_val = np.einsum(middle, [0, 1], V_B_c, [1])
        second_summand_val = np.einsum(V_B, [0, 1],
111
                                       second_summand_val, [0])
112 113 114
#        # second_summand_val *= -1
#        second_summand = first_summand.copy_empty()
#        second_summand.val = second_summand_val
115

116 117
        result_1 = np.vdot(c, first_summand_val)
        result_2 = -np.vdot(c, second_summand_val)
118 119 120

        # compute regularizing determinant of the covariance
        # det(A + UV^T) =  det(A) det(I + V^T A^-1 U)
121
        if self.use_determinant:
122 123
            log_det = np.sum(np.log(data_covariance +
                                    np.sum((obs_val-obs_mean)**2, axis=0)/k))/n
124
        else:
125
            log_det = 0.
126

127
        result = -0.5*(result_1 + result_2 + log_det)
128

129
        self.logger.info("Calculated (%s): -1/2(%g + %g + %g) = %g" %
130
                         (self.observable_name,
131
                          result_1, result_2, log_det, result))
132 133
#        result_array[i] = result
#        total_result = result_array.mean()
134 135

        return result