Merge branch 'simple_likelihood_test' into 'master'

Simple likelihood test

See merge request !2

imagine/__init__.py

@@ -6,6 +6,8 @@ from magnetic_fields import *
 from observers import *
 from priors import *
 
+from pymultinest import pymultinest
+
 from pipeline import Pipeline
 
 import nifty

imagine/carrier_mapper.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 
-def carrier_mapper(x, a=-np.inf, m=0, b=np.inf):
+def infinity_mapper(x, a=-np.inf, m=0, b=np.inf):
     """
     Maps x from [-inf, inf] into the interval [a, b], where x=0 -> m
     """
@@ -31,3 +31,11 @@ def carrier_mapper(x, a=-np.inf, m=0, b=np.inf):
     # strech y to the interval [a,b]
     y = y*(b-a) + a
     return y
+
+
+def unity_mapper(x, a=0, b=1):
+    """
+    Maps x from [0, 1] into the interval [a, b]
+    """
+    # rescale and shift
+    return x * (b-a) + a

imagine/likelihoods/__init__.py

 # -*- coding: utf-8 -*-
 
 from likelihood import Likelihood
+from ensemble_likelihood import EnsembleLikelihood

imagine/likelihoods/ensemble_likelihood/__init__.py

+
+from ensemble_likelihood import EnsembleLikelihood

imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py

 # -*- coding: utf-8 -*-
 
+import numpy as np
+
 from imagine.likelihoods.likelihood import Likelihood
 
+
 class EnsembleLikelihood(Likelihood):
-    def __init__(self, measured_data, data_covariance):
+    def __init__(self, measured_data, data_covariance_operator):
         self.measured_data = measured_data
-        self.data_covariance = data_covariance
+        self.data_covariance_operator = data_covariance_operator
 
     def __call__(self, observable):
-        mean = observable.mean(spaces=0)
+        # 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
+
+        k = observable.shape[0]
+
+        A = self.data_covariance_operator
+        obs_val = observable.val.get_full_data()
+        obs_mean = observable.mean(spaces=0).val.get_full_data()
+
+        u_val = observable.val.get_full_data() - obs_mean
+        U = observable.copy_empty()
+        U.val = u_val
+        a_u = A.inverse_times(U, spaces=1)
+
+        # 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]))
+        middle = np.linalg.inv(middle)
+        result_array = np.zeros(k)
+        for i in xrange(k):
+            c = self.measured_data - obs_val[i]
+
+            # assuming that A == A^dagger, this can be shortend
+            # a_c = A.inverse_times(c)
+            # u_a_c = a_c.dot(U, spaces=1)
+            # u_a_c = u_a_c.conjugate()
+
+            # and: double conjugate shouldn't make a difference
+            # u_a_c = c.conjugate().dot(a_u, spaces=1).conjugate()
+            u_a_c = c.dot(a_u, spaces=1)
+            u_a_c_val = u_a_c.val.get_full_data()
+
+            first_summand = A.inverse_times(c)
+
+            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 = first_summand.copy_empty()
+            second_summand.val = second_summand_val
+
+            result = c.dot(first_summand - second_summand)
+            result_array[i] = result
+
+        return -result_array.mean()

imagine/likelihoods/simple_likelihood/__init__.py

+# -*- coding: utf-8 -*-
+

imagine/likelihoods/simple_likelihood/simple_likelihood.py

+# -*- coding: utf-8 -*-
+
+import numpy as np
+
+from imagine.likelihoods.likelihood import Likelihood
+
+
+class SimpleLikelihood(Likelihood):
+    def __init__(self, measured_data):
+        self.measured_data = measured_data
+
+    def __call__(self, observable):
+        shape = observable.shape
+        data = self.measured_data.val.get_full_data()
+        obs = observable.val.get_full_data()
+
+        quadratic_diff = ((data - obs).conjugate() * (data - obs)).sum()
+        quadratic_diff /= np.prod(shape)
+
+        return -quadratic_diff

imagine/magnetic_fields/magnetic_field/magnetic_field_factory.py

@@ -4,7 +4,7 @@ import numpy as np
 
 from keepers import Loggable
 
-from imagine.carrier_mapper import carrier_mapper
+from imagine.carrier_mapper import unity_mapper
 
 from nifty import FieldArray, RGSpace
 
@@ -96,10 +96,13 @@ class MagneticFieldFactory(Loggable, object):
         for variable_name in variables:
             if variable_name in self.variable_to_parameter_mappings:
                 mapping = self.variable_to_parameter_mappings[variable_name]
-                mapped_variable = carrier_mapper(variables[variable_name],
-                                                 a=mapping[0],
-                                                 m=mapping[1],
-                                                 b=mapping[2])
+                mapped_variable = unity_mapper(variables[variable_name],
+                                               a=mapping[0],
+                                               b=mapping[2])
+#                mapped_variable = carrier_mapper(variables[variable_name],
+#                                                 a=mapping[0],
+#                                                 m=mapping[1],
+#                                                 b=mapping[2])
             else:
                 mapped_variable = np.float(variables[variable_name])
             parameter_dict[variable_name] = mapped_variable

imagine/priors/__init__.py

 # -*- coding: utf-8 -*-
 
 from prior import Prior
+from flat_prior import FlatPrior

imagine/priors/flat_prior/__init__.py

+# -*- coding: utf-8 -*-
+
+from flat_prior import FlatPrior

imagine/priors/flat_prior/flat_prior.py

+# -*- coding: utf-8 -*-
+
+from imagine.priors.prior import Prior
+
+
+class FlatPrior(Prior):
+    def __call__(self, cube, ndim, nparams):
+        return cube

imagine/priors/prior/__init__.py

+# -*- coding: utf-8 -*-
+
+from prior import Prior

imagine/priors/prior.py → imagine/priors/prior/prior.py

@@ -7,5 +7,5 @@ from keepers import Loggable
 
 class Prior(Loggable, object):
     @abc.abstractmethod
-    def __call__(self, parameters):
+    def __call__(self, cube, ndim, nparams):
         raise NotImplemented

imagine/pymultinest/__init__.py

+# -*- coding: utf-8 -*-
+
+from pymultinest import pymultinest

imagine/pymultinest/mpi4py_shadow.py

+# -*- coding: utf-8 -*-
+
+
+class _MPI(object):
+    def __init__(self):
+        self.COMM_WORLD = _COMM_WORLD()
+
+
+class _COMM_WORLD(object):
+    def Get_size(self):
+        return 1
+
+MPI = _MPI()

imagine/pymultinest/pymultinest.py

+# -*- coding: utf-8 -*-
+
+import sys
+import mpi4py
+import mpi4py_shadow
+
+mpi4py_backup = sys.modules['mpi4py']
+sys.modules['mpi4py'] = mpi4py_shadow
+
+import pymultinest
+
+sys.modules['mpi4py'] = mpi4py_backup