From d5859a78d5c648fc397ba3559f300e61d29a4b82 Mon Sep 17 00:00:00 2001
From: Theo Steininger <theo.steininger@ultimanet.de>
Date: Fri, 8 Dec 2017 22:49:12 +0100
Subject: [PATCH] Refactored EnsembleLikelihood and adapted ParameterFile

---
 .../ensemble_likelihood/ensemble_likelihood.py    | 13 ++++---------
 .../hammurapy/input/default_parameters.xml        | 15 +++++++++------
 2 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py b/imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py
index 88af268..28f8483 100644
--- a/imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py
+++ b/imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from nifty import DiagonalOperator, FieldArray, Field
+from nifty import FieldArray, Field
 
 from imagine.likelihoods.likelihood import Likelihood
 
@@ -30,7 +30,7 @@ class EnsembleLikelihood(Likelihood):
                                           self.data_covariance_operator)
 
     def _process_simple_field(self, observable, measured_data,
-                              data_covariance_operator):
+                              data_covariance):
         # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization
         # B = A^{-1} + U U^dagger
         # A = data_covariance
@@ -50,7 +50,7 @@ class EnsembleLikelihood(Likelihood):
         self.logger.debug("mu: %f" % mu)
 
         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
         denominator = (k + 1 - 2./n) * (alpha - ((mu*n)**2)/n)
@@ -65,12 +65,7 @@ class EnsembleLikelihood(Likelihood):
         # rescale U half/half
         u_val *= np.sqrt(1-rho) / np.sqrt(k)
 
-        # 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_diagonal_val = data_covariance_operator.diagonal(bare=False).val
+        A_diagonal_val = data_covariance
         self.logger.info(('rho*mu', rho*mu,
                           'rho', rho,
                           'mu', mu,
diff --git a/imagine/observers/hammurapy/input/default_parameters.xml b/imagine/observers/hammurapy/input/default_parameters.xml
index bb8b11a..b6fc903 100644
--- a/imagine/observers/hammurapy/input/default_parameters.xml
+++ b/imagine/observers/hammurapy/input/default_parameters.xml
@@ -6,7 +6,6 @@
 <root>
     <!-- observable output -->
     <Obsout>
-
     </Obsout>
     <!-- physical field grid in/out -->
     <!-- resolution defined in './Grid/Box' -->
@@ -74,7 +73,7 @@
     <!-- magnetic fields -->
     <MagneticField>
         <!-- regular fields -->
-        <Regular type="WMAP">
+        <Regular type="Verify">
             <!-- WMAP LSA -->
             <WMAP>
                 <b0 value="1.2"/> <!-- microG -->
@@ -266,9 +265,11 @@
             <alpha value="3.0"/>
             <beta value="0.0"/>
             <theta value="0.0"/>
-            <hr value="5.0"/> <!-- kpc -->
-            <hz value="1.0"/> <!-- kpc -->
-            <je value="0.25"/> <!-- local CRE flux norm factor @ 10GeV -->
+            <r0 value="5.0"/> <!-- kpc -->
+            <z0 value="1.0"/> <!-- kpc -->
+            <!-- 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>
         
         <Numeric type="2D">
@@ -297,7 +298,9 @@
         <Verify>
             <alpha value="3.0"/>
             <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>
     </CRE>
 </root>
-- 
GitLab