ensemble_likelihood.py 4.92 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
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
20 21
from nifty4 import Field
from .likelihood import Likelihood
22

23 24

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

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

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

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

54
        U = obs_val - obs_mean
55 56
        U *= np.sqrt(n)

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

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

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

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

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

81
        V_B = V/B
Theo Steininger's avatar
Theo Steininger committed
82 83 84

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

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

106 107 108
        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],
109
                                       second_summand_val, [0])
110 111 112
#        # second_summand_val *= -1
#        second_summand = first_summand.copy_empty()
#        second_summand.val = second_summand_val
113

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

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

125
        result = -0.5*(result_1 + result_2 + log_det)
126

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

        return result