Commit 4a64eb39 authored by Theo Steininger's avatar Theo Steininger
Browse files

Changed convention from `sqrt(power_spectrum)` to `power_spectrum`.

parent ccacfc73
Pipeline #13278 passed with stage
in 6 minutes and 39 seconds
...@@ -283,7 +283,6 @@ class Field(Loggable, Versionable, object): ...@@ -283,7 +283,6 @@ class Field(Loggable, Versionable, object):
spaces : int *optional* spaces : int *optional*
The subspace for which the powerspectrum shall be computed The subspace for which the powerspectrum shall be computed
(default : None). (default : None).
if spaces==None : Tries to synthesize for the whole domain
logarithmic : boolean *optional* logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning. True if the output PowerSpace should use logarithmic binning.
{default : False} {default : False}
...@@ -336,23 +335,31 @@ class Field(Loggable, Versionable, object): ...@@ -336,23 +335,31 @@ class Field(Loggable, Versionable, object):
# 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 spaces is None:
if len(self.domain) == 1: spaces = range(len(self.domain))
spaces = (0,)
else:
raise ValueError(
"Field has multiple spaces as domain "
"but `spaces` is None.")
if len(spaces) == 0: if len(spaces) == 0:
raise ValueError( raise ValueError(
"No space for analysis specified.") "No space for analysis specified.")
elif len(spaces) > 1:
raise ValueError(
"Conversion of only one space at a time is allowed.")
space_index = spaces[0] work_field = abs(self)
work_field = work_field*work_field
for space_index in spaces:
work_field = self._single_power_analyze(
work_field=work_field,
space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds,
keep_phase_information=keep_phase_information)
return work_field
@classmethod
def _single_power_analyze(cls, work_field, space_index, logarithmic, nbin,
binbounds, keep_phase_information):
if not self.domain[space_index].harmonic: if not work_field.domain[space_index].harmonic:
raise ValueError( raise ValueError(
"The analyzed space must be harmonic.") "The analyzed space must be harmonic.")
...@@ -363,10 +370,10 @@ class Field(Loggable, Versionable, object): ...@@ -363,10 +370,10 @@ class Field(Loggable, Versionable, object):
# If it was complex, all the power is put into a real power spectrum. # If it was complex, all the power is put into a real power spectrum.
distribution_strategy = \ distribution_strategy = \
self.val.get_axes_local_distribution_strategy( work_field.val.get_axes_local_distribution_strategy(
self.domain_axes[space_index]) work_field.domain_axes[space_index])
harmonic_domain = self.domain[space_index] harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain, power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy, distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin, logarithmic=logarithmic, nbin=nbin,
...@@ -379,32 +386,32 @@ class Field(Loggable, Versionable, object): ...@@ -379,32 +386,32 @@ class Field(Loggable, Versionable, object):
if keep_phase_information: if keep_phase_information:
hermitian_part, anti_hermitian_part = \ hermitian_part, anti_hermitian_part = \
harmonic_domain.hermitian_decomposition( harmonic_domain.hermitian_decomposition(
self.val, work_field.val,
axes=self.domain_axes[space_index]) axes=work_field.domain_axes[space_index])
[hermitian_power, anti_hermitian_power] = \ [hermitian_power, anti_hermitian_power] = \
[self._calculate_power_spectrum( [cls._calculate_power_spectrum(
x=part, field_val=part,
pindex=pindex, pindex=pindex,
rho=rho, rho=rho,
axes=self.domain_axes[space_index]) axes=work_field.domain_axes[space_index])
for part in [hermitian_part, anti_hermitian_part]] for part in [hermitian_part, anti_hermitian_part]]
power_spectrum = hermitian_power + 1j * anti_hermitian_power power_spectrum = hermitian_power + 1j * anti_hermitian_power
else: else:
power_spectrum = self._calculate_power_spectrum( power_spectrum = cls._calculate_power_spectrum(
x=self.val, field_val=work_field.val,
pindex=pindex, pindex=pindex,
rho=rho, rho=rho,
axes=self.domain_axes[space_index]) axes=work_field.domain_axes[space_index])
# create the result field and put power_spectrum into it # create the result field and put power_spectrum into it
result_domain = list(self.domain) result_domain = list(work_field.domain)
result_domain[space_index] = power_domain result_domain[space_index] = power_domain
result_dtype = power_spectrum.dtype result_dtype = power_spectrum.dtype
result_field = self.copy_empty( result_field = work_field.copy_empty(
domain=result_domain, domain=result_domain,
dtype=result_dtype, dtype=result_dtype,
distribution_strategy=power_spectrum.distribution_strategy) distribution_strategy=power_spectrum.distribution_strategy)
...@@ -412,17 +419,16 @@ class Field(Loggable, Versionable, object): ...@@ -412,17 +419,16 @@ class Field(Loggable, Versionable, object):
return result_field return result_field
def _calculate_power_spectrum(self, x, pindex, rho, axes=None): @classmethod
fieldabs = abs(x) def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None):
fieldabs **= 2
if axes is not None: if axes is not None:
pindex = self._shape_up_pindex( pindex = cls._shape_up_pindex(
pindex=pindex, pindex=pindex,
target_shape=x.shape, target_shape=field_val.shape,
target_strategy=x.distribution_strategy, target_strategy=field_val.distribution_strategy,
axes=axes) axes=axes)
power_spectrum = pindex.bincount(weights=fieldabs, power_spectrum = pindex.bincount(weights=field_val,
axis=axes) axis=axes)
if axes is not None: if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape) new_rho_shape = [1, ] * len(power_spectrum.shape)
...@@ -430,10 +436,10 @@ class Field(Loggable, Versionable, object): ...@@ -430,10 +436,10 @@ class Field(Loggable, Versionable, object):
rho = rho.reshape(new_rho_shape) rho = rho.reshape(new_rho_shape)
power_spectrum /= rho power_spectrum /= rho
power_spectrum **= 0.5
return power_spectrum return power_spectrum
def _shape_up_pindex(self, pindex, target_shape, target_strategy, axes): @staticmethod
def _shape_up_pindex(pindex, target_shape, target_strategy, axes):
if pindex.distribution_strategy not in \ if pindex.distribution_strategy not in \
DISTRIBUTION_STRATEGIES['global']: DISTRIBUTION_STRATEGIES['global']:
raise ValueError("pindex's distribution strategy must be " raise ValueError("pindex's distribution strategy must be "
...@@ -548,6 +554,8 @@ class Field(Loggable, Versionable, object): ...@@ -548,6 +554,8 @@ class Field(Loggable, Versionable, object):
# components # components
spec = self.val.get_full_data() spec = self.val.get_full_data()
spec = np.sqrt(spec)
for power_space_index in spaces: for power_space_index in spaces:
spec = self._spec_to_rescaler(spec, result_list, power_space_index) spec = self._spec_to_rescaler(spec, result_list, power_space_index)
local_rescaler = spec local_rescaler = spec
...@@ -1147,8 +1155,8 @@ class Field(Loggable, Versionable, object): ...@@ -1147,8 +1155,8 @@ class Field(Loggable, Versionable, object):
""" """
return_field = self.copy_empty()
new_val = abs(self.get_val(copy=False)) new_val = abs(self.get_val(copy=False))
return_field = self.copy_empty(dtype=new_val.dtype)
return_field.set_val(new_val, copy=False) return_field.set_val(new_val, copy=False)
return return_field return return_field
......
...@@ -32,7 +32,7 @@ def create_power_operator(domain, power_spectrum, dtype=None, ...@@ -32,7 +32,7 @@ def create_power_operator(domain, power_spectrum, dtype=None,
Parameters Parameters
---------- ----------
domain : DomainObject domain : DomainObject
Domain over which the power operator shall live. Domain over which the power operator shall live.
power_spectrum : (array-like, method) power_spectrum : (array-like, method)
An array-like object, or a method that implements the square root An array-like object, or a method that implements the square root
of a power spectrum as a function of k. of a power spectrum as a function of k.
...@@ -54,7 +54,6 @@ def create_power_operator(domain, power_spectrum, dtype=None, ...@@ -54,7 +54,6 @@ def create_power_operator(domain, power_spectrum, dtype=None,
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
fp = Field(power_domain, val=power_spectrum, dtype=dtype, fp = Field(power_domain, val=power_spectrum, dtype=dtype,
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
fp **= 2
f = fp.power_synthesize(mean=1, std=0, real_signal=False) f = fp.power_synthesize(mean=1, std=0, real_signal=False)
return DiagonalOperator(domain, diagonal=f, bare=True) return DiagonalOperator(domain, diagonal=f, bare=True)
...@@ -32,20 +32,20 @@ from itertools import product, chain ...@@ -32,20 +32,20 @@ from itertools import product, chain
from d2o.config import dependency_injector as gdi from d2o.config import dependency_injector as gdi
HARMONIC_SPACES = [RGSpace((8,), harmonic=True), HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
RGSpace((7,), harmonic=True,zerocenter=True), RGSpace((7,), harmonic=True,zerocenter=True),
RGSpace((8,), harmonic=True,zerocenter=True), RGSpace((8,), harmonic=True,zerocenter=True),
RGSpace((7,8), harmonic=True), RGSpace((7,8), harmonic=True),
RGSpace((7,8), harmonic=True, zerocenter=True), RGSpace((7,8), harmonic=True, zerocenter=True),
RGSpace((6,6), harmonic=True, zerocenter=True), RGSpace((6,6), harmonic=True, zerocenter=True),
RGSpace((7,5), harmonic=True, zerocenter=True), RGSpace((7,5), harmonic=True, zerocenter=True),
RGSpace((5,5), harmonic=True), RGSpace((5,5), harmonic=True),
RGSpace((4,5,7), harmonic=True), RGSpace((4,5,7), harmonic=True),
RGSpace((4,5,7), harmonic=True, zerocenter=True), RGSpace((4,5,7), harmonic=True, zerocenter=True),
LMSpace(6), LMSpace(6),
LMSpace(9)] LMSpace(9)]
#Try all sensible kinds of combinations of spaces, distributuion strategy and #Try all sensible kinds of combinations of spaces, distributuion strategy and
#binning parameters #binning parameters
_maybe_fftw = ["fftw"] if ('pyfftw' in gdi) else [] _maybe_fftw = ["fftw"] if ('pyfftw' in gdi) else []
...@@ -138,13 +138,13 @@ class PowerSpaceConsistencyCheck(unittest.TestCase): ...@@ -138,13 +138,13 @@ class PowerSpaceConsistencyCheck(unittest.TestCase):
binbounds=binbounds) binbounds=binbounds)
assert_equal(p.pindex.flatten()[p.pundex],np.arange(p.dim), assert_equal(p.pindex.flatten()[p.pundex],np.arange(p.dim),
err_msg='pundex is not right-inverse of pindex!') err_msg='pundex is not right-inverse of pindex!')
@expand(CONSISTENCY_CONFIGS) @expand(CONSISTENCY_CONFIGS)
def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy, def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy,
binbounds, nbin,logarithmic): binbounds, nbin,logarithmic):
assert_equal(p.pindex.flatten().bincount(), p.rho, assert_equal(p.pindex.flatten().bincount(), p.rho,
err_msg='rho is not equal to pindex degeneracy') err_msg='rho is not equal to pindex degeneracy')
class PowerSpaceFunctionalityTest(unittest.TestCase): class PowerSpaceFunctionalityTest(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS) @expand(CONSISTENCY_CONFIGS)
def test_constructor(self, harmonic_partner, distribution_strategy, def test_constructor(self, harmonic_partner, distribution_strategy,
......
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