Commit 3a42a026 authored by Jait Dixit's avatar Jait Dixit
Browse files

Merge branch 'feature/field_multiple_space' into smooth_operator

parents 871c8d12 31a4e7e3
...@@ -80,8 +80,8 @@ variable_default_field_dtype = keepers.Variable( ...@@ -80,8 +80,8 @@ variable_default_field_dtype = keepers.Variable(
_dtype_validator, _dtype_validator,
genus='str') genus='str')
variable_default_datamodel = keepers.Variable( variable_default_distribution_strategy = keepers.Variable(
'default_datamodel', 'default_distribution_strategy',
['fftw', 'equal'], ['fftw', 'equal'],
lambda z: (('pyfftw' in dependency_injector) lambda z: (('pyfftw' in dependency_injector)
if z == 'fftw' else True), if z == 'fftw' else True),
...@@ -95,7 +95,7 @@ nifty_configuration = keepers.get_Configuration( ...@@ -95,7 +95,7 @@ nifty_configuration = keepers.get_Configuration(
variable_use_libsharp, variable_use_libsharp,
variable_verbosity, variable_verbosity,
variable_default_field_dtype, variable_default_field_dtype,
variable_default_datamodel, variable_default_distribution_strategy,
], ],
path=os.path.expanduser('~') + "/.nifty/nifty_config") path=os.path.expanduser('~') + "/.nifty/nifty_config")
######## ########
......
...@@ -6,7 +6,6 @@ from d2o import distributed_data_object,\ ...@@ -6,7 +6,6 @@ from d2o import distributed_data_object,\
from nifty.config import about,\ from nifty.config import about,\
nifty_configuration as gc,\ nifty_configuration as gc,\
dependency_injector as gdi
from nifty.field_types import FieldType from nifty.field_types import FieldType
...@@ -17,15 +16,11 @@ import nifty.nifty_utilities as utilities ...@@ -17,15 +16,11 @@ import nifty.nifty_utilities as utilities
from nifty.random import Random from nifty.random import Random
POINT_DISTRIBUTION_STRATEGIES = DISTRIBUTION_STRATEGIES['global']
COMM = getattr(gdi[gc['mpi_module']], gc['default_comm'])
class Field(object): class Field(object):
# ---Initialization methods--- # ---Initialization methods---
def __init__(self, domain=None, val=None, dtype=None, field_type=None, def __init__(self, domain=None, val=None, dtype=None, field_type=None,
datamodel=None, copy=False): distribution_strategy=None, copy=False):
self.domain = self._parse_domain(domain=domain, val=val) self.domain = self._parse_domain(domain=domain, val=val)
self.domain_axes = self._get_axes_tuple(self.domain) self.domain_axes = self._get_axes_tuple(self.domain)
...@@ -44,8 +39,9 @@ class Field(object): ...@@ -44,8 +39,9 @@ class Field(object):
domain=self.domain, domain=self.domain,
field_type=self.field_type) field_type=self.field_type)
self.datamodel = self._parse_datamodel(datamodel=datamodel, self.distribution_strategy = self._parse_distribution_strategy(
val=val) distribution_strategy=distribution_strategy,
val=val)
self.set_val(new_val=val, copy=copy) self.set_val(new_val=val, copy=copy)
...@@ -111,28 +107,29 @@ class Field(object): ...@@ -111,28 +107,29 @@ class Field(object):
return dtype return dtype
def _parse_datamodel(self, datamodel, val): def _parse_distribution_strategy(self, distribution_strategy, val):
if datamodel is None: if distribution_strategy is None:
if isinstance(val, distributed_data_object): if isinstance(val, distributed_data_object):
datamodel = val.distribution_strategy distribution_strategy = val.distribution_strategy
elif isinstance(val, Field): elif isinstance(val, Field):
datamodel = val.datamodel distribution_strategy = val.distribution_strategy
else: else:
about.warnings.cprint("WARNING: Datamodel set to default!") about.warnings.cprint("WARNING: Datamodel set to default!")
datamodel = gc['default_datamodel'] distribution_strategy = gc['default_distribution_strategy']
elif datamodel not in DISTRIBUTION_STRATEGIES['all']: elif distribution_strategy not in DISTRIBUTION_STRATEGIES['global']:
raise ValueError(about._errors.cstring( raise ValueError(about._errors.cstring(
"ERROR: Invalid datamodel!")) "ERROR: distribution_strategy must be a global-type "
return datamodel "strategy."))
return distribution_strategy
# ---Factory methods--- # ---Factory methods---
@classmethod @classmethod
def from_random(cls, random_type, domain=None, dtype=None, field_type=None, def from_random(cls, random_type, domain=None, dtype=None, field_type=None,
datamodel=None, **kwargs): distribution_strategy=None, **kwargs):
# create a initially empty field # create a initially empty field
f = cls(domain=domain, dtype=dtype, field_type=field_type, f = cls(domain=domain, dtype=dtype, field_type=field_type,
datamodel=datamodel) distribution_strategy=distribution_strategy)
# now use the processed input in terms of f in order to parse the # now use the processed input in terms of f in order to parse the
# random arguments # random arguments
...@@ -229,7 +226,7 @@ class Field(object): ...@@ -229,7 +226,7 @@ class Field(object):
harmonic_domain = self.domain[space_index] harmonic_domain = self.domain[space_index]
power_domain = PowerSpace(harmonic_domain=harmonic_domain, power_domain = PowerSpace(harmonic_domain=harmonic_domain,
datamodel=distribution_strategy, distribution_strategy=distribution_strategy,
log=log, nbin=nbin, binbounds=binbounds, log=log, nbin=nbin, binbounds=binbounds,
dtype=power_dtype) dtype=power_dtype)
...@@ -314,17 +311,116 @@ class Field(object): ...@@ -314,17 +311,116 @@ class Field(object):
return result_obj return result_obj
def power_synthesize(self, spaces=None, real_signal=True): def power_synthesize(self, spaces=None, real_signal=True):
# assert that all spaces in `self.domain` are eiher of signal-type or # assert that all spaces in `self.domain` are either of signal-type or
# power_space instances # power_space instances
for sp in self.domain: for sp in self.domain:
if sp.harmonic and not isinstance(sp, PowerSpace): if not sp.harmonic and not isinstance(sp, PowerSpace):
raise AttributeError( raise AttributeError(
"ERROR: Field has a space in `domain` which is neither " "ERROR: Field has a space in `domain` which is neither "
"harmonic nor a PowerSpace.") "harmonic nor a PowerSpace.")
pass
# synthesize random fields in harmonic domain using # check if the `spaces` input is valid
# np.random.multivariate_normal(mean=[0,0], cov=[[0.5,0],[0,0.5]], size=shape) spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
if len(self.domain) == 1:
spaces = (0,)
else:
raise ValueError(about._errors.cstring(
"ERROR: Field has multiple spaces as domain "
"but `spaces` is None."))
if len(spaces) == 0:
raise ValueError(about._errors.cstring(
"ERROR: No space for synthesis specified."))
elif len(spaces) > 1:
raise ValueError(about._errors.cstring(
"ERROR: Conversion of only one space at a time is allowed."))
power_space_index = spaces[0]
power_domain = self.domain[power_space_index]
if not isinstance(power_domain, PowerSpace):
raise ValueError(about._errors.cstring(
"ERROR: A PowerSpace is needed for field synthetization."))
# create the result domain
result_domain = list(self.domain)
harmonic_domain = power_domain.harmonic_domain
result_domain[power_space_index] = harmonic_domain
# create random samples: one or two, depending on whether the
# power spectrum is real or complex
if issubclass(power_domain.dtype.type, np.complexfloating):
result_list = [None, None]
else:
result_list = [None]
result_list = [self.__class__.from_random(
'normal',
result_domain,
dtype=harmonic_domain.dtype,
field_type=self.field_type,
distribution_strategy=self.distribution_strategy)
for x in result_list]
# from now on extract the values from the random fields for further
# processing without killing the fields.
# if the signal-space field should be real, hermitianize the field
# components
if real_signal:
result_val_list = [harmonic_domain.hermitian_decomposition(
x.val,
axes=x.domain_axes[power_space_index])[0]
for x in result_list]
else:
result_val_list = [x.val for x in result_list]
# weight the random fields with the power spectrum
# therefore get the pindex from the power space
pindex = power_domain.pindex
# take the local data from pindex. This data must be compatible to the
# local data of the field given the slice of the PowerSpace
local_distribution_strategy = \
result_list[0].val.get_axes_local_distribution_strategy(
result_list[0].domain_axes[power_space_index])
if pindex.distribution_strategy is not local_distribution_strategy:
about.warnings.cprint(
"WARNING: The distribution_stragey of pindex does not fit the "
"slice_local distribution strategy of the synthesized field.")
# Now use numpy advanced indexing in order to put the entries of the
# power spectrum into the appropriate places of the pindex array.
# Do this for every 'pindex-slice' in parallel using the 'slice(None)'s
local_pindex = pindex.get_local_data(copy=False)
local_spec = self.val.get_local_data(copy=False)
local_blow_up = [slice(None)]*len(self.shape)
local_blow_up[self.domain_axes[power_space_index][0]] = local_pindex
# here, the power_spectrum is distributed into the new shape
local_rescaler = local_spec[local_blow_up]
# apply the rescaler to the random fields
result_val_list[0].apply_scalar_function(
lambda x: x * local_rescaler.real,
inplace=True)
if issubclass(power_domain.dtype.type, np.complexfloating):
result_val_list[1].apply_scalar_function(
lambda x: x * local_rescaler.imag,
inplace=True)
# store the result into the fields
[x.set_val(new_val=y, copy=False) for x, y in
zip(result_list, result_val_list)]
if issubclass(power_domain.dtype.type, np.complexfloating):
result = result_list[0] + 1j*result_list[1]
else:
result = result_list[0]
return result
# ---Properties--- # ---Properties---
...@@ -413,24 +509,26 @@ class Field(object): ...@@ -413,24 +509,26 @@ class Field(object):
if dtype is None: if dtype is None:
dtype = self.dtype dtype = self.dtype
return_x = distributed_data_object(global_shape=self.shape, return_x = distributed_data_object(
dtype=dtype, global_shape=self.shape,
distribution_strategy=self.datamodel) dtype=dtype,
distribution_strategy=self.distribution_strategy)
return_x.set_full_data(x, copy=False) return_x.set_full_data(x, copy=False)
return return_x return return_x
def copy(self, domain=None, dtype=None, field_type=None, def copy(self, domain=None, dtype=None, field_type=None,
datamodel=None): distribution_strategy=None):
copied_val = self.get_val(copy=True) copied_val = self.get_val(copy=True)
new_field = self.copy_empty(domain=domain, new_field = self.copy_empty(
dtype=dtype, domain=domain,
field_type=field_type, dtype=dtype,
datamodel=datamodel) field_type=field_type,
distribution_strategy=distribution_strategy)
new_field.set_val(new_val=copied_val, copy=False) new_field.set_val(new_val=copied_val, copy=False)
return new_field return new_field
def copy_empty(self, domain=None, dtype=None, field_type=None, def copy_empty(self, domain=None, dtype=None, field_type=None,
datamodel=None): distribution_strategy=None):
if domain is None: if domain is None:
domain = self.domain domain = self.domain
else: else:
...@@ -446,8 +544,8 @@ class Field(object): ...@@ -446,8 +544,8 @@ class Field(object):
else: else:
field_type = self._parse_field_type(field_type) field_type = self._parse_field_type(field_type)
if datamodel is None: if distribution_strategy is None:
datamodel = self.datamodel distribution_strategy = self.distribution_strategy
fast_copyable = True fast_copyable = True
try: try:
...@@ -463,13 +561,13 @@ class Field(object): ...@@ -463,13 +561,13 @@ class Field(object):
fast_copyable = False fast_copyable = False
if (fast_copyable and dtype == self.dtype and if (fast_copyable and dtype == self.dtype and
datamodel == self.datamodel): distribution_strategy == self.distribution_strategy):
new_field = self._fast_copy_empty() new_field = self._fast_copy_empty()
else: else:
new_field = Field(domain=domain, new_field = Field(domain=domain,
dtype=dtype, dtype=dtype,
field_type=field_type, field_type=field_type,
datamodel=datamodel) distribution_strategy=distribution_strategy)
return new_field return new_field
def _fast_copy_empty(self): def _fast_copy_empty(self):
......
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from field_array import FieldType,\ from field_type import FieldType
FieldArray from field_array import FieldArray
\ No newline at end of file
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from base_field_type import FieldType from field_type import FieldType
class FieldArray(FieldType): class FieldArray(FieldType):
......
...@@ -16,26 +16,63 @@ class DiagonalOperator(EndomorphicOperator): ...@@ -16,26 +16,63 @@ class DiagonalOperator(EndomorphicOperator):
# ---Overwritten properties and methods--- # ---Overwritten properties and methods---
def __init__(self, domain=(), field_type=(), implemented=False, def __init__(self, domain=(), field_type=(), implemented=False,
diagonal=None, bare=False, copy=True, datamodel=None): diagonal=None, bare=False, copy=True,
distribution_strategy=None):
super(DiagonalOperator, self).__init__(domain=domain, super(DiagonalOperator, self).__init__(domain=domain,
field_type=field_type, field_type=field_type,
implemented=implemented) implemented=implemented)
self._implemented = bool(implemented) self._implemented = bool(implemented)
if datamodel is None: if distribution_strategy is None:
if isinstance(diagonal, distributed_data_object): if isinstance(diagonal, distributed_data_object):
datamodel = diagonal.distribution_strategy distribution_strategy = diagonal.distribution_strategy
elif isinstance(diagonal, Field): elif isinstance(diagonal, Field):
datamodel = diagonal.datamodel distribution_strategy = diagonal.distribution_strategy
self.datamodel = self._parse_datamodel(datamodel=datamodel, self.distribution_strategy = self._parse_distribution_strategy(
val=diagonal) distribution_strategy=distribution_strategy,
val=diagonal)
self.set_diagonal(diagonal=diagonal, bare=bare, copy=copy) self.set_diagonal(diagonal=diagonal, bare=bare, copy=copy)
def _times(self, x, spaces, types): def _times(self, x, spaces, types):
pass # if the distribution_strategy of self is sub-slice compatible to
# the one of x, reshape the local data of self and apply it directly
active_axes = []
if spaces is None:
for axes in x.domain_axes:
active_axes += axes
else:
for space_index in spaces:
active_axes += x.domain_axes[space_index]
if types is None:
for axes in x.field_type_axes:
active_axes += axes
else:
for type_index in types:
active_axes += x.field_type_axes[type_index]
if x.val.get_axes_local_distribution_strategy(active_axes) == \
self.distribution_strategy:
local_data = self._diagonal.val.get_local_data(copy=False)
# check if domains match completely
# -> multiply directly
# check if axes_local_distribution_strategy matches.
# If yes, extract local data of self.diagonal and x and use numpy
# reshape.
# assert that indices in spaces and types are striktly increasing
# otherwise a wild transpose would be necessary
# build new shape (1,1,x,1,y,1,1,z)
# copy self.diagonal into new shape
# apply reshaped array to x
def _adjoint_times(self, x, spaces, types): def _adjoint_times(self, x, spaces, types):
pass pass
...@@ -97,29 +134,29 @@ class DiagonalOperator(EndomorphicOperator): ...@@ -97,29 +134,29 @@ class DiagonalOperator(EndomorphicOperator):
# ---Added properties and methods--- # ---Added properties and methods---
@property @property
def datamodel(self): def distribution_strategy(self):
return self._datamodel return self._distribution_strategy
def _parse_datamodel(self, datamodel, val): def _parse_distribution_strategy(self, distribution_strategy, val):
if datamodel is None: if distribution_strategy is None:
if isinstance(val, distributed_data_object): if isinstance(val, distributed_data_object):
datamodel = val.distribution_strategy distribution_strategy = val.distribution_strategy
elif isinstance(val, Field): elif isinstance(val, Field):
datamodel = val.datamodel distribution_strategy = val.distribution_strategy
else: else:
about.warnings.cprint("WARNING: Datamodel set to default!") about.warnings.cprint("WARNING: Datamodel set to default!")
datamodel = gc['default_datamodel'] distribution_strategy = gc['default_distribution_strategy']
elif datamodel not in DISTRIBUTION_STRATEGIES['all']: elif distribution_strategy not in DISTRIBUTION_STRATEGIES['all']:
raise ValueError(about._errors.cstring( raise ValueError(about._errors.cstring(
"ERROR: Invalid datamodel!")) "ERROR: Invalid distribution_strategy!"))
return datamodel return distribution_strategy
def set_diagonal(self, diagonal, bare=False, copy=True): def set_diagonal(self, diagonal, bare=False, copy=True):
# use the casting functionality from Field to process `diagonal` # use the casting functionality from Field to process `diagonal`
f = Field(domain=self.domain, f = Field(domain=self.domain,
val=diagonal, val=diagonal,
field_type=self.field_type, field_type=self.field_type,
datamodel=self.datamodel, distribution_strategy=self.distribution_strategy,
copy=copy) copy=copy)
# weight if the given values were `bare` and `implemented` is True # weight if the given values were `bare` and `implemented` is True
......
...@@ -39,7 +39,10 @@ class FFTOperator(LinearOperator): ...@@ -39,7 +39,10 @@ class FFTOperator(LinearOperator):
def _times(self, x, spaces, types): def _times(self, x, spaces, types):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain)) spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
if spaces is None: if spaces is None:
axes = None # this case means that x lives on only one space, which is
# identical to the space in the domain of `self`. Otherwise the
# input check of LinearOperator would have failed.
axes = x.domain_axes[0]
else: else:
axes = x.domain_axes[spaces[0]] axes = x.domain_axes[spaces[0]]
...@@ -59,11 +62,13 @@ class FFTOperator(LinearOperator): ...@@ -59,11 +62,13 @@ class FFTOperator(LinearOperator):
def _inverse_times(self, x, spaces, types): def _inverse_times(self, x, spaces, types):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain)) spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
if spaces is None: if spaces is None:
axes = None # this case means that x lives on only one space, which is
# identical to the space in the domain of `self`. Otherwise the
# input check of LinearOperator would have failed.
axes = x.domain_axes[0]
else: else:
axes = x.domain_axes[spaces[0]] axes = x.domain_axes[spaces[0]]
axes = x.domain_axes[spaces[0]]
new_val = self._inverse_transformation.transform(x.val, axes=axes) new_val = self._inverse_transformation.transform(x.val, axes=axes)
if spaces is None: if spaces is None:
......
...@@ -157,8 +157,7 @@ class LinearOperator(object): ...@@ -157,8 +157,7 @@ class LinearOperator(object):
# cases: # cases:
# 1. Case: # 1. Case:
# The user specifies with `spaces` that the operators domain should # The user specifies with `spaces` that the operators domain should
# be applied to a certain domain in the domain-tuple of x. This is # be applied to certain spaces in the domain-tuple of x.
# only valid if len(self.domain)==1.
# 2. Case: # 2. Case:
# The domains of self and x match completely. # The domains of self and x match completely.
...@@ -175,16 +174,8 @@ class LinearOperator(object): ...@@ -175,16 +174,8 @@ class LinearOperator(object):
"ERROR: The operator's and and field's domains don't " "ERROR: The operator's and and field's domains don't "
"match.")) "match."))
else: else:
if len(self_domain) > 1: for i, space_index in enumerate(spaces):
raise ValueError(about._errors.cstring( if x.domain[space_index] != self_domain[i]:
"ERROR: Specifying `spaces` for operators with multiple "
"domain spaces is not valid."))
elif len(spaces) != len(self_domain):
raise ValueError(about._errors.cstring(
"ERROR: Length of `spaces` does not match the number of "
"spaces in the operator's domain."))
elif len(spaces) == 1:
if x.domain[spaces[0]] != self_domain[0]:
raise ValueError(about._errors.cstring( raise ValueError(about._errors.cstring(
"ERROR: The operator's and and field's domains don't " "ERROR: The operator's and and field's domains don't "
"match."))