There is a maintenance of MPCDF Gitlab on Thursday, April 22st 2020, 9:00 am CEST - Expect some service interruptions during this time

Commit 1cad6829 authored by theos's avatar theos

RGPowerSpace is now able to compute power spectra.

parent d380119a
......@@ -220,3 +220,20 @@ def field_map(ishape, function, *args):
)
# result[ii] = function(*map(lambda z: z[ii], args))
return result
def cast_axis_to_tuple(axis, length):
if axis is None:
return None
try:
axis = tuple(int(item) for item in axis)
except(TypeError):
if np.isscalar(axis):
axis = (int(axis), )
else:
raise TypeError(
"ERROR: Could not convert axis-input to tuple of ints")
# shift negative indices to positive ones
axis = tuple(item if (item >= 0) else (item + length) for item in axis)
return axis
\ No newline at end of file
......@@ -563,7 +563,7 @@ class RGPowerIndices(PowerIndices):
#######
# rho #
#######
global_rho = global_pindex.bincount()
global_rho = global_pindex.bincount().get_full_data()
##########
# pundex #
......
# -*- coding: utf-8 -*-
import numpy as np
from d2o import STRATEGIES
from nifty.power import PowerSpace
from nifty.nifty_paradict import rg_power_space_paradict
from power_index_factory import RGPowerIndexFactory
......@@ -30,6 +32,51 @@ class RGPowerSpace(PowerSpace):
self.harmonic = True
self.discrete = False
@property
def shape(self):
return tuple(self.paradict['shape'])
def calculate_power_spectrum(self, x, axes=None):
fieldabs = abs(x)**2
# need a bincount with axes function here.
pindex = self.power_indices['pindex']
if axes is not None:
pindex = self._shape_up_pindex(
pindex=pindex,
target_shape=x.shape,
target_strategy=x.distribution_strategy,
axes=axes)
power_spectrum = pindex.bincount(weights=fieldabs,
axis=axes)
rho = self.power_indices['rho']
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
new_rho_shape[axes[0]] = len(rho)
rho = rho.reshape(new_rho_shape)
power_spectrum /= rho
return power_spectrum
def _shape_up_pindex(self, pindex, target_shape, target_strategy, axes):
if pindex.distribution_strategy not in STRATEGIES['global']:
raise ValueError("ERROR: pindex's distribution strategy must be "
"global-type")
if pindex.distribution_strategy in STRATEGIES['slicing']:
if ((0 not in axes) or
(target_strategy is not pindex.distribution_strategy)):
raise ValueError(
"ERROR: A slicing distributor shall not be reshaped to "
"something non-sliced.")
semiscaled_shape = [1, ] * len(target_shape)
for i in axes:
semiscaled_shape[i] = target_shape[i]
local_data = pindex.get_local_data(copy=False)
semiscaled_local_data = local_data.reshape(semiscaled_shape)
result_obj = pindex.copy_empty(global_shape=target_shape,
distribution_strategy=target_strategy)
result_obj.set_full_data(semiscaled_local_data, copy=False)
return result_obj
......@@ -954,7 +954,7 @@ class RGSpace(Space):
rho = kwargs.get('rho', power_indices['rho'])
fieldabs = abs(x)**2
power_spectrum = np.zeros(rho.shape)
#power_spectrum = np.zeros(rho.shape)
power_spectrum = pindex.bincount(weights=fieldabs)
......
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