Commit d5859a78 authored by Theo Steininger's avatar Theo Steininger

Refactored EnsembleLikelihood and adapted ParameterFile

parent fd6d782c
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
import numpy as np import numpy as np
from nifty import DiagonalOperator, FieldArray, Field from nifty import FieldArray, Field
from imagine.likelihoods.likelihood import Likelihood from imagine.likelihoods.likelihood import Likelihood
...@@ -30,7 +30,7 @@ class EnsembleLikelihood(Likelihood): ...@@ -30,7 +30,7 @@ class EnsembleLikelihood(Likelihood):
self.data_covariance_operator) self.data_covariance_operator)
def _process_simple_field(self, observable, measured_data, def _process_simple_field(self, observable, measured_data,
data_covariance_operator): data_covariance):
# https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization
# B = A^{-1} + U U^dagger # B = A^{-1} + U U^dagger
# A = data_covariance # A = data_covariance
...@@ -50,7 +50,7 @@ class EnsembleLikelihood(Likelihood): ...@@ -50,7 +50,7 @@ class EnsembleLikelihood(Likelihood):
self.logger.debug("mu: %f" % mu) self.logger.debug("mu: %f" % mu)
alpha = (np.einsum(u_val, [0, 1], u_val, [2, 1])**2).sum() alpha = (np.einsum(u_val, [0, 1], u_val, [2, 1])**2).sum()
alpha /= k*2 alpha /= k**2
numerator = (1 - 2./n)*alpha + (mu*n)**2 numerator = (1 - 2./n)*alpha + (mu*n)**2
denominator = (k + 1 - 2./n) * (alpha - ((mu*n)**2)/n) denominator = (k + 1 - 2./n) * (alpha - ((mu*n)**2)/n)
...@@ -65,12 +65,7 @@ class EnsembleLikelihood(Likelihood): ...@@ -65,12 +65,7 @@ class EnsembleLikelihood(Likelihood):
# rescale U half/half # rescale U half/half
u_val *= np.sqrt(1-rho) / np.sqrt(k) u_val *= np.sqrt(1-rho) / np.sqrt(k)
# we assume that data_covariance_operator is a DiagonalOperator A_diagonal_val = data_covariance
if not isinstance(data_covariance_operator, DiagonalOperator):
raise TypeError("data_covariance_operator must be a NIFTY "
"DiagonalOperator.")
A_diagonal_val = data_covariance_operator.diagonal(bare=False).val
self.logger.info(('rho*mu', rho*mu, self.logger.info(('rho*mu', rho*mu,
'rho', rho, 'rho', rho,
'mu', mu, 'mu', mu,
......
...@@ -6,7 +6,6 @@ ...@@ -6,7 +6,6 @@
<root> <root>
<!-- observable output --> <!-- observable output -->
<Obsout> <Obsout>
</Obsout> </Obsout>
<!-- physical field grid in/out --> <!-- physical field grid in/out -->
<!-- resolution defined in './Grid/Box' --> <!-- resolution defined in './Grid/Box' -->
...@@ -74,7 +73,7 @@ ...@@ -74,7 +73,7 @@
<!-- magnetic fields --> <!-- magnetic fields -->
<MagneticField> <MagneticField>
<!-- regular fields --> <!-- regular fields -->
<Regular type="WMAP"> <Regular type="Verify">
<!-- WMAP LSA --> <!-- WMAP LSA -->
<WMAP> <WMAP>
<b0 value="1.2"/> <!-- microG --> <b0 value="1.2"/> <!-- microG -->
...@@ -266,9 +265,11 @@ ...@@ -266,9 +265,11 @@
<alpha value="3.0"/> <alpha value="3.0"/>
<beta value="0.0"/> <beta value="0.0"/>
<theta value="0.0"/> <theta value="0.0"/>
<hr value="5.0"/> <!-- kpc --> <r0 value="5.0"/> <!-- kpc -->
<hz value="1.0"/> <!-- kpc --> <z0 value="1.0"/> <!-- kpc -->
<je value="0.25"/> <!-- local CRE flux norm factor @ 10GeV --> <!-- by default, we choose AMS02 20.6GeV -->
<E0 value="20.6"/> <!-- CRE kinetic energy reference, GeV -->
<j0 value="0.0217"/> <!-- local CRE flux norm factor @ E0 -->
</Analytic> </Analytic>
<Numeric type="2D"> <Numeric type="2D">
...@@ -297,7 +298,9 @@ ...@@ -297,7 +298,9 @@
<Verify> <Verify>
<alpha value="3.0"/> <alpha value="3.0"/>
<r0 value="3.0"/> <!-- cutoff radius, kpc --> <r0 value="3.0"/> <!-- cutoff radius, kpc -->
<je value="0.25"/> <!-- local CRE flux norm factor @ 10GeV --> <!-- by default, we choose AMS02 20.6GeV -->
<E0 value="20.6"/> <!-- CRE kinetic energy reference, GeV -->
<j0 value="0.0217"/> <!-- local CRE flux norm factor @ E0 -->
</Verify> </Verify>
</CRE> </CRE>
</root> </root>
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