Commit 42951308 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'multi_power_synth' into 'master'

Multi power synth

See merge request !78
parents 9ba59e8c 1f5a6c6c
Pipeline #11887 passed with stages
in 16 minutes and 6 seconds
...@@ -90,9 +90,11 @@ class Field(Loggable, Versionable, object): ...@@ -90,9 +90,11 @@ class Field(Loggable, Versionable, object):
try: try:
dtype = val.dtype dtype = val.dtype
except AttributeError: except AttributeError:
if val is not None: try:
if val is None:
raise TypeError
dtype = np.result_type(val) dtype = np.result_type(val)
else: except(TypeError):
dtype = np.dtype(gc['default_field_dtype']) dtype = np.dtype(gc['default_field_dtype'])
else: else:
dtype = np.dtype(dtype) dtype = np.dtype(dtype)
...@@ -301,59 +303,41 @@ class Field(Loggable, Versionable, object): ...@@ -301,59 +303,41 @@ class Field(Loggable, Versionable, object):
return result_obj return result_obj
def power_synthesize(self, spaces=None, real_signal=True, def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
mean=None, std=None): mean=None, std=None):
# check if all spaces in `self.domain` are either of signal-type or
# power_space instances
for sp in self.domain:
if not sp.harmonic and not isinstance(sp, PowerSpace):
self.logger.info(
"Field has a space in `domain` which is neither "
"harmonic nor a PowerSpace.")
# check if the `spaces` input is valid # check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain)) spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
if len(self.domain) == 1:
spaces = (0,)
else:
raise ValueError(
"Field has multiple spaces as domain "
"but `spaces` is None.")
if len(spaces) == 0: if spaces is None:
raise ValueError( spaces = range(len(self.domain))
"No space for synthesis specified.")
elif len(spaces) > 1:
raise ValueError(
"Conversion of only one space at a time is allowed.")
power_space_index = spaces[0] for power_space_index in spaces:
power_domain = self.domain[power_space_index] power_space = self.domain[power_space_index]
if not isinstance(power_domain, PowerSpace): if not isinstance(power_space, PowerSpace):
raise ValueError( raise ValueError("A PowerSpace is needed for field "
"A PowerSpace is needed for field synthetization.") "synthetization.")
# create the result domain # create the result domain
result_domain = list(self.domain) result_domain = list(self.domain)
harmonic_domain = power_domain.harmonic_domain for power_space_index in spaces:
power_space = self.domain[power_space_index]
harmonic_domain = power_space.harmonic_domain
result_domain[power_space_index] = harmonic_domain result_domain[power_space_index] = harmonic_domain
# create random samples: one or two, depending on whether the # create random samples: one or two, depending on whether the
# power spectrum is real or complex # power spectrum is real or complex
if real_power:
if issubclass(self.dtype.type, np.complexfloating):
result_list = [None, None]
else:
result_list = [None] result_list = [None]
else:
result_list = [None, None]
result_list = [self.__class__.from_random( result_list = [self.__class__.from_random(
'normal', 'normal',
mean=mean, mean=mean,
std=std, std=std,
domain=result_domain, domain=result_domain,
dtype=self.dtype, dtype=np.complex,
distribution_strategy=self.distribution_strategy) distribution_strategy=self.distribution_strategy)
for x in result_list] for x in result_list]
...@@ -361,18 +345,51 @@ class Field(Loggable, Versionable, object): ...@@ -361,18 +345,51 @@ class Field(Loggable, Versionable, object):
# processing without killing the fields. # processing without killing the fields.
# if the signal-space field should be real, hermitianize the field # if the signal-space field should be real, hermitianize the field
# components # components
spec = self.val.get_full_data()
for power_space_index in spaces:
spec = self._spec_to_rescaler(spec, result_list, power_space_index)
local_rescaler = spec
result_val_list = [x.val for x in result_list]
# apply the rescaler to the random fields
result_val_list[0].apply_scalar_function(
lambda x: x * local_rescaler.real,
inplace=True)
if not real_power:
result_val_list[1].apply_scalar_function(
lambda x: x * local_rescaler.imag,
inplace=True)
if real_signal: if real_signal:
for power_space_index in spaces:
harmonic_domain = result_domain[power_space_index]
result_val_list = [harmonic_domain.hermitian_decomposition( result_val_list = [harmonic_domain.hermitian_decomposition(
x.val, result_val,
axes=x.domain_axes[power_space_index], axes=result.domain_axes[power_space_index],
preserve_gaussian_variance=True)[0] preserve_gaussian_variance=True)[0]
for x in result_list] for (result, result_val)
in zip(result_list, result_val_list)]
# 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 real_power:
result = result_list[0]
else: else:
result_val_list = [x.val for x in result_list] result = result_list[0] + 1j*result_list[1]
return result
def _spec_to_rescaler(self, spec, result_list, power_space_index):
power_space = self.domain[power_space_index]
# weight the random fields with the power spectrum # weight the random fields with the power spectrum
# therefore get the pindex from the power space # therefore get the pindex from the power space
pindex = power_domain.pindex pindex = power_space.pindex
# take the local data from pindex. This data must be compatible to the # 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 data of the field given the slice of the PowerSpace
local_distribution_strategy = \ local_distribution_strategy = \
...@@ -388,34 +405,12 @@ class Field(Loggable, Versionable, object): ...@@ -388,34 +405,12 @@ class Field(Loggable, Versionable, object):
# power spectrum into the appropriate places of the pindex array. # power spectrum into the appropriate places of the pindex array.
# Do this for every 'pindex-slice' in parallel using the 'slice(None)'s # Do this for every 'pindex-slice' in parallel using the 'slice(None)'s
local_pindex = pindex.get_local_data(copy=False) local_pindex = pindex.get_local_data(copy=False)
full_spec = self.val.get_full_data()
local_blow_up = [slice(None)]*len(self.shape) local_blow_up = [slice(None)]*len(self.shape)
local_blow_up[self.domain_axes[power_space_index][0]] = local_pindex local_blow_up[self.domain_axes[power_space_index][0]] = local_pindex
# here, the power_spectrum is distributed into the new shape # here, the power_spectrum is distributed into the new shape
local_rescaler = full_spec[local_blow_up] local_rescaler = spec[local_blow_up]
return local_rescaler
# apply the rescaler to the random fields
result_val_list[0].apply_scalar_function(
lambda x: x * local_rescaler.real,
inplace=True)
if issubclass(self.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(self.dtype.type, np.complexfloating):
result = result_list[0] + 1j*result_list[1]
else:
result = result_list[0]
return result
# ---Properties--- # ---Properties---
......
...@@ -10,7 +10,7 @@ class Figure2D(FigureFromPlot): ...@@ -10,7 +10,7 @@ class Figure2D(FigureFromPlot):
# TODO: add sanitization of plots input # TODO: add sanitization of plots input
if isinstance(plots[0], Heatmap) and not width and not height: if isinstance(plots[0], Heatmap) and not width and not height:
(x, y) = self.plots[0].data.shape (x, y) = plots[0].data.shape
if x > y: if x > y:
width = 500 width = 500
...@@ -18,7 +18,7 @@ class Figure2D(FigureFromPlot): ...@@ -18,7 +18,7 @@ class Figure2D(FigureFromPlot):
else: else:
height = 500 height = 500
width = int(500 * y / x) width = int(500 * y / x)
if isinstance(self.plots[0], Mollweide): if isinstance(plots[0], Mollweide):
if not xaxis: if not xaxis:
xaxis = False xaxis = False
if not yaxis: if not yaxis:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment