Commit 7e02469c authored by Jiaxin Wang's avatar Jiaxin Wang
Browse files

update wrt hamx

parent 9476279e
......@@ -18,8 +18,8 @@
from imagine.magnetic_fields.magnetic_field import MagneticField
class ConstantMagneticField(MagneticField):
@property
def parameter_list(self):
parameter_list = ['b_x', 'b_y', 'b_z']
......
......@@ -17,14 +17,11 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from imagine.magnetic_fields.magnetic_field.magnetic_field_factory \
import MagneticFieldFactory
from imagine.magnetic_fields.magnetic_field.magnetic_field_factory import MagneticFieldFactory
from constant_magnetic_field import ConstantMagneticField
class ConstantMagneticFieldFactory(MagneticFieldFactory):
@property
def magnetic_field_class(self):
return ConstantMagneticField
......
......@@ -17,9 +17,7 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from mpi4py import MPI
import simplejson as json
import numpy as np
from nifty import Field, FieldArray, RGSpace
......@@ -29,17 +27,16 @@ comm = MPI.COMM_WORLD
size = comm.size
rank = comm.rank
class MagneticField(Field):
def __init__(self, parameters={}, domain=None, val=None, dtype=None,
distribution_strategy=None, copy=False, random_seed=None):
super(MagneticField, self).__init__(
domain=domain,
val=val,
dtype=dtype,
distribution_strategy=distribution_strategy,
copy=copy)
super(MagneticField, self).__init__(domain=domain,
val=val,
dtype=dtype,
distribution_strategy=distribution_strategy,
copy=copy)
assert(len(self.domain) == 3)
assert(isinstance(self.domain[0], FieldArray))
......@@ -47,19 +44,17 @@ class MagneticField(Field):
assert(isinstance(self.domain[2], FieldArray))
self._parameters = {}
for p in self.descriptor_lookup:
for p in self.descriptor_lookup:#{
self._parameters[p] = np.float(parameters[p])
casted_random_seed = distributed_data_object(
global_shape=(self.shape[0],),
dtype=np.int,
distribution_strategy=distribution_strategy)
if random_seed is None:
#}
casted_random_seed = distributed_data_object(global_shape=(self.shape[0],),
dtype=np.int,
distribution_strategy=distribution_strategy)
if random_seed is None:#{
random_seed = np.random.randint(np.uint32(-1)/3,
size=self.shape[0])
random_seed = comm.bcast(random_seed, root=0)
#}
casted_random_seed[:] = random_seed
self.random_seed = casted_random_seed
......@@ -72,9 +67,9 @@ class MagneticField(Field):
return self._parameters
def set_val(self, new_val=None, copy=False):
if new_val is not None:
raise RuntimeError("Setting the field values explicitly is not "
"supported by MagneticField.")
if new_val is not None:#{
raise RuntimeError("Setting the field values explicitly is not supported by MagneticField.")
#}
self._val = self._create_field()
def _to_hdf5(self, hdf5_group):
......
......@@ -17,16 +17,11 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from keepers import Loggable
from imagine.carrier_mapper import unity_mapper
from nifty import FieldArray, RGSpace
from magnetic_field import MagneticField
class MagneticFieldFactory(Loggable, object):
def __init__(self, box_dimensions, resolution):
......@@ -34,9 +29,7 @@ class MagneticFieldFactory(Loggable, object):
self.box_dimensions = box_dimensions
self.resolution = resolution
self._parameter_defaults = self._initial_parameter_defaults
self._variable_to_parameter_mappings = \
self._initial_variable_to_parameter_mappings
self._variable_to_parameter_mappings = self._initial_variable_to_parameter_mappings
distances = np.array(self.box_dimensions) / np.array(self.resolution)
self._grid_space = RGSpace(shape=self.resolution,
distances=distances)
......@@ -44,9 +37,9 @@ class MagneticFieldFactory(Loggable, object):
self._ensemble_cache = {}
def _get_ensemble(self, ensemble_size):
if ensemble_size not in self._ensemble_cache:
self._ensemble_cache[ensemble_size] = \
FieldArray(shape=(ensemble_size,))
if ensemble_size not in self._ensemble_cache:#{
self._ensemble_cache[ensemble_size] = FieldArray(shape=(ensemble_size,))
#}
return self._ensemble_cache[ensemble_size]
@property
......@@ -56,8 +49,9 @@ class MagneticFieldFactory(Loggable, object):
@box_dimensions.setter
def box_dimensions(self, box_dimensions):
dim = tuple(np.array(box_dimensions, dtype=np.float))
if len(dim) != 3:
if len(dim) != 3:#{
raise ValueError("Input of box_dimensions must have length three.")
#}
self._box_dimensions = dim
@property
......@@ -67,8 +61,9 @@ class MagneticFieldFactory(Loggable, object):
@resolution.setter
def resolution(self, resolution):
resolution = tuple(np.array(resolution, dtype=np.int))
if len(resolution) != 3:
if len(resolution) != 3:#{
raise ValueError("Input for resolution must have length three.")
#}
self._resolution = resolution
@property
......@@ -104,10 +99,11 @@ class MagneticFieldFactory(Loggable, object):
@property
def variable_defaults(self):
variable_defaults = {}
for parameter in self.parameter_defaults:
for parameter in self.parameter_defaults:#{
low, high = self.variable_to_parameter_mappings[parameter]
default = self.parameter_defaults[parameter]
variable_defaults[parameter] = (default - low)/(high - low)
#}
return variable_defaults
@property
......@@ -122,17 +118,17 @@ class MagneticFieldFactory(Loggable, object):
value: [min, max]
"""
for k, v in new_mapping.items():
if k in self._variable_to_parameter_mappings:
if k in self._variable_to_parameter_mappings:#{
key = str(k)
value = [np.float(v[0]), np.float(v[1])]
self._variable_to_parameter_mappings.update({key: value})
self.logger.debug("Updated variable_to_parameter_mapping %s "
"to %s" % (key, str(value)))
self.logger.debug("Updated variable_to_parameter_mapping %s to %s" % (key, str(value)))
#}
def _map_variables_to_parameters(self, variables):
parameter_dict = {}
for variable_name in variables:
if variable_name in self.variable_to_parameter_mappings:
if variable_name in self.variable_to_parameter_mappings:#{
mapping = self.variable_to_parameter_mappings[variable_name]
mapped_variable = unity_mapper(variables[variable_name],
a=mapping[0],
......@@ -141,8 +137,10 @@ class MagneticFieldFactory(Loggable, object):
# a=mapping[0],
# m=mapping[1],
# b=mapping[2])
else:
#}
else:#{
mapped_variable = np.float(variables[variable_name])
#}
parameter_dict[variable_name] = mapped_variable
return parameter_dict
......@@ -150,16 +148,12 @@ class MagneticFieldFactory(Loggable, object):
mapped_variables = self._map_variables_to_parameters(variables)
work_parameters = self.parameter_defaults.copy()
work_parameters.update(mapped_variables)
domain = (self._get_ensemble(ensemble_size),
self._grid_space,
self._vector)
result_magnetic_field = self.magnetic_field_class(
domain=domain,
parameters=work_parameters,
distribution_strategy='equal',
random_seed=random_seed)
self.logger.debug("Generated magnetic field with work-parameters %s" %
work_parameters)
result_magnetic_field = self.magnetic_field_class(domain=domain,
parameters=work_parameters,
distribution_strategy='equal',
random_seed=random_seed)
self.logger.debug("Generated magnetic field with work-parameters %s" % work_parameters)
return result_magnetic_field
......@@ -18,20 +18,21 @@
from imagine.magnetic_fields.magnetic_field import MagneticField
class WMAP3yrMagneticField(MagneticField):
@property
def descriptor_lookup(self):
lookup = \
{'b0': ['./MagneticField/Regular/WMAP/b0', 'value'],
'psi0': ['./MagneticField/Regular/WMAP/psi0', 'value'],
'psi1': ['./MagneticField/Regular/WMAP/psi1', 'value'],
'chi0': ['./MagneticField/Regular/WMAP/chi0', 'value'],
'random_rms': ['./MagneticField/Random/Global/rms', 'value'],
'random_rho': ['./MagneticField/Random/Global/rho', 'value'],
'random_a0': ['./MagneticField/Random/Global/a0', 'value']
}
lookup = {'b0': ['./MagneticField/Regular/WMAP/b0', 'value'],
'psi0': ['./MagneticField/Regular/WMAP/psi0', 'value'],
'psi1': ['./MagneticField/Regular/WMAP/psi1', 'value'],
'chi0': ['./MagneticField/Regular/WMAP/chi0', 'value'],
'random_rms': ['./MagneticField/Random/Global/ES/rms', 'value'],
'random_rho': ['./MagneticField/Random/Global/ES/rho', 'value'],
'random_a0': ['./MagneticField/Random/Global/ES/a0', 'value'],
'random_k0': ['./MagneticField/Random/Global/ES/k0', 'value'],
'random_r0': ['./MagneticField/Random/Global/ES/r0', 'value'],
'random_z0': ['./MagneticField/Random/Global/ES/z0', 'value']
}
return lookup
def _create_field(self):
......
......@@ -16,13 +16,11 @@
# IMAGINE is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from imagine.magnetic_fields.magnetic_field.magnetic_field_factory \
import MagneticFieldFactory
from imagine.magnetic_fields.magnetic_field.magnetic_field_factory import MagneticFieldFactory
from wmap3yr_magnetic_field import WMAP3yrMagneticField
class WMAP3yrMagneticFieldFactory(MagneticFieldFactory):
@property
def magnetic_field_class(self):
return WMAP3yrMagneticField
......@@ -33,9 +31,13 @@ class WMAP3yrMagneticFieldFactory(MagneticFieldFactory):
'psi0': 27.0,
'psi1': 0.9,
'chi0': 25.0,
# using ES random field generator in hammurabiX
'random_rms': 2.0,
'random_rho': 0.5,
'random_a0': 1.7}
'random_a0': 1.7,
'random_k0': 0.5,
'random_r0': 8.0,
'random_z0': 1.0}
return defaults
@property
......@@ -43,13 +45,15 @@ class WMAP3yrMagneticFieldFactory(MagneticFieldFactory):
return self._generate_variable_to_parameter_mapping_defaults(n=3)
def _generate_variable_to_parameter_mapping_defaults(self, n):
defaults = {
'b0': self._positive_interval(6.0, 1.9, n), # b0 astro-ph/0603450
'psi0': self._positive_interval(27.0, 7.0, n), # psi0 astro-ph/0603450
'psi1': self._positive_interval(0.9, 5.0, n), # psi1 astro-ph/0603450
'chi0': self._positive_interval(25, 8.0, n), # xsi0 astro-ph/0603450
'random_rms': self._positive_interval(2.0, 0.6, n),
'random_rho': self._positive_interval(0.5, 0.166, n),
'random_a0': self._positive_interval(1.7, 0.5, n)
}
defaults = {'b0': self._positive_interval(6.0, 1.9, n), # b0 astro-ph/0603450
'psi0': self._positive_interval(27.0, 7.0, n), # psi0 astro-ph/0603450
'psi1': self._positive_interval(0.9, 5.0, n), # psi1 astro-ph/0603450
'chi0': self._positive_interval(25, 8.0, n), # xsi0 astro-ph/0603450
'random_rms': self._positive_interval(2.0, 0.6, n),
'random_rho': self._positive_interval(0.5, 0.166, n),
'random_a0': self._positive_interval(1.7, 0.5, n),
'random_k0': self._positive_interval(0.5, 0.15, n),
'random_r0': self._positive_interval(8.0, 2.0, n),
'random_z0': self._positive_interval(1.0, 0.3, n)
}
return defaults
......@@ -18,18 +18,15 @@
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)
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))
......@@ -42,8 +39,8 @@ class Observable(Field):
return self._ensemble_mean
def _to_hdf5(self, hdf5_group):
if hasattr(self, '_ensemble_mean'):
return_dict = {'ensemble_mean': self._ensemble_mean}
if hasattr(self, "_ensemble_mean"):
return_dict = {"ensemble_mean": self._ensemble_mean}
else:
return_dict = {}
return_dict.update(
......@@ -55,7 +52,7 @@ class Observable(Field):
new_field = super(Observable, cls)._from_hdf5(hdf5_group=hdf5_group,
repository=repository)
try:
observable_mean = repository.get('ensemble_mean', hdf5_group)
observable_mean = repository.get("ensemble_mean", hdf5_group)
new_field._observable_mean = observable_mean
except(KeyError):
pass
......
......@@ -20,11 +20,9 @@ import os
import tempfile
import subprocess
import xml.etree.ElementTree as et
import numpy as np
from d2o import distributed_data_object
from nifty import HPSpace
from imagine.observers.observer import Observer
......@@ -32,35 +30,30 @@ from .observable_mixins import ObservableMixin
from .model_mixins import MagneticFieldModel
from imagine.observables import Observable
class Hammurapy(Observer):
def __init__(self, hammurabi_executable, magnetic_field_model, observables,
input_directory='./input', working_directory_base='.',
nside=64):
input_directory='./input', working_directory_base='.',nside=64):
self.hammurabi_executable = os.path.abspath(hammurabi_executable)
if not isinstance(magnetic_field_model, MagneticFieldModel):
raise TypeError("magnetic_field_model must be an instance of the "
"MagneticField class.")
if not isinstance(magnetic_field_model, MagneticFieldModel):#{
raise TypeError("magnetic_field_model must be an instance of the MagneticField class.")
#}
self.magnetic_field_model = magnetic_field_model
if not isinstance(observables, list):
if not isinstance(observables, list):#{
if isinstance(observables, tuple):
observables = list(observables)
else:
observables = [observables]
for obs in observables:
#}
for obs in observables:#{
if not isinstance(obs, ObservableMixin):
raise TypeError("observables must be an instance of the "
"ObservableMixin class.")
raise TypeError("observables must be an instance of the ObservableMixin class.")
#}
self.observables = observables
self.input_directory = os.path.abspath(input_directory)
self.working_directory_base = os.path.abspath(working_directory_base)
self.nside = int(nside)
self._hpSpace = HPSpace(nside=self.nside)
self.last_call_log = ""
@property
......@@ -75,7 +68,7 @@ class Hammurapy(Observer):
# Try multiple times in order to overcome 'Directory not empty' errors
# Hopefully open file handles get closed in the meantime
n = 0
while n < 10:
while n < 10:#{
temp_process = subprocess.Popen(['rm', '-rf', path],
stderr=subprocess.PIPE)
# wait until process is finished
......@@ -84,21 +77,22 @@ class Hammurapy(Observer):
if errlog == '':
self.logger.debug("Successfully removed temporary folder.")
break
#}
else:
self.logger.warning('Could not delete %s' % path)
def _call_hammurabi(self, path):
temp_process = subprocess.Popen(
[self.hammurabi_executable, 'parameters.xml'],
stdout=subprocess.PIPE,
cwd=path)
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, magnetic_field):
observable_dict = {}
for observable in self.observables:
observable_dict = self._initialize_observable_dict_helper(
observable_dict, magnetic_field, observable)
observable_dict = self._initialize_observable_dict_helper(observable_dict,
magnetic_field,
observable)
return observable_dict
def _initialize_observable_dict_helper(self, observable_dict,
......@@ -108,39 +102,31 @@ class Hammurapy(Observer):
# It is important to initialize the Observables with an explicit
# value. Otherwise the d2o will not instantaneuosly be created
# (c.f. lazy object creation).
observable_dict[component] = Observable(
val=0,
domain=(ensemble_space, self._hpSpace),
distribution_strategy='equal')
observable_dict[component] = Observable(val=0,
domain=(ensemble_space, self._hpSpace),
distribution_strategy='equal')
return observable_dict
def _write_parameter_xml(self, magnetic_field, local_ensemble_index,
working_directory):
# load the default xml
try:
try:#{
default_parameters_xml = os.path.join(self.input_directory,
'default_parameters.xml')
tree = et.parse(default_parameters_xml)
except IOError:
#}
except IOError:#{
import imagine
module_path = os.path.split(
imagine.observers.hammurapy.__file__)[0]
default_parameters_xml = os.path.join(
module_path, 'input/default_parameters.xml')
module_path = os.path.split(imagine.observers.hammurapy.__file__)[0]
default_parameters_xml = os.path.join(module_path, 'input/default_parameters.xml')
tree = et.parse(default_parameters_xml)
#}
root = tree.getroot()
# modify the default_parameters.xml
custom_parameters = [
['./Grid/Integration/shell/auto/nside_min', 'value',
self.nside],
['./Grid/Integration/nside_sim', 'value', self.nside],
# ['./Interface/fe_grid', 'read', '1'],
# ['./Interface/fe_grid', 'filename',
# os.path.join(self.input_directory, 'fe_grid.bin')]
custom_parameters = [['./Grid/Shell/layer/auto/nside_min', 'value', self.nside],
['./Grid/Shell/nside_sim', 'value', self.nside],
]
# access the magnetic-field's random-eed d2o directly, since we
# know that the distribution strategy is the same for the
# randam samples and the magnetic field itself
......@@ -156,27 +142,24 @@ class Hammurapy(Observer):
# set up grid parameters
grid_space = magnetic_field.domain[1]
lx, ly, lz = \
np.array(grid_space.shape)*np.array(grid_space.distances)/2.
lx, ly, lz = np.array(grid_space.shape)*np.array(grid_space.distances)/2.
nx, ny, nz = grid_space.shape
custom_parameters += [['./Grid/Box/x_min', 'value', -lx],
['./Grid/Box/x_max', 'value', lx],
['./Grid/Box/y_min', 'value', -ly],
['./Grid/Box/y_max', 'value', ly],
['./Grid/Box/z_min', 'value', -lz],
['./Grid/Box/z_max', 'value', lz],
['./Grid/Box/nx', 'value', nx],
['./Grid/Box/ny', 'value', ny],
['./Grid/Box/nz', 'value', nz]]
for parameter in custom_parameters:
custom_parameters += [['./Grid/Box_GMF/x_min', 'value', -lx],
['./Grid/Box_GMF/x_max', 'value', lx],
['./Grid/Box_GMF/y_min', 'value', -ly],
['./Grid/Box_GMF/y_max', 'value', ly],
['./Grid/Box_GMF/z_min', 'value', -lz],
['./Grid/Box_GMF/z_max', 'value', lz],
['./Grid/Box_GMF/nx', 'value', nx],
['./Grid/Box_GMF/ny', 'value', ny],
['./Grid/Box_GMF/nz', 'value', nz]]
for parameter in custom_parameters:#{
root.find(parameter[0]).set(parameter[1], str(parameter[2]))
#}
self.magnetic_field_model.update_parameter_xml(root)
for observable in self.observables:
for observable in self.observables:#{
observable.update_parameter_xml(root)
#}
parameters_file_path = os.path.join(working_directory,
'parameters.xml')
tree.write(parameters_file_path)
......@@ -184,54 +167,41 @@ class Hammurapy(Observer):
def _fill_observable_dict(self, observable_dict, working_directory,
local_ensemble_index):
for observable in self.observables:
observable.fill_observable_dict(
observable_dict=observable_dict,
working_directory=working_directory,
local_ensemble_index=local_ensemble_index,
nside=self.nside)
observable.fill_observable_dict(observable_dict=observable_dict,
working_directory=working_directory,
local_ensemble_index=local_ensemble_index,
nside=self.nside)
return observable_dict
def __call__(self, magnetic_field):
if not isinstance(magnetic_field, self.magnetic_field_class):
raise ValueError("Given magnetic field is not a subclass of" +
" %s" % str(self.magnetic_field_class))
observable_dict = self._initialize_observable_dict(
magnetic_field=magnetic_field)
raise ValueError("Given magnetic field is not a subclass of %s" % str(self.magnetic_field_class))
observable_dict = self._initialize_observable_dict(magnetic_field=magnetic_field)
# iterate over ensemble and put result into result_observable
# get the local shape by creating a dummy d2o
ensemble_number = magnetic_field.shape[0]
dummy = distributed_data_object(global_shape=(ensemble_number,),
distribution_strategy='equal',
dtype=np.float)
local_length = dummy.distributor.local_length
for local_ensemble_index in xrange(local_length):
self.logger.debug("Processing local_ensemble_index %i." %
local_ensemble_index)
self.logger.debug("Processing local_ensemble_index %i." % local_ensemble_index)
# create a temporary folder
working_directory = self._make_temp_folder()
self._write_parameter_xml(