Skip to content
GitLab
Commit 86b44136 authored by Theo Steininger's avatar Theo Steininger
Browse files

First steps towards multi power synthesize.

parent 7967f2d5
Pipeline #11764 passed with stage
in 10 minutes and 20 seconds
Hide whitespace changes
Inline Side-by-side
nifty/field.py
View file @ 86b44136
...@@ -303,59 +303,38 @@ class Field(Loggable, Versionable, object): ...@@ -303,59 +303,38 @@ 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: for power_space_index in spaces:
raise ValueError( power_space = self.domain[power_space_index]
"No space for synthesis specified.") if not isinstance(power_space, PowerSpace):
elif len(spaces) > 1: raise ValueError("A PowerSpace is needed for field "
raise ValueError( "synthetization.")
"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(
"A PowerSpace is needed for field 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:
result_domain[power_space_index] = harmonic_domain power_space = self.domain[power_space_index]
harmonic_domain = power_space.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(power_domain.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=harmonic_domain.dtype, dtype=np.complex,
distribution_strategy=self.distribution_strategy) distribution_strategy=self.distribution_strategy)
for x in result_list] for x in result_list]
...@@ -363,18 +342,54 @@ class Field(Loggable, Versionable, object): ...@@ -363,18 +342,54 @@ 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:
result_val_list = [harmonic_domain.hermitian_decomposition( for power_space_index in spaces:
x.val, harmonic_domain = result_domain[power_space_index]
axes=x.domain_axes[power_space_index], result_val_list = \
preserve_gaussian_variance=True)[0] [harmonic_domain.hermitian_decomposition(
for x in result_list] result_val_list[0],
axes=result_list[0].domain_axes[power_space_index],
preserve_gaussian_variance=True)[0],
harmonic_domain.hermitian_decomposition(
result_val_list[1],
axes=result_list[1].domain_axes[power_space_index],
preserve_gaussian_variance=True)[1]]
# 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 = \
...@@ -390,34 +405,12 @@ class Field(Loggable, Versionable, object): ...@@ -390,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(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---
......
nifty/plotting/figures/figure_2D.py
View file @ 86b44136
...@@ -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:
......
Supports Markdown
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