Refactored structure of observable_mixins. Ensemble_mean is now calculated individually.

Add first code snippets for galmag model.

imagine/__init__.py

@@ -3,6 +3,7 @@ from .version import __version__
 
 from likelihoods import *
 from magnetic_fields import *
+from observables import *
 from observers import *
 from priors import *
 

imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py

@@ -37,7 +37,7 @@ class EnsembleLikelihood(Likelihood):
 
         A = data_covariance_operator
         obs_val = observable.val.get_full_data()
-        obs_mean = observable.mean(spaces=0).val.get_full_data()
+        obs_mean = observable.ensemble_mean().get_full_data()
 
         # divide out profile
         obs_val /= profile
@@ -87,8 +87,8 @@ class EnsembleLikelihood(Likelihood):
         result_1 = -c.dot(first_summand)
         result_2 = -c.dot(second_summand)
         result = result_1 + result_2
-        self.logger.debug("Calculated %i of %i: %f + %f = %f" %
-                          (i, k, result_1, result_2, result))
+        self.logger.debug("Calculated: %f + %f = %f" %
+                          (result_1, result_2, result))
 #        result_array[i] = result
 #        total_result = result_array.mean()
         total_result = result

imagine/magnetic_fields/galmag_field/__init__.py

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

imagine/magnetic_fields/galmag_field/galmag_field.py

+# -*- coding: utf-8 -*-
+
+from imagine.magnetic_fields.magnetic_field import MagneticField
+
+try:
+    import galmag
+except ImportError:
+    galmag_available = False
+else:
+    galmag_available = True
+
+class GalMagField(MagneticField):
+    def __init__(self, parameters=[], domain=None, val=None, dtype=None,
+                 distribution_strategy=None, copy=False, random_seed=None):
+        if not galmag_available:
+            raise ImportError("GalMag module needed but not available.")
+
+        super(MagneticField, self).__init__(
+                                parameters=parameters,
+                                domain=domain,
+                                val=val,
+                                dtype=dtype,
+                                distribution_strategy=distribution_strategy,
+                                copy=copy,
+                                random_seed=random_seed)
+
+    @property
+    def parameter_list(self):
+        parameter_list = ['b_x', 'b_y', 'b_z']
+        return parameter_list
+
+    def _create_field(self):
+        grid_domain = self.domain[1]
+        box_resolution = grid_domain.shape
+        bl = np.array(grid_domain.distances) * box_resolution
+        box_limits = [[-bl[0], bl[0]],
+                      [-bl[1], bl[1]],
+                      [-bl[2], bl[2]]]
+
+        ensemble_size = self.shape[0]
+        val = self.cast(None)
+        for i in ensemble_size:
+            b = galmag.B_field.B_field(box_limits, box_resolution)
+            val[i, :, :, :, 0] = b.x
+            val[i, :, :, :, 1] = b.y
+            val[i, :, :, :, 2] = b.z
+
+        return val

imagine/magnetic_fields/galmag_field/galmag_field_factory.py

+# -*- coding: utf-8 -*-
+
+from imagine.magnetic_fields.magnetic_field.magnetic_field_factory \
+    import MagneticFieldFactory
+
+from galmag_field import GalMagField
+
+
+class GalMagFieldFactory(MagneticFieldFactory):
+    @property
+    def magnetic_field_class(self):
+        return GalMagField
+
+    @property
+    def _initial_parameter_defaults(self):
+        defaults = {'b_x': 0,
+                    'b_y': 0,
+                    'b_z': 0}
+        return defaults
+
+    @property
+    def _initial_variable_to_parameter_mappings(self):
+        return {'b_x': [-100, 0, 100],
+                'b_y': [-100, 0, 100],
+                'b_z': [-100, 0, 100]}

imagine/observables/__init__.py

+# -*- coding: utf-8 -*-
+
+from observable import Observable

imagine/observables/observable.py

+# -*- coding: utf-8 -*-
+
+
+from nifty import Field, FieldArray
+
+
+class Observable(Field):
+    def __init__(self, domain=None, val=None, dtype=None,
+                 distribution_strategy=None, copy=False):
+
+        super(Observable, self).__init__(
+                                domain=domain,
+                                val=val,
+                                dtype=dtype,
+                                distribution_strategy=distribution_strategy,
+                                copy=copy)
+
+        assert(len(self.domain) == 2)
+        assert(isinstance(self.domain[0], FieldArray))
+
+    def ensemble_mean(self):
+        try:
+            self._ensemble_mean
+        except(AttributeError):
+            self._ensemble_mean = self.mean(spaces=0)
+        finally:
+            return self._ensemble_mean

