diff --git a/imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py b/imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py index 88af268732b455bbe80670218f1b8ffd26c99c45..28f84839207e3c6a9314f204635ef819aa99e921 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 bb8b11a5c273458d8a935502b6f24f9a6749c4b7..b6fc9033e12bc0d930b089fddbac11e9e60d2c09 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>