From 0c38434c0a3b31cfea8030e3877f7505dae25896 Mon Sep 17 00:00:00 2001 From: Marco Selig Date: Wed, 13 Mar 2013 13:31:37 +0100 Subject: [PATCH] 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. --- nifty_core.py | 61 ++++++++++++++--- nifty_power.py | 51 +++++++++++++- powerspectrum.py | 171 +++++++++++++++++++++++++++++------------------ 3 files changed, 207 insertions(+), 76 deletions(-) diff --git a/nifty_core.py b/nifty_core.py index 4cf168a5..cc27ec38 100644 --- a/nifty_core.py +++ b/nifty_core.py @@ -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] ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/nifty_power.py b/nifty_power.py index ea4126d3..58b3c3bf 100644 --- a/nifty_power.py +++ b/nifty_power.py @@ -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): diff --git a/powerspectrum.py b/powerspectrum.py index 823f1337..cbee5fd6 100644 --- a/powerspectrum.py +++ b/powerspectrum.py @@ -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 -- GitLab