Removed volume factors from EnsembleLikelihood

parent dfb2cb45
......@@ -5,7 +5,6 @@ import numpy as np
from nifty import DiagonalOperator, FieldArray, Field
from imagine.likelihoods.likelihood import Likelihood
from imagine.create_ring_profile import create_ring_profile
class EnsembleLikelihood(Likelihood):
......@@ -14,12 +13,6 @@ class EnsembleLikelihood(Likelihood):
self.observable_name = observable_name
self.measured_data = self._strip_data(measured_data)
self.data_covariance_operator = data_covariance_operator
self.data_covariance_includes_profile = False
if profile is None:
profile = create_ring_profile(
self.measured_data.val.get_full_data())
self.profile = profile
def _strip_data(self, data):
# if the first element in the domain tuple is a FieldArray we must
......@@ -34,41 +27,30 @@ class EnsembleLikelihood(Likelihood):
field = observable[self.observable_name]
return self._process_simple_field(field,
self.measured_data,
self.data_covariance_operator,
self.profile)
self.data_covariance_operator)
def _process_simple_field(self, observable, measured_data,
data_covariance_operator, profile):
data_covariance_operator):
# 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
weight = observable.domain[1].weight(1)
k = observable.shape[0]
n = observable.shape[1]
obs_val = observable.val.get_full_data()
obs_mean = observable.ensemble_mean().val.get_full_data()
# divide out profile
obs_val /= profile
obs_mean /= profile
measured_data = measured_data / profile
u_val = obs_val - obs_mean
# compute quantities for OAS estimator
mu = np.vdot(u_val, u_val)*weight/k/n
mu = np.vdot(u_val, u_val)/k/n
self.logger.debug("mu: %f" % mu)
alpha = (np.einsum(u_val, [0, 1], u_val, [2, 1])**2).sum()
alpha /= k*2
# correct the volume factor: one factor comes from the internal scalar
# product and one from the trace
alpha *= weight**2
numerator = (1 - 2./n)*alpha + (mu*n)**2
denominator = (k + 1 - 2./n) * (alpha - ((mu*n)**2)/n)
......@@ -82,38 +64,31 @@ class EnsembleLikelihood(Likelihood):
# rescale U half/half
u_val *= np.sqrt(1-rho) / np.sqrt(k)
U = observable.copy_empty()
U.val = u_val
# 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)
if not self.data_covariance_includes_profile:
A_bare_diagonal *= (profile**2)
A_bare_diagonal.val += rho*mu
A = DiagonalOperator(
domain=data_covariance_operator.domain,
diagonal=A_bare_diagonal,
bare=True, copy=False,
default_spaces=data_covariance_operator.default_spaces)
A_diagonal_val = data_covariance_operator.diagonal(bare=False).val
self.logger.info(('rho*mu', rho*mu,
'rho', rho,
'mu', mu,
'alhpa', alpha))
A_diagonal_val += rho*mu
a_u = A.inverse_times(U, spaces=1)
a_u_val = u_val/A_diagonal_val
# 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])*weight)
a_u_val, [2, 1]))
middle = np.linalg.inv(middle)
c = measured_data - obs_mean
# If the data was incomplete, i.e. contains np.NANs, set those values
# to zero.
c.val.data = np.nan_to_num(c.val.data)
# assuming that A == A^dagger, this can be shortend
# a_c = A.inverse_times(c)
# u_a_c = a_c.dot(U, spaces=1)
......@@ -125,20 +100,21 @@ class EnsembleLikelihood(Likelihood):
# 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])
c_val = c.val.get_full_data()
u_a_c_val = np.einsum(c_val, [1], a_u_val, [0, 1])
first_summand = A.inverse_times(c)
first_summand_val = c_val/A_diagonal_val
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
# # second_summand_val *= -1
# second_summand = first_summand.copy_empty()
# second_summand.val = second_summand_val
result_1 = c.vdot(first_summand)
result_2 = -c.vdot(second_summand)
result_1 = np.vdot(c_val, first_summand_val)
result_2 = -np.vdot(c_val, second_summand_val)
result = -(result_1 + result_2)
print (result_1, result_2, result)
self.logger.info("Calculated (%s): -(%f + %f) = %f" %
(self.observable_name, result_1, result_2, result))
# result_array[i] = result
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment