Commit 1300fa47 authored by Marco Selig's avatar Marco Selig
Browse files

get_power_indices sketched; gfft dependence from smoothing removed.

parent 3c7c1a1f
...@@ -1080,76 +1080,13 @@ class space(object): ...@@ -1080,76 +1080,13 @@ class space(object):
""" """
""" """
pass pindex = self.get_power_index(irreducible=False)
kindex,rho = self.get_power_index(irreducible=True)
# factor = 1
# lnfactor = np.log(factor) # pindex,kindex,rho = gp.reset_power_indices(pindex,kindex,rho)
#
# pindex = self.get_power_index(irreducible=False) pundex = list(np.unravel_index(np.unique(pindex,return_index=True,return_inverse=False)[1],pindex.shape,order='C'))
# kindex, rho = self.get_power_index(irreducible=True) return pindex,pundex,kindex,rho ## FIXME: order
#
# k = kindex[pindex]
# logmax = np.log(kindex.max()/kindex[1] + 1.)
# # lmax = np.ceil((logmax - np.log(2.))/lnfactor) + 1
# binedge = 0.
# ls_binned = np.zeros(self.dim(split=True),dtype=int)
# ls_binned[np.where(pindex==0)[0]] = 0
# kindex_binned = [0.]
# rho_binned = [(pindex==0).sum()]#*lm.get_meta_volume()/np.prod(lm.vol)
# i = 0
# while binedge < logmax:
# i += 1
# binedgeold = binedge
# binedge = lnfactor + binedgeold
# indices = np.where((np.log(kindex/kindex[1] + 1.) > binedgeold) & (np.log(kindex/kindex[1] + 1.) <= binedge))[0]
# if len(indices) == 0:
# binedge = np.log(np.exp(binedgeold) + np.min(self.vol))
# indices = np.where((np.log(kindex/kindex[1] + 1.) > binedgeold) & (np.log(kindex/kindex[1] + 1.) <= binedge))[0]
# if len(indices) == 0:
# raise NameError('empty bin')
# kindex_binned.append(np.exp(np.log(kindex[indices]).mean()))
# indices = np.array([np.where((np.log(k/kindex[1] + 1.) > binedgeold) & (np.log(k/kindex[1] + 1.) <= binedge))[k] for k in range(2)]).transpose()
# # indices = np.where((np.log(k/kindex[1] + 1.) > binedgeold) & (np.log(k/kindex[1] + 1.) <= binedge))[k] for k in range(2)]).transpose()
# print indices
# rho_binned.append(len(indices))
# ls_binned[indices[:,0],indices[:,1]] = i
# print ls_binned
# raw_input()
#
# return ls_binned, kindex, kindex_binned, rho_binned, lnfactor
# factor = 1
# lnfactor = np.log(factor)
# ls = lm.get_power_index()
# ls_irr, rho = lm.get_power_index(irreducible=True)
# ls_irr_long = ls_irr[ls]
# logmax = np.log(ls_irr.max()/ls_irr[1] + 1.)
# # lmax = np.ceil((logmax - np.log(2.))/lnfactor) + 1
# binedge = 0.
# ls_binned = np.zeros(lm.dim(split=True),dtype=int)
# ls_binned[np.where(ls==0)[0]] = 0
# ls_irr_binned = [0.]
# rho_binned = [(ls==0).sum()]#*lm.get_meta_volume()/np.prod(lm.vol)
# i = 0
# while binedge < logmax:
# i += 1
# binedgeold = binedge
# binedge = lnfactor + binedgeold
# indices = np.where((np.log(ls_irr/ls_irr[1] + 1.) > binedgeold) & (np.log(ls_irr/ls_irr[1] + 1.) <= binedge))[0]
# if len(indices) == 0:
# binedge = np.log(np.exp(binedgeold) + np.min(lm.vol))
# indices = np.where((np.log(ls_irr/ls_irr[1] + 1.) > binedgeold) & (np.log(ls_irr/ls_irr[1] + 1.) <= binedge))[0]
# if len(indices) == 0:
# raise NameError('empty bin')
# ls_irr_binned.append(np.exp(np.log(ls_irr[indices]).mean()))
# indices = np.array([np.where((np.log(ls_irr_long/ls_irr[1] + 1.) > binedgeold) & (np.log(ls_irr_long/ls_irr[1] + 1.) <= binedge))[k] for k in range(2)]).transpose()
# # indices = np.where((np.log(ls_irr_long/ls_irr[1] + 1.) > binedgeold) & (np.log(ls_irr_long/ls_irr[1] + 1.) <= binedge))[k] for k in range(2)]).transpose()
# print indices
# rho_binned.append(len(indices))
# ls_binned[indices[:,0],indices[:,1]] = i
# print ls_binned
# raw_input()
# return ls_binned, ls_irr, ls_irr_binned, rho_binned, lnfactor
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -8499,7 +8436,7 @@ class diagonal_operator(operator): ...@@ -8499,7 +8436,7 @@ class diagonal_operator(operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_random_field(self,domain=None,target=None,**kwargs): ## TODO: remove **kwargs for downward compatibility in future version def get_random_field(self,domain=None,target=None,**kwargs):
""" """
Generates a Gaussian random field with variance equal to the Generates a Gaussian random field with variance equal to the
diagonal. diagonal.
...@@ -8519,6 +8456,8 @@ class diagonal_operator(operator): ...@@ -8519,6 +8456,8 @@ class diagonal_operator(operator):
Random field. Random field.
""" """
if(len(kwargs)): ## TODO: remove **kwargs in future version
about.warnings.cprint("WARNING: deprecated keyword(s).")
## weight if ... ## weight if ...
if(not self.domain.discrete): if(not self.domain.discrete):
diag = self.domain.calc_weight(self.val,power=-1) diag = self.domain.calc_weight(self.val,power=-1)
...@@ -8776,7 +8715,7 @@ class power_operator(diagonal_operator): ...@@ -8776,7 +8715,7 @@ class power_operator(diagonal_operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_power(self,bare=True,pundex=None,pindex=None,**kwargs): ## TODO: remove **kwargs for downward compatibility in future version def get_power(self,bare=True,pundex=None,pindex=None,**kwargs):
""" """
Computes the power spectrum Computes the power spectrum
...@@ -8798,6 +8737,8 @@ class power_operator(diagonal_operator): ...@@ -8798,6 +8737,8 @@ class power_operator(diagonal_operator):
spec : ndarray spec : ndarray
The power spectrum The power spectrum
""" """
if(len(kwargs)): ## TODO: remove **kwargs in future version
about.warnings.cprint("WARNING: deprecated keyword(s).")
## weight if ... ## weight if ...
if(not self.domain.discrete)and(bare): if(not self.domain.discrete)and(bare):
diag = np.real(self.domain.calc_weight(self.val,power=-1)) diag = np.real(self.domain.calc_weight(self.val,power=-1))
...@@ -8810,45 +8751,6 @@ class power_operator(diagonal_operator): ...@@ -8810,45 +8751,6 @@ class power_operator(diagonal_operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# def get_random_field(self,domain=None,target=None,**kwargs): ## TODO: remove in future version
# """
# Generates a Gaussian random field with variance equal to the power spectrum
#
# Parameters
# ----------
# domain : space
# The space wherein the field lives (default: None,
# indicates to use self.domain)
# target : space
# The space wherein the transform of the field lives
# (default: None, indicates to use self.domain.target)
# pindex : numpy.ndarray, *optional*
# Indexing array giving the power spectrum index of each band
# (default: None).
# kindex : numpy.ndarray, *optional*
# Scale of each irreducible band (default: None).
#
# Returns
# -------
# x : field
# The random field defined on domain
# """
# if(np.any(self.val==0)):
# return super(power_operator,self).get_random_field(domain=domain,target=target,**kwargs)
#
# if(domain is None)or(domain==self.domain):
# if(np.any(self.val==0)):
# return super(power_operator,self).get_random_field(domain=self.domain,target=target,**kwargs)
# else:
# return field(self.domain,val=None,target=target,random="syn",spec=self.get_power(),**kwargs)
# else:
# if(np.any(self.val==0)):
# return super(power_operator,self).get_random_field(domain=self.domain,target=domain,**kwargs).transform(target=domain,overwrite=False)
# else:
# return field(self.domain,val=None,target=domain,random="syn",spec=self.get_power(),**kwargs).transform(target=domain,overwrite=False)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_projection_operator(self,pindex=None): def get_projection_operator(self,pindex=None):
""" """
Generates a spectral projection operator Generates a spectral projection operator
......
...@@ -20,9 +20,10 @@ ...@@ -20,9 +20,10 @@
## along with this program. If not, see <http://www.gnu.org/licenses/>. ## along with this program. If not, see <http://www.gnu.org/licenses/>.
## TODO: optimize ## TODO: optimize
## TODO: doc strings
from __future__ import division from __future__ import division
import gfft as gf #import gfft as gf
import numpy as np import numpy as np
...@@ -207,11 +208,16 @@ def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \ ...@@ -207,11 +208,16 @@ def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \
tfield = val tfield = val
vol = 1/np.array(val.shape)/vol vol = 1/np.array(val.shape)/vol
else: else:
#
tfield = gf.gfft(val, ftmachine='fft', \ if(zero_center):
in_zero_center=zero_center, out_zero_center=True, \ tfield = np.fft.fftshift(np.fft.fftn(np.fft.fftshift(val)))
enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \ else:
alpha=-1, verbose=False) tfield = np.fft.fftshift(np.fft.fftn(val))
# # Transform back to the signal space using GFFT.
# tfield = gf.gfft(val, ftmachine='fft', \
# in_zero_center=zero_center, out_zero_center=True, \
# enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \
# alpha=-1, verbose=False)
# Construct the Fourier transformed smoothing kernel # Construct the Fourier transformed smoothing kernel
tkernel = gaussian_kernel(val.shape, vol, smooth_length) tkernel = gaussian_kernel(val.shape, vol, smooth_length)
...@@ -221,12 +227,16 @@ def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \ ...@@ -221,12 +227,16 @@ def smooth_field(val, fourier, zero_center, enforce_hermitian_symmetry, vol, \
if(fourier): if(fourier):
sfield = tfield sfield = tfield
else: else:
#
# Transform back to the signal space using GFFT. if(zero_center):
sfield = gf.gfft(tfield, ftmachine='ifft', \ sfield = np.fft.ifftshift(np.fft.ifftn(np.fft.ifftshift(tfield)))
in_zero_center=True, out_zero_center=zero_center, \ else:
enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \ sfield = np.fft.ifftn(np.fft.ifftshift(tfield))
alpha=-1, verbose=False) # # Transform back to the signal space using GFFT.
# sfield = gf.gfft(tfield, ftmachine='ifft', \
# in_zero_center=True, out_zero_center=zero_center, \
# enforce_hermitian_symmetry=enforce_hermitian_symmetry, W=-1, \
# alpha=-1, verbose=False)
return sfield return sfield
......
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