Skip to content
GitLab
Commit 6c8d5721 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'imagine_on_hamx' into 'master' 

Imagine on hamx

See merge request !4
parents 7d179bd2 1a2cb14e
Hide whitespace changes
Inline Side-by-side
imagine/likelihoods/ensemble_likelihood/ensemble_likelihood.py
View file @ 6c8d5721
......@@ -40,7 +40,7 @@ class EnsembleLikelihood(Likelihood):
n = observable.shape[1]
obs_val = observable.val.get_full_data()
obs_mean = observable.ensemble_mean().get_full_data()
obs_mean = observable.ensemble_mean().val.get_full_data()
# divide out profile
obs_val /= profile
......
imagine/magnetic_fields/__init__.py
View file @ 6c8d5721
......@@ -6,11 +6,5 @@ from magnetic_field import MagneticField, \
from constant_magnetic_field import ConstantMagneticField, \
ConstantMagneticFieldFactory
from jaffe13_magnetic_field import Jaffe13MagneticField, \
Jaffe13MagneticFieldFactory
from jf12_magnetic_field import JF12MagneticField, \
JF12MagneticFieldFactory
from wmap3yr_magnetic_field import WMAP3yrMagneticField, \
WMAP3yrMagneticFieldFactory
imagine/magnetic_fields/galmag_field/galmag_field.py deleted 100644 → 0
View file @ 7d179bd2
# -*- 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 deleted 100644 → 0
View file @ 7d179bd2
# -*- 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/magnetic_fields/jaffe13_magnetic_field/__init__.py deleted 100644 → 0
View file @ 7d179bd2
# -*- coding: utf-8 -*-
from jaffe13_magnetic_field import Jaffe13MagneticField
from jaffe13_magnetic_field_factory import Jaffe13MagneticFieldFactory
\ No newline at end of file
imagine/magnetic_fields/jaffe13_magnetic_field/jaffe13_magnetic_field.py deleted 100644 → 0
View file @ 7d179bd2
# -*- coding: utf-8 -*-
from imagine.magnetic_fields.magnetic_field import MagneticField
class Jaffe13MagneticField(MagneticField):
@property
def parameter_list(self):
parameter_list = ['B_f_ord', 'B_field_alpha', 'B_field_cutoff',
'B_ran_b2', 'B_ran_h_d', 'B_ran_h_d2', 'B_ran_h_r',
'B_ran_h_r2', 'bb_amps_0', 'bb_amps_1',
'bb_amps_2', 'bb_amps_3', 'bb_amps_4', 'bb_bar_a',
'bb_bar_boa', 'bb_bar_phi0_deg', 'bb_cr0_coh',
'bb_cr0_iso', 'bb_cr0_ord', 'bb_d0_iso', 'bb_d0',
'bb_delta_phi_iso_deg', 'bb_delta_phi_ord_deg',
'bb_disk_b0', 'bb_disk_h_d', 'bb_halo_b0',
'bb_halo_h_d', 'bb_phi0_deg', 'bb_pitch_biso',
'bb_pitch', 'bb_r_compconst', 'bb_r_innercut',
'bb_r_peak', 'bb_r_scale', 'bb_rmax_arms',
'bb_spiral_cpow', 'bb_spiral_h_d']
return parameter_list
def _create_field(self):
raise NotImplementedError
imagine/magnetic_fields/jaffe13_magnetic_field/jaffe13_magnetic_field_factory.py deleted 100644 → 0
View file @ 7d179bd2
# -*- coding: utf-8 -*-
from imagine.magnetic_fields.magnetic_field.magnetic_field_factory \
import MagneticFieldFactory
from jaffe13_magnetic_field import Jaffe13MagneticField
class Jaffe13MagneticFieldFactory(MagneticFieldFactory):
@property
def magnetic_field_class(self):
return Jaffe13MagneticField
@property
def _initial_parameter_defaults(self):
defaults = {'B_f_ord': 0.5,
'B_field_RMS_uG': 1.,
'B_field_alpha': -2.37,
'B_field_cutoff': 5.,
'B_ran_b2': 0.1,
'B_ran_h_d': 4.,
'B_ran_h_d2': 1.,
'B_ran_h_r': 15.,
'B_ran_h_r2': 15.,
'bb_amps_0': 2.,
'bb_amps_1': 0.133,
'bb_amps_2': -3.78,
'bb_amps_3': 0.32,
'bb_amps_4': -0.023,
'bb_bar_a': 0.,
'bb_bar_boa': 1.,
'bb_bar_phi0_deg': 45.,
'bb_cr0_coh': 0.5,
'bb_cr0_iso': 0.3,
'bb_cr0_ord': 0.5,
'bb_d0_iso': 0.3,
'bb_d0': 0.3,
'bb_delta_phi_iso_deg': 0.,
'bb_delta_phi_ord_deg': 0.,
'bb_disk_b0': 0.167,
'bb_disk_h_d': 0.1,
'bb_halo_b0': 1.38,
'bb_halo_h_d': 3.,
'bb_phi0_deg': 70.,
'bb_pitch_biso': -11.,
'bb_pitch': -11.5,
'bb_r_compconst': 12.,
'bb_r_innercut': 0.5,
'bb_r_peak': -1.,
'bb_r_scale': 20.,
'bb_rmax_arms': 30.,
'bb_spiral_cpow': 3.,
'bb_spiral_h_d': 0.1,
}
return defaults
@property
def _initial_variable_to_parameter_mappings(self):
return self._generate_variable_to_parameter_mapping_defaults(n=3)
def _generate_variable_to_parameter_mapping_defaults(self, n):
defaults = {'B_f_ord': self._interval(0.5, 1, n),
'B_field_RMS_uG': self._interval(1., 0.3, n),
'B_field_alpha': self._interval(-2.37, 1., n),
'B_field_cutoff': self._interval(5., 1., n),
'B_ran_b2': self._interval(0.1, 0.3, n),
'B_ran_h_d': self._interval(4., 1., n),
'B_ran_h_d2': self._interval(1., 1., n),
'B_ran_h_r': self._interval(15., 4., n),
'B_ran_h_r2': self._interval(15., 4., n),
'bb_amps_0': self._interval(2., 1., n),
'bb_amps_1': self._interval(0.133, 1., n),
'bb_amps_2': self._interval(-3.78, 1., n),
'bb_amps_3': self._interval(0.32, 1., n),
'bb_amps_4': self._interval(-0.023, 1., n),
'bb_bar_a': self._interval(0., 1., n),
'bb_bar_boa': self._interval(1., 1., n),
'bb_bar_phi0_deg': self._interval(45., 15., n),
'bb_cr0_coh': self._interval(0.5, 1., n),
'bb_cr0_iso': self._interval(0.3, 1., n),
'bb_cr0_ord': self._interval(0.5, 1., n),
'bb_d0_iso': self._interval(0.3, 0.1, n),
'bb_d0': self._interval(0.3, 0.1, n),
'bb_delta_phi_iso_deg': self._interval(0., 1., n),
'bb_delta_phi_ord_deg': self._interval(0., 1., n),
'bb_disk_b0': self._interval(0.167, 1., n),
'bb_disk_h_d': self._interval(0.1, 1., n),
'bb_halo_b0': self._interval(1.38, 1., n),
'bb_halo_h_d': self._interval(3., 0.9, n),
'bb_phi0_deg': self._interval(70., 20., n),
'bb_pitch_biso': self._interval(-11, 5., n),
'bb_pitch': self._interval(-11.5, 5., n),
'bb_r_compconst': self._interval(12., 3., n),
'bb_r_innercut': self._interval(0.5, 1., n),
'bb_r_peak': self._interval(-1., 1., n),
'bb_r_scale': self._interval(20., 4., n),
'bb_rmax_arms': self._interval(30., 5., n),
'bb_spiral_cpow': self._interval(3., 1., n),
'bb_spiral_h_d': self._interval(0.1, 1., n),
}
return defaults
imagine/magnetic_fields/jf12_magnetic_field/__init__.py deleted 100644 → 0
View file @ 7d179bd2
# -*- coding: utf-8 -*-
from jf12_magnetic_field import JF12MagneticField
from jf12_magnetic_field_factory import JF12MagneticFieldFactory
\ No newline at end of file
imagine/magnetic_fields/jf12_magnetic_field/jf12_magnetic_field.py deleted 100644 → 0
View file @ 7d179bd2
# -*- coding: utf-8 -*-
from imagine.magnetic_fields.magnetic_field import MagneticField
class JF12MagneticField(MagneticField):
@property
def parameter_list(self):
parameter_list = ['b51_ran_b1', 'b51_ran_b2', 'b51_ran_b3',
'b51_ran_b4', 'b51_ran_b5', 'b51_ran_b6',
'b51_ran_b7', 'b51_ran_b8', 'b51_coh_b1',
'b51_coh_b2', 'b51_coh_b3', 'b51_coh_b4',
'b51_coh_b5', 'b51_coh_b6', 'b51_coh_b7',
'b51_z0_spiral', 'b51_z0_smooth', 'b51_r0_smooth',
'b51_b0_smooth', 'b51_b0_x', 'b51_Xtheta',
'b51_r0_x', 'b51_h_disk', 'b51_Bn', 'b51_Bs',
'b51_z0_halo', 'b51_b_ring', 'b51_b0_interior',
'b51_shift', 'B_analytic_beta']
return parameter_list
def _create_field(self):
raise NotImplementedError
imagine/magnetic_fields/jf12_magnetic_field/jf12_magnetic_field_factory.py deleted 100644 → 0
View file @ 7d179bd2
# -*- coding: utf-8 -*-
from imagine.magnetic_fields.magnetic_field.magnetic_field_factory \
import MagneticFieldFactory
from jf12_magnetic_field import JF12MagneticField
class JF12MagneticFieldFactory(MagneticFieldFactory):
@property
def magnetic_field_class(self):
return JF12MagneticField
@property
def _initial_parameter_defaults(self):
defaults = {'b51_ran_b1': 10.8,
'b51_ran_b2': 6.96,
'b51_ran_b3': 9.59,
'b51_ran_b4': 6.96,
'b51_ran_b5': 1.96,
'b51_ran_b6': 16.34,
'b51_ran_b7': 37.29,
'b51_ran_b8': 10.35,
'b51_coh_b1': 0.1,
'b51_coh_b2': 3.,
'b51_coh_b3': -0.9,
'b51_coh_b4': -0.8,
'b51_coh_b5': -2.,
'b51_coh_b6': -4.2,
'b51_coh_b7': 0.,
'b51_z0_spiral': 0.61,
'b51_z0_smooth': 2.84,
'b51_r0_smooth': 10.97,
'b51_b0_smooth': 4.68,
'b51_b0_x': 4.6,
'b51_Xtheta': 49.,
'b51_r0_x': 2.9,
'b51_h_disk': 0.4,
'b51_Bn': 1.4,
'b51_Bs': -1.1,
'b51_z0_halo': 5.3,
'b51_b_ring': 0.1,
'b51_b0_interior': 7.63,
'b51_shift': 0.,
'B_analytic_beta': 1.36}
return defaults
@property
def _initial_variable_to_parameter_mappings(self):
return self._generate_variable_to_parameter_mapping_defaults(n=3)
def _generate_variable_to_parameter_mapping_defaults(self, n):
defaults = {
'b51_ran_b1': self._interval(10.8, 2.33, n), # b_1, 1210.7820
'b51_ran_b2': self._interval(6.96, 1.58, n), # b_2, 1210.7820
'b51_ran_b3': self._interval(9.59, 1.10, n), # b_3, 1210.7820
'b51_ran_b4': self._interval(6.96, 0.87, n), # b_4 1210.7820
'b51_ran_b5': self._interval(1.96, 1.32, n), # b_5, 1210.7820
'b51_ran_b6': self._interval(16.34, 2.53, n), # b_6, 1210.7820
'b51_ran_b7': self._interval(37.29, 2.39, n), # b_7, 1210.7820
'b51_ran_b8': self._interval(10.35, 4.43, n), # b_8, 1210.7820
'b51_b0_interior': self._interval(7.63, 1.39, n), # b_int, 1210.7820
'b51_z0_spiral': self._positive_interval(0.61, 0.04, n), # z_0^disk, 1210.7820
'b51_b0_smooth': self._interval(4.68, 1.39, n), # B_0, 1210.7820
'b51_r0_smooth': self._positive_interval(10.97, 3.80, n), # r_0, 1210.7820
'b51_z0_smooth': self._positive_interval(2.84, 1.30, n), # z_0, 1210.7820
# beta is missing (1210.7820)
'b51_coh_b1': self._interval(0.1, 1.8, n), # b_1, 1204.3662
'b51_coh_b2': self._interval(3.0, 0.6, n), # b_2, 1204.3662
'b51_coh_b3': self._interval(-0.9, 0.8, n), # b_3, 1204.3662
'b51_coh_b4': self._interval(-0.8, 0.3, n), # b_4, 1204.3662
'b51_coh_b5': self._interval(-2.0, 0.1, n), # b_5, 1204.3662
'b51_coh_b6': self._interval(-4.2, 0.5, n), # b_6, 1204.3662
'b51_coh_b7': self._interval(0.0, 1.8, n), # b_7, 1204.3662
'b51_b_ring': self._interval(0.1, 0.1, n), # b_ring, 1204.3662
'b51_h_disk': self._positive_interval(0.4, 0.03, n), # b_disk, 1204.3662
# w_disk is missing, 1204.3662
# B_n_disk is missing, 1204.3662
'b51_Bs': self._interval(-1.1, 0.1, n), # B_s, 1204.3662
# r_n is missing, 1204.3662
# r_s is missing, 1204.3662
# w_h is missing, 1204.3662
'b51_z0_halo': self._positive_interval(5.3, 1.6, n), # z_0, 1204.3662
'b51_b0_x': self._interval(4.6, 0.3, n), # B_X, 1204.3662
'b51_Xtheta': self._positive_interval(49., 1., n), # Theta_X^0, 1204.3662
# r_x^c is missing, 1204.3662
'b51_r0_x': self._positive_interval(2.9, 0.1, n), # r_x, 1204.3662
# gamma is missing (1204.3662)
'b51_shift': [-2, 0, 2],
'B_analytic_beta': self._interval(1.36, 0.4, n),
}
return defaults
imagine/magnetic_fields/magnetic_field/magnetic_field.py
View file @ 6c8d5721
......@@ -31,7 +31,7 @@ class MagneticField(Field):
assert(isinstance(self.domain[2], FieldArray))
self._parameters = {}
for p in self.parameter_list:
for p in self.descriptor_lookup:
self._parameters[p] = np.float(parameters[p])
casted_random_seed = distributed_data_object(
......@@ -48,8 +48,8 @@ class MagneticField(Field):
self.random_seed = casted_random_seed
@property
def parameter_list(self):
return []
def descriptor_lookup(self):
return {}
@property
def parameters(self):
......
imagine/magnetic_fields/wmap3yr_magnetic_field/wmap3yr_magnetic_field.py
View file @ 6c8d5721
......@@ -4,14 +4,18 @@ from imagine.magnetic_fields.magnetic_field import MagneticField
class WMAP3yrMagneticField(MagneticField):
@property
def parameter_list(self):
parameter_list = ['B_field_b0',
'B_field_psi0_deg',
'B_field_psi1_deg',
'B_field_xsi0_deg',
'B_field_RMS_uG']
return parameter_list
def descriptor_lookup(self):
lookup = \
{'b0': ('./Galaxy/MagneticField/Regular/WMAP/b0', 'value'),
'psi0': ('./Galaxy/MagneticField/Regular/WMAP/psi0', 'value'),
'psi1': ('./Galaxy/MagneticField/Regular/WMAP/psi1', 'value'),
'chi0': ('./Galaxy/MagneticField/Regular/WMAP/chi0', 'value'),
'random_rms': ('./Galaxy/MagneticField/Random/Iso/rms', 'value'),
'random_rho': ('./Galaxy/MagneticField/Random/Anisoglob/rho',
'value')}
return lookup
def _create_field(self):
raise NotImplementedError
imagine/magnetic_fields/wmap3yr_magnetic_field/wmap3yr_magnetic_field_factory.py
View file @ 6c8d5721
......@@ -13,13 +13,12 @@ class WMAP3yrMagneticFieldFactory(MagneticFieldFactory):
@property
def _initial_parameter_defaults(self):
defaults = {
'B_field_b0': 6,
'B_field_psi0_deg': 27,
'B_field_psi1_deg': 0.9,
'B_field_xsi0_deg': 25,
'B_field_RMS_uG': 1.0,
}
defaults = {'b0': 6,
'psi0': 27,
'psi1': 0.9,
'chi0': 25,
'random_rms': 1.0,
'random_rho': 0.5}
return defaults
@property
......@@ -28,10 +27,11 @@ class WMAP3yrMagneticFieldFactory(MagneticFieldFactory):
def _generate_variable_to_parameter_mapping_defaults(self, n):
defaults = {
'B_field_b0': self._positive_interval(6.0, 2.0, n), # b0 astro-ph/0603450
'B_field_psi0_deg': self._positive_interval(27.0, 5.0, n), # psi0 astro-ph/0603450
'B_field_psi1_deg': self._positive_interval(0.9, 5.0, n), # psi1 astro-ph/0603450
'B_field_xsi0_deg': self._positive_interval(25, 5.0, n), # xsi0 astro-ph/0603450
'B_field_RMS_uG': self._positive_interval(1.0, 2.0, n)
'b0': self._positive_interval(6.0, 2.0, n), # b0 astro-ph/0603450
'psi0': self._positive_interval(27.0, 5.0, n), # psi0 astro-ph/0603450
'psi1': self._positive_interval(0.9, 5.0, n), # psi1 astro-ph/0603450
'chi0': self._positive_interval(25, 5.0, n), # xsi0 astro-ph/0603450
'random_rms': self._positive_interval(1.0, 2.0, n),
'random_rho': self._positive_interval(0.5, 1/6., n)
}
return defaults
imagine/observers/hammurapy/hammurapy.py
View file @ 6c8d5721
......@@ -4,6 +4,7 @@ import abc
import os
import tempfile
import subprocess
import xml.etree.ElementTree as et
import healpy
import numpy as np
......@@ -15,48 +16,16 @@ from imagine.magnetic_fields.magnetic_field import MagneticField
class Hammurapy(Observer):
def __init__(self, hammurabi_executable, conf_directory='./confs',
working_directory_base='.', nside=128,
analytic_ensemble_mean=False):
def __init__(self, hammurabi_executable, input_directory='./input',
working_directory_base='.', nside=64):
self.hammurabi_executable = os.path.abspath(hammurabi_executable)
self.conf_directory = os.path.abspath(conf_directory)
self.input_directory = os.path.abspath(input_directory)
self.working_directory_base = os.path.abspath(working_directory_base)
self.analytic_ensemble_mean = bool(analytic_ensemble_mean)
self.nside = int(nside)
self.last_call_log = ""
sync_template_fname = os.path.join(self.conf_directory,
'IQU_sync.fits')
dust_template_fname = os.path.join(self.conf_directory,
'IQU_dust.fits')
self.basic_parameters = {'B_ran_mem_lim': '6',
'obs_shell_index_numb': '1',
'total_shell_numb': '1',
'vec_size_R': '400',
'max_radius': '30',
'max_z': '10',
'B_field_transform_lon': '-999',
'B_field_transform_lat': '-999',
'TE_grid_filename': 'negrid_n400.bin',
'TE_nx': '400',
'TE_ny': '400',
'TE_nz': '80',
'TE_interp': 'T',
'do_sync_emission': 'F',
'do_rm': 'F',
'do_dm': 'F',
'do_dust': 'F',
'do_tau': 'F',
'do_ff': 'F',
'obs_freq_GHz': '23',
'sync_template_fname': sync_template_fname,
'dust_template_fname': dust_template_fname
}
@abc.abstractproperty
def magnetic_field_class(self):
return MagneticField
......@@ -89,7 +58,7 @@ class Hammurapy(Observer):
try:
loaded_map = healpy.read_map(map_path, verbose=False,
field=i)
loaded_map = healpy.ud_grade(loaded_map, nside_out=nside)
# loaded_map = healpy.ud_grade(loaded_map, nside_out=nside)
result_list += [loaded_map]
i += 1
except IndexError:
......@@ -97,54 +66,73 @@ class Hammurapy(Observer):
return result_list
def _call_hammurabi(self, path):
parameters_file_path = os.path.join(path, 'parameters.txt')
temp_process = subprocess.Popen([self.hammurabi_executable,
parameters_file_path],
stdout=subprocess.PIPE,
cwd=self.conf_directory)
temp_process = subprocess.Popen(
[self.hammurabi_executable, 'parameters.xml'],
stdout=subprocess.PIPE,
cwd=path)
self.last_call_log = temp_process.communicate()[0]
def _initialize_observable_dict(self, observable_dict, magnetic_field):
pass
def _build_parameter_dict(self, parameter_dict, magnetic_field,
working_directory, local_ensemble_index):
local_ensemble_index):
parameter_dict.update(
{('./Interface/fe_grid', 'read'): '1',
('./Interface/fe_grid', 'filename'):
os.path.join(self.input_directory, 'fe_grid.bin'),
})
# access the magnetic-field's random-seed d2o directly, since we
# know that the distribution strategy is the same for the
# randam samples and the magnetic field itself
random_seed = magnetic_field.random_seed.data[local_ensemble_index]
parameter_dict.update(
{('./Galaxy/MagneticField/Random', 'seed'): random_seed})
for key, value in magnetic_field.parameters.iteritems():
large_key = magnetic_field.descriptor_lookup[key]
parameter_dict[large_key] = value
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
if local_ensemble_index is not None:
# access the magnetic-field's random-seed d2o directly, since we
# know that the distribution strategy is the same for the
# randam samples and the magnetic field itself
random_seed = magnetic_field.random_seed.data[local_ensemble_index]
parameter_dict.update({'B_field_seed': random_seed})
parameter_dict.update({'B_field_lx': lx,
'B_field_ly': ly,
'B_field_lz': lz,
'B_field_nx': nx,
'B_field_ny': ny,
'B_field_nz': nz,
})
parameter_dict.update({'obs_NSIDE': self.nside})
parameter_dict.update({('./Grid/Box/lx', 'value'): lx,
('./Grid/Box/ly', 'value'): ly,
('./Grid/Box/lz', 'value'): lz,