imagine/observers/hammurapy/hammurapy.py

@@ -30,7 +30,8 @@ class Hammurapy(Observer):
         dust_template_fname = os.path.join(self.conf_directory,
                                            'IQU_dust.fits')
 
-        self.basic_parameters = {'obs_shell_index_numb': '1',
+        self.basic_parameters = {'B_ran_mem_lim': '4',
+                                 'obs_shell_index_numb': '1',
                                  'total_shell_numb': '3',
                                  'vec_size_R': '500',
                                  'max_radius': '35',
@@ -106,7 +107,9 @@ class Hammurapy(Observer):
         grid_space = magnetic_field.domain[1]
         lx, ly, lz = np.array(grid_space.shape)*np.array(grid_space.distances)
         nx, ny, nz = grid_space.shape
-        random_seed = magnetic_field.random_seed[local_ensemble_index]
+        if local_ensemble_index is not None:
+            random_seed = magnetic_field.random_seed[local_ensemble_index]
+            parameter_dict.update({'B_field_seed': random_seed})
 
         parameter_dict.update({'B_field_lx': lx,
                                'B_field_ly': ly,
@@ -116,7 +119,6 @@ class Hammurapy(Observer):
                                'B_field_nz': nz,
                                })
         parameter_dict.update({'obs_NSIDE': self.nside})
-        parameter_dict.update({'B_field_seed': random_seed})
 
     def _write_parameter_dict(self, parameter_dict, working_directory):
         parameters_string = ''
@@ -132,6 +134,9 @@ class Hammurapy(Observer):
                               ensemble_index):
         return observable_dict
 
+    def _add_ensemble_mean(self, observable_dict, working_directory):
+        pass
+
     def __call__(self, magnetic_field):
 
         if not isinstance(magnetic_field, self.magnetic_field_class):
@@ -158,6 +163,11 @@ class Hammurapy(Observer):
 
             # create dictionary for parameter file
             parameter_dict = self.basic_parameters.copy()
+
+            # set the parameters for a numerical run
+            parameter_dict['B_field_interp'] = 'T'
+            parameter_dict['use_B_analytic'] = 'F'
+
             self._build_parameter_dict(
                                     parameter_dict=parameter_dict,
                                     magnetic_field=magnetic_field,
@@ -182,4 +192,44 @@ class Hammurapy(Observer):
             finally:
                 self._remove_folder(working_directory)
 
+        # compute ensemble_mean on the last node and add it to the observables
+        rank = dummy.comm.rank
+        size = dummy.comm.size
+
+        if rank + 1 == size:
+            self.logger.debug("Processing ensemble mean.")
+            # create a temporary folder
+            working_directory = self._make_temp_folder()
+
+            # create dictionary for parameter file
+            parameter_dict = self.basic_parameters.copy()
+
+            # set the parameters for an analytic run
+            parameter_dict['use_B_analytic'] = 'T'
+
+            self._build_parameter_dict(
+                                    parameter_dict=parameter_dict,
+                                    magnetic_field=magnetic_field,
+                                    working_directory=working_directory,
+                                    local_ensemble_index=None)
+
+            self._write_parameter_dict(parameter_dict=parameter_dict,
+                                       working_directory=working_directory)
+
+            # call hammurabi
+            self._call_hammurabi(working_directory)
+
+            # if hammurabi failed, _add_ensemble_mean will fail
+            try:
+                self._add_ensemble_mean(observable_dict, working_directory)
+            except:
+                self.logger.critical("Hammurabi failed! Last call log:\n" +
+                                     self.last_call_log)
+                raise
+            finally:
+                self._remove_folder(working_directory)
+        else:
+            working_directory = None
+            self._add_ensemble_mean(observable_dict, working_directory)
+
         return observable_dict

imagine/observers/hammurapy/model_mixins/galmag_mixin.py

+# -*- coding: utf-8 -*-
+
+import os
+
+from imagine.magnetic_fields.galmag_field import GalMagField
+
+
+class GalMagMixin(object):
+    def __init__(self, hammurabi_executable, conf_directory='./confs',
+                 working_directory_base='.', nside=128):
+        self.__parameter_dict = {'B_field_type': '6'}
+        super(GalMagMixin, self).__init__(hammurabi_executable,
+                                          conf_directory,
+                                          working_directory_base,
+                                          nside)
+
+    @property
+    def magnetic_field_class(self):
+        return GalMagField
+
+    def _build_parameter_dict(self, parameter_dict, magnetic_field,
+                              working_directory, local_ensemble_index):
+        parameter_dict.update(self.__parameter_dict)
+
+        parameter_dict.update(magnetic_field.parameters)
+
+        # write the field array to disk
+        array_file_name = os.path.join(working_directory, 'b_field.arr')
+        field_array = magnetic_field.val.get_full_data()
+        # hammurabi reads the array like B[z,y,x]
+        # field_array = field_array.transpose((1, 2, 3, 0))
+        field_array.tofile(array_file_name)
+
+        super(GalMagMixin, self)._build_parameter_dict(parameter_dict,
+                                                       magnetic_field,
+                                                       working_directory,
+                                                       local_ensemble_index)

imagine/observers/hammurapy/model_mixins/jaffe13_mixin.py

@@ -10,9 +10,6 @@ class Jaffe13Mixin(object):
                                  'B_field_do_random': 'T',
                                  'B_analytic_beta': '0',
                                  'B_field_RMS_uG': '3.5',
-                                 'B_field_interp': 'T',
-                                 'use_B_analytic': 'F',
-                                 'B_ran_mem_lim': '4',
                                  'bb_molr_aniso': 'T',
                                  'bb_ord_interarm': 'T',
                                  'bb_scale_coh_amps': 'T',

imagine/observers/hammurapy/model_mixins/jf12_mixin.py

@@ -8,10 +8,7 @@ class JF12Mixin(object):
                  working_directory_base='.', nside=128):
         self.__parameter_dict = {'B_field_type': '7',
                                  'B_field_do_random': 'T',
-                                 'B_analytic_beta': '1.36',
-                                 'B_field_interp': 'T',
-                                 'use_B_analytic': 'F',
-                                 'B_ran_mem_lim': '4'}
+                                 'B_analytic_beta': '1.36'}
         super(JF12Mixin, self).__init__(hammurabi_executable,
                                         conf_directory,
                                         working_directory_base,

imagine/observers/hammurapy/observable_mixins/dm_mixin.py

 # -*- coding: utf-8 -*-
 
-import os
-
-from nifty import Field, HPSpace
-
-
-class DMMixin(object):
-    def __init__(self, hammurabi_executable, conf_directory='./confs',
-                 working_directory_base='.', nside=128):
-        self.__hpSpace = HPSpace(nside=int(nside))
-        super(DMMixin, self).__init__(hammurabi_executable,
-                                      conf_directory,
-                                      working_directory_base,
-                                      nside)
-
-    def _initialize_observable_dict(self, observable_dict, magnetic_field):
-        ensemble_space = magnetic_field.domain[0]
-
-        observable_dict['dm'] = Field(domain=(ensemble_space, self.__hpSpace),
-                                      distribution_strategy='equal')
-        super(DMMixin, self)._initialize_observable_dict(observable_dict,
-                                                         magnetic_field)
-
-    def _build_parameter_dict(self, parameter_dict, magnetic_field,
-                              working_directory, local_ensemble_index):
-        obs_DM_file_name = os.path.join(working_directory, 'dm.fits')
-        parameter_dict['do_dm'] = 'T'
-        parameter_dict['obs_DM_file_name'] = obs_DM_file_name
-        super(DMMixin, self)._build_parameter_dict(parameter_dict,
-                                                   magnetic_field,
-                                                   working_directory,
-                                                   local_ensemble_index)
-
-    def _fill_observable_dict(self, observable_dict, working_directory,
-                              local_ensemble_index):
-        self.logger.debug('Reading DM-map.')
-        [dm_map] = self._read_fits_file(path=working_directory,
-                                        name='dm.fits',
-                                        nside=self.nside)
-
-        dm_field = observable_dict['dm']
-        dm_field.val.data[local_ensemble_index] = dm_map
-
-        super(DMMixin, self)._fill_observable_dict(observable_dict,
-                                                   working_directory,
-                                                   local_ensemble_index)
+from mixin_base import MixinBase
+
+
+class DMMixin(MixinBase):
+
+    @property
+    def __config_dict(self):
+        return {'obs_name': 'dm',
+                'component_names': ['dm'],
+                'parameter_dict_update': {'do_dm': 'T'},
+                'filename_key': 'obs_DM_file_name',
+                }
+
+    def _initialize_observable_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._initialize_observable_dict_helper(*args, **fat_kwargs)
+        super(DMMixin, self)._initialize_observable_dict(*args, **kwargs)
+
+    def _build_parameter_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._build_parameter_dict_helper(*args, **fat_kwargs)
+        super(DMMixin, self)._build_parameter_dict(*args, **kwargs)
+
+    def _fill_observable_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._fill_observable_dict_helper(*args, **fat_kwargs)
+        super(DMMixin, self)._fill_observable_dict(*args, **kwargs)
+
+    def _add_ensemble_mean(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._add_ensemble_mean_helper(*args, **fat_kwargs)
+        super(DMMixin, self)._add_ensemble_mean(*args, **kwargs)

imagine/observers/hammurapy/observable_mixins/dust_mixin.py

 # -*- coding: utf-8 -*-
 
-import os
-
-from nifty import Field, HPSpace
-
-
-class DustMixin(object):
-    def __init__(self, hammurabi_executable, conf_directory='./confs',
-                 working_directory_base='.', nside=128):
-        self.__hpSpace = HPSpace(nside=int(nside))
-        super(DustMixin, self).__init__(hammurabi_executable,
-                                        conf_directory,
-                                        working_directory_base,
-                                        nside)
-
-    def _initialize_observable_dict(self, observable_dict, magnetic_field):
-        ensemble_space = magnetic_field.domain[0]
-
-        for name in ['dust_I', 'dust_Q', 'dust_U']:
-            observable_dict[name] = Field(domain=(ensemble_space,
-                                                  self.__hpSpace),
-                                          distribution_strategy='equal')
-
-        super(DustMixin, self)._initialize_observable_dict(observable_dict,
-                                                           magnetic_field)
-
-    def _build_parameter_dict(self, parameter_dict, magnetic_field,
-                              working_directory, local_ensemble_index):
-        obs_dust_file_name = os.path.join(working_directory, 'IQU_dust.fits')
-        parameter_dict['do_dust'] = 'T'
-        parameter_dict['obs_dust_file_name'] = obs_dust_file_name
-        super(DustMixin, self)._build_parameter_dict(parameter_dict,
-                                                     magnetic_field,
-                                                     working_directory,
-                                                     local_ensemble_index)
-
-    def _fill_observable_dict(self, observable_dict, working_directory,
-                              local_ensemble_index):
-        self.logger.debug('Reading Dust-map.')
-        [dust_I, dust_Q, dust_U] = self._read_fits_file(path=working_directory,
-                                                        name='IQU_dust.fits',
-                                                        nside=self.nside)
-
-        observable_dict['dust_I'].val.data[local_ensemble_index, :] = dust_I
-        observable_dict['dust_Q'].val.data[local_ensemble_index, :] = dust_Q
-        observable_dict['dust_U'].val.data[local_ensemble_index, :] = dust_U
-
-        super(DustMixin, self)._fill_observable_dict(observable_dict,
-                                                     working_directory,
-                                                     local_ensemble_index)
+from mixin_base import MixinBase
+
+
+class DustMixin(MixinBase):
+
+    @property
+    def __config_dict(self):
+        return {'obs_name': 'dust',
+                'component_names': ['dust_I', 'dust_Q', 'dust_U'],
+                'parameter_dict_update': {'do_dust': 'T'},
+                'filename_key': 'obs_dust_file_name',
+                }
+
+    def _initialize_observable_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._initialize_observable_dict_helper(*args, **fat_kwargs)
+        super(DustMixin, self)._initialize_observable_dict(*args, **kwargs)
+
+    def _build_parameter_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._build_parameter_dict_helper(*args, **fat_kwargs)
+        super(DustMixin, self)._build_parameter_dict(*args, **kwargs)
+
+    def _fill_observable_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._fill_observable_dict_helper(*args, **fat_kwargs)
+        super(DustMixin, self)._fill_observable_dict(*args, **kwargs)
+
+    def _add_ensemble_mean(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._add_ensemble_mean_helper(*args, **fat_kwargs)
+        super(DustMixin, self)._add_ensemble_mean(*args, **kwargs)

imagine/observers/hammurapy/observable_mixins/ff_mixin.py

 # -*- coding: utf-8 -*-
 
-import os
-
-from nifty import Field, HPSpace
-
-
-class FFMixin(object):
-    def __init__(self, hammurabi_executable, conf_directory='./confs',
-                 working_directory_base='.', nside=128):
-        self.__hpSpace = HPSpace(nside=int(nside))
-        super(FFMixin, self).__init__(hammurabi_executable,
-                                      conf_directory,
-                                      working_directory_base,
-                                      nside)
-    def _initialize_observable_dict(self, observable_dict, magnetic_field):
-        ensemble_space = magnetic_field.domain[0]
-
-        observable_dict['ff'] = Field(domain=(ensemble_space, self.__hpSpace),
-                                      distribution_strategy='equal')
-        super(FFMixin, self)._initialize_observable_dict(observable_dict,
-                                                         magnetic_field)
-
-    def _build_parameter_dict(self, parameter_dict, magnetic_field,
-                              working_directory, local_ensemble_index):
-        obs_ff_file_name = os.path.join(working_directory, 'free.fits')
-        parameter_dict['do_ff'] = 'T'
-        parameter_dict['obs_ff_file_name'] = obs_ff_file_name
-        super(FFMixin, self)._build_parameter_dict(parameter_dict,
-                                                   magnetic_field,
-                                                   working_directory,
-                                                   local_ensemble_index)
-
-    def _fill_observable_dict(self, observable_dict, working_directory,
-                              local_ensemble_index):
-        self.logger.debug('Reading FF-map.')
-        [ff_map] = self._read_fits_file(path=working_directory,
-                                        name='free.fits',
-                                        nside=self.nside)
-
-        ff_field = observable_dict['ff']
-        ff_field.val.data[local_ensemble_index] = ff_map
-
-        super(FFMixin, self)._fill_observable_dict(observable_dict,
-                                                   working_directory,
-                                                   local_ensemble_index)
+from mixin_base import MixinBase
+
+
+class FFMixin(MixinBase):
+
+    @property
+    def __config_dict(self):
+        return {'obs_name': 'ff',
+                'component_names': ['ff'],
+                'parameter_dict_update': {'do_ff': 'T'},
+                'filename_key': 'obs_ff_file_name',
+                }
+
+    def _initialize_observable_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._initialize_observable_dict_helper(*args, **fat_kwargs)
+        super(FFMixin, self)._initialize_observable_dict(*args, **kwargs)
+
+    def _build_parameter_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._build_parameter_dict_helper(*args, **fat_kwargs)
+        super(FFMixin, self)._build_parameter_dict(*args, **kwargs)
+
+    def _fill_observable_dict(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._fill_observable_dict_helper(*args, **fat_kwargs)
+        super(FFMixin, self)._fill_observable_dict(*args, **kwargs)
+
+    def _add_ensemble_mean(self, *args, **kwargs):
+        fat_kwargs = kwargs.copy()
+        fat_kwargs.update({'config_dict': self.__config_dict})
+        self._add_ensemble_mean_helper(*args, **fat_kwargs)
+        super(FFMixin, self)._add_ensemble_mean(*args, **kwargs)

imagine/observers/hammurapy/observable_mixins/mixin_base.py

+# -*- coding: utf-8 -*-
+
+import os
+
+import numpy as np
+
+from mpi4py import MPI
+
+from d2o import distributed_data_object
+
+from nifty import HPSpace
+
+from imagine.observables import Observable
+
+
+class MixinBase(object):
+    def __init__(self, hammurabi_executable, conf_directory='./confs',
+                 working_directory_base='.', nside=128):
+        self.__hpSpace = HPSpace(nside=int(nside))
+        super(MixinBase, self).__init__(hammurabi_executable,
+                                        conf_directory,
+                                        working_directory_base,
+                                        nside)
+
+    @property
+    def __config_dict(self):
+        return {'obs_name': '',
+                'component_names': [],
+                'parameter_dict_update': {},
+                'filename_key': '',