Commit e3707e06 authored by Theo Steininger's avatar Theo Steininger

Merge branch 'advanced_ensemble_likelihood' into 'master'

Advanced ensemble likelihood

See merge request !3
parents 398c77f2 a16031fa
......@@ -3,7 +3,7 @@ FROM ubuntu:latest
RUN apt-get update
RUN apt-get install -y build-essential python python-pip python-dev git gfortran autoconf gsl-bin libgsl-dev wget unzip
RUN pip install numpy scipy cython astropy ipython
RUN pip install numpy scipy cython astropy ipython==5.3.0
RUN mkdir /home/Downloads
WORKDIR /home/Downloads
......
......@@ -50,7 +50,7 @@ CXX = g++
CXXFLAGS = -m64 -O2 -fopenmp -Wall
#
# With OSX g++-mp-4.3
#CXXFLAGS = -fopenmp -O2 -fno-inline-functions -Wall -Wextra -Wno-unknown-pragmas
#CXXFLAGS = -fopenmp -O2 -g -fno-inline-functions -Wall -Wextra -Wno-unknown-pragmas -ansi
#
# Which Fortran compiler you are using? Only if compiling with NE2001 and/or Galprop
......@@ -203,5 +203,8 @@ tarfile:
tar cvzf hammurabi.tgz *cpp *h Makefile *.f README hampy
unittar:
tar cvzf hammurabi_unit_test.tgz unit_test/inputs unit_test/GALDEF unit_test/FITS unit_test/negrid_n400.bin unit_test/ref.mini unit_test/ref unit_test/ref.big
tar cvzf hammurabi_unit_test_inputs.tgz unit_test/inputs unit_test/GALDEF unit_test/FITS unit_test/negrid_n400.bin
tar cvzf hammurabi_unit_test_ref.tgz unit_test/ref unit_test/ref.big
tar cvzf hammurabi_unit_test_ref.mini.tgz unit_test/ref.mini
tar cvzf hammurabi_unit_test_ref.big.tgz unit_test/ref.big
......@@ -2,6 +2,8 @@
import numpy as np
from nifty import DiagonalOperator
from imagine.likelihoods.likelihood import Likelihood
from imagine.create_ring_profile import create_ring_profile
......@@ -24,18 +26,17 @@ class EnsembleLikelihood(Likelihood):
self.data_covariance_operator,
self.profile)
def _process_simple_field(self, field, measured_data,
def _process_simple_field(self, observable, measured_data,
data_covariance_operator, profile):
# 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
observable = field
k = observable.shape[0]
n = observable.shape[1]
A = data_covariance_operator
obs_val = observable.val.get_full_data()
obs_mean = observable.ensemble_mean().get_full_data()
......@@ -45,8 +46,37 @@ class EnsembleLikelihood(Likelihood):
measured_data = measured_data / profile
u_val = obs_val - obs_mean
# compute quantities for OAS estimator
mu = np.vdot(u_val, u_val)/n
alpha = (np.einsum(u_val, [0, 1], u_val, [2, 0])**2).sum()
numerator = alpha + mu**2
denominator = (k + 1) / (alpha - (mu**2)/n)
if denominator == 0:
rho = 1
else:
rho = np.min(1, numerator/denominator)
# rescale U half/half
u_val *= np.sqrt(1-rho)
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)
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_u = A.inverse_times(U, spaces=1)
# build middle-matrix (kxk)
......
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