Commit 0c38434c authored by Marco Selig's avatar Marco Selig

powerspectrum.calc_ps_fast optimized; nifty_core adjusted;...

powerspectrum.calc_ps_fast optimized; nifty_core adjusted; power_operator.get_power optimized by introduction of a power un(in)dex; weight_power added to nifty_power; docstrings revised.
parent b4477a95
......@@ -1036,6 +1036,35 @@ class space(object):
"""
raise AttributeError(about._errors.cstring("ERROR: no generic instance method 'get_power_index'."))
def get_power_undex(self,pindex=None):
"""
Provides the unindexing list for an indexed power spectrum.
Parameters
----------
pindex : numpy.ndarray, *optional*
Indexing array giving the power spectrum index for each
represented mode.
Returns
-------
pundex : list
Unindexing list undoing power indexing.
Notes
-----
Indexing with the unindexing list undoes the indexing with the
indexing array; i.e., ``x == x[pindex][pundex]``.
See also
--------
get_power_index
"""
if(pindex is None):
pindex = self.get_power_index(irreducible=False)
return list(np.unravel_index(np.unique(pindex,return_index=True,return_inverse=False)[1],pindex.shape,order='C'))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def enforce_shape(self,x):
......@@ -1379,6 +1408,9 @@ class space(object):
pindex : numpy.ndarray, *optional*
Indexing array assigning the input array components to
components of the power spectrum (default: None).
kindex : numpy.ndarray, *optional*
Scale corresponding to each band in the power spectrum
(default: None).
rho : numpy.ndarray, *optional*
Number of degrees of freedom per band (default: None).
"""
......@@ -2607,9 +2639,12 @@ class rg_space(space):
Other parameters
----------------
ind : numpy.ndarray, *optional*
pindex : numpy.ndarray, *optional*
Indexing array assigning the input array components to
components of the power spectrum (default: None).
kindex : numpy.ndarray, *optional*
Scale corresponding to each band in the power spectrum
(default: None).
rho : numpy.ndarray, *optional*
Number of degrees of freedom per band (default: None).
"""
......@@ -2618,10 +2653,11 @@ class rg_space(space):
if(not self.fourier):
x = self.calc_weight(x,power=1)
## power spectrum
if(kwargs.has_key("ind"))and(kwargs.has_key("rho")):
return gp.calc_ps_fast(x,self.para[:(np.size(self.para)-1)//2],self.vol,[kwargs.get("ind"),kwargs.get("rho")],self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=self.fourier)
else:
return gp.calc_ps(x,self.para[:(np.size(self.para)-1)//2],self.vol,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=self.fourier)
return gp.calc_ps_fast(x,self.para[:(np.size(self.para)-1)//2],self.vol,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=self.fourier,**kwargs)
# if(kwargs.has_key("pindex"))and(kwargs.has_key("rho")):
# return gp.calc_ps_fast(x,self.para[:(np.size(self.para)-1)//2],self.vol,[kwargs.get("pindex"),kwargs.get("rho")],self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=self.fourier)
# else:
# return gp.calc_ps(x,self.para[:(np.size(self.para)-1)//2],self.vol,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=self.fourier)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -6075,13 +6111,15 @@ class field(object):
pindex : ndarray
Specifies the indexing array for the distribution of
indices in conjugate space (default: None).
kindex : numpy.ndarray, *optional*
Scale corresponding to each band in the power spectrum
(default: None).
rho : scalar
Number of degrees of freedom per irreducible band
(default=None).
iter : scalar
Number of iterations (default: 0)
Returns
-------
spec : ndarray
......@@ -8638,7 +8676,7 @@ class power_operator(diagonal_operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_power(self,bare=True,**kwargs):
def get_power(self,bare=True,pundex=None,pindex=None,**kwargs): ## **kwargs for downward compatibility
"""
Computes the power spectrum
......@@ -8648,11 +8686,12 @@ class power_operator(diagonal_operator):
whether the entries are `bare` or not
(mandatory for the correct incorporation of volume weights)
(default: False)
pundex : ndarray, *optional*
unindexing array, obtainable from domain.get_power_undex
(default: None)
pindex : ndarray, *optional*
indexing array, obtainable from domain.get_power_index
(default: None)
rho : ndarray, *optional*
number of degrees of freedom per irreducible band (default: None)
Returns
-------
......@@ -8665,7 +8704,9 @@ class power_operator(diagonal_operator):
else:
diag = self.val
return self.domain.calc_power(np.sqrt(diag),**kwargs) ## TODO: check efficiency
if(pundex is None):
pundex = self.domain.get_power_undex(pindex=pindex)
return diag[pundex]
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......
......@@ -45,10 +45,59 @@
"""
from __future__ import division
import numpy as np
from nifty_core import *
#import numpy as np
import smoothing as gs
##-----------------------------------------------------------------------------
def weight_power(domain,spec,power=1,pindex=None,pundex=None):
"""
Weights a given power spectrum with the corresponding pixel volumes
to a given power.
Parameters
----------
domain : space
The space wherein valid arguments live.
spec : {scalar, ndarray, field}
The power spectrum. A scalars is interpreted as a constant
spectrum.
pindex : ndarray
Indexing array giving the power spectrum index for each
represented mode.
pundex : list
Unindexing list undoing power indexing.
Returns
-------
spev : ndarray
Weighted power spectrum.
Raises
------
TypeError
If `domain` is no space.
ValueError
If `domain` is no harmonic space.
"""
## check domain
if(not isinstance(domain,space)):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
if(pindex is None):
try:
pindex = domain.get_power_index(irreducible=False)
except(AttributeError):
raise ValueError(about._errors.cstring("ERROR: invalid input."))
if(pundex is None):
pundex = domain.get_power_undex(pindex=pindex)
return np.real(domain.calc_weight(domain.enforce_power(spec,size=len(set(pindex.flatten(order='C'))))[pindex],power=power)[pundex])
##-----------------------------------------------------------------------------
def smooth_power(power,kindex,exclude=1,sigma=-1):
......
......@@ -123,57 +123,58 @@ def draw_vector_nd(axes,dgrid,ps,symtype=0,fourier=False,zerocentered=False):
return vector
def calc_ps(field,axes,dgrid,zerocentered=False,fourier=False):
"""
Calculates the power spectrum of a given field assuming that the field
is statistically homogenous and isotropic.
Parameters
----------
field : ndarray
The input field from which the power spectrum should be determined.
axes : ndarray
An array with the length of each axis.
dgrid : ndarray
An array with the pixel length of each axis.
zerocentered : bool : *optional*
Whether the output array should be zerocentered, i.e. starting with
negative Fourier modes going over the zero mode to positive modes,
or not zerocentered, where zero, positive and negative modes are
simpy ordered consecutively.
fourier : bool : *optional*
Whether the output should be in Fourier space or not
(default=False).
"""
## field absolutes
if(not fourier):
foufield = np.fft.fftshift(np.fft.fftn(field))
elif(np.any(zerocentered==False)):
foufield = np.fft.fftshift(field, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
else:
foufield = field
fieldabs = np.abs(foufield)**2
kdict = nkdict(axes,dgrid,fourier)
klength = nklength(kdict)
ps = np.zeros(klength.size)
rho = np.zeros(klength.size)
for ii in np.ndindex(kdict.shape):
position = np.searchsorted(klength,kdict[ii])
rho[position] += 1
ps[position] += fieldabs[ii]
ps = np.divide(ps,rho)
return ps
def calc_ps_fast(field,axes,dgrid,kpack,zerocentered=False,fourier=False): ## kpack = [powerind,rho]
#def calc_ps(field,axes,dgrid,zerocentered=False,fourier=False):
#
# """
# Calculates the power spectrum of a given field assuming that the field
# is statistically homogenous and isotropic.
#
# Parameters
# ----------
# field : ndarray
# The input field from which the power spectrum should be determined.
#
# axes : ndarray
# An array with the length of each axis.
#
# dgrid : ndarray
# An array with the pixel length of each axis.
#
# zerocentered : bool : *optional*
# Whether the output array should be zerocentered, i.e. starting with
# negative Fourier modes going over the zero mode to positive modes,
# or not zerocentered, where zero, positive and negative modes are
# simpy ordered consecutively.
#
# fourier : bool : *optional*
# Whether the output should be in Fourier space or not
# (default=False).
#
# """
#
# ## field absolutes
# if(not fourier):
# foufield = np.fft.fftshift(np.fft.fftn(field))
# elif(np.any(zerocentered==False)):
# foufield = np.fft.fftshift(field, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
# else:
# foufield = field
# fieldabs = np.abs(foufield)**2
#
# kdict = nkdict(axes,dgrid,fourier)
# klength = nklength(kdict)
#
# ## power spectrum
# ps = np.zeros(klength.size)
# rho = np.zeros(klength.size)
# for ii in np.ndindex(kdict.shape):
# position = np.searchsorted(klength,kdict[ii])
# rho[position] += 1
# ps[position] += fieldabs[ii]
# ps = np.divide(ps,rho)
# return ps
def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,kindex=None,rho=None):
"""
Calculates the power spectrum of a given field faster assuming that the
......@@ -190,15 +191,6 @@ def calc_ps_fast(field,axes,dgrid,kpack,zerocentered=False,fourier=False): ## kp
dgrid : ndarray
An array with the pixel length of each axis.
kpack : list {powerind,rho}
Specifies pre-calculated properties of the Fourier grid.
- powerind gives the index of the Fourier grid points in a numpy
array ordered following the zerocentered flag.
- rho is the degeneracy of the Fourier grid, indicating how many
k-vectors in Fourier space have the same length.
zerocentered : bool : *optional*
Whether the output array should be zerocentered, i.e. starting with
negative Fourier modes going over the zero mode to positive modes,
......@@ -209,6 +201,16 @@ def calc_ps_fast(field,axes,dgrid,kpack,zerocentered=False,fourier=False): ## kp
Whether the output should be in Fourier space or not
(default=False).
pindex : ndarray
Index of the Fourier grid points in a numpy.ndarray ordered
following the zerocentered flag (default=None).
kindex : ndarray
Array of all k-vector lengths (default=None).
rho : ndarray
Degeneracy of the Fourier grid, indicating how many k-vectors in
Fourier space have the same length (default=None).
"""
......@@ -220,11 +222,50 @@ def calc_ps_fast(field,axes,dgrid,kpack,zerocentered=False,fourier=False): ## kp
else:
foufield = field
fieldabs = np.abs(foufield)**2
## power spectrum
ps = np.zeros(kpack[1].size)
for ii in np.ndindex(kpack[0].shape):
ps[kpack[0][ii]] += fieldabs[ii]
ps = np.divide(ps,kpack[1])
if(rho is None):
if(pindex is None):
## kdict
kdict = nkdict(axes,dgrid,fourier)
## klength
if(kindex is None):
klength = nklength(kdict)
else:
klength = kindex
## power spectrum
ps = np.zeros(klength.size)
rho = np.zeros(klength.size)
for ii in np.ndindex(kdict.shape):
position = np.searchsorted(klength,kdict[ii])
ps[position] += fieldabs[ii]
rho[position] += 1
else:
## power spectrum
ps = np.zeros(len(set(pindex.flatten())))
rho = np.zeros(ps.size)
for ii in np.ndindex(pindex.shape):
ps[pindex[ii]] += fieldabs[ii]
rho[pindex[ii]] += 1
elif(pindex is None):
## kdict
kdict = nkdict(axes,dgrid,fourier)
## klength
if(kindex is None):
klength = nklength(kdict)
else:
klength = kindex
## power spectrum
ps = np.zeros(klength.size)
for ii in np.ndindex(kdict.shape):
position = np.searchsorted(klength,kdict[ii])
ps[position] += fieldabs[ii]
else:
## power spectrum
ps = np.zeros(rho.size)
for ii in np.ndindex(pindex.shape):
ps[pindex[ii]] += fieldabs[ii]
ps = np.divide(ps,rho)
return ps
......
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