Commit b47f67fd authored by M Selig's avatar M Selig

Merge pull request #11 from mselig/develop

verison update.
parents ca8d66e9 a5f253a1
...@@ -13,6 +13,7 @@ ...@@ -13,6 +13,7 @@
build build
demos/* demos/*
!demos/__init__.py
!demos/demo_faraday.py !demos/demo_faraday.py
!demos/demo_faraday_map.npy !demos/demo_faraday_map.npy
!demos/demo_excaliwir.py !demos/demo_excaliwir.py
......
...@@ -99,7 +99,7 @@ Requirements ...@@ -99,7 +99,7 @@ Requirements
Download Download
........ ........
The latest release is tagged **v0.8.0** and is available as a source package The latest release is tagged **v0.9.0** and is available as a source package
at `<https://github.com/mselig/nifty/tags>`_. The current version can be at `<https://github.com/mselig/nifty/tags>`_. The current version can be
obtained by cloning the repository:: obtained by cloning the repository::
......
...@@ -26,7 +26,7 @@ from nifty_power import * ...@@ -26,7 +26,7 @@ from nifty_power import *
from nifty_tools import * from nifty_tools import *
from nifty_explicit import * from nifty_explicit import *
from nifty_power_conversion import * #from nifty_power_conversion import *
##----------------------------------------------------------------------------- ##-----------------------------------------------------------------------------
......
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2014 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
## See the GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
pass
...@@ -237,4 +237,43 @@ class ncmap(object): ...@@ -237,4 +237,43 @@ class ncmap(object):
return cm("Plus Minus", segmentdata, N=int(ncolors), gamma=1.0) return cm("Plus Minus", segmentdata, N=int(ncolors), gamma=1.0)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@staticmethod
def planck(ncolors=256):
"""
Returns a color map similar to the one used for the "Planck CMB Map".
Parameters
----------
ncolors : int, *optional*
Number of color segments (default: 256).
Returns
-------
cmap : matplotlib.colors.LinearSegmentedColormap instance
Linear segmented color map.
"""
segmentdata = {"red": [(0.0, 0.00, 0.00), (0.1, 0.00, 0.00),
(0.2, 0.00, 0.00), (0.3, 0.00, 0.00),
(0.4, 0.00, 0.00), (0.5, 1.00, 1.00),
(0.6, 1.00, 1.00), (0.7, 1.00, 1.00),
(0.8, 0.83, 0.83), (0.9, 0.67, 0.67),
(1.0, 0.50, 0.50)],
"green": [(0.0, 0.00, 0.00), (0.1, 0.00, 0.00),
(0.2, 0.00, 0.00), (0.3, 0.30, 0.30),
(0.4, 0.70, 0.70), (0.5, 1.00, 1.00),
(0.6, 0.70, 0.70), (0.7, 0.30, 0.30),
(0.8, 0.00, 0.00), (0.9, 0.00, 0.00),
(1.0, 0.00, 0.00)],
"blue": [(0.0, 0.50, 0.50), (0.1, 0.67, 0.67),
(0.2, 0.83, 0.83), (0.3, 1.00, 1.00),
(0.4, 1.00, 1.00), (0.5, 1.00, 1.00),
(0.6, 0.00, 0.00), (0.7, 0.00, 0.00),
(0.8, 0.00, 0.00), (0.9, 0.00, 0.00),
(1.0, 0.00, 0.00)]}
return cm("Planck-like", segmentdata, N=int(ncolors), gamma=1.0)
##----------------------------------------------------------------------------- ##-----------------------------------------------------------------------------
...@@ -514,7 +514,7 @@ class _about(object): ## nifty support class for global settings ...@@ -514,7 +514,7 @@ class _about(object): ## nifty support class for global settings
""" """
## version ## version
self._version = "0.8.4" self._version = "0.9.0"
## switches and notifications ## switches and notifications
self._errors = notification(default=True,ccode=notification._code) self._errors = notification(default=True,ccode=notification._code)
...@@ -2485,7 +2485,7 @@ class rg_space(space): ...@@ -2485,7 +2485,7 @@ class rg_space(space):
config = {"binbounds":kwargs.get("binbounds",None),"log":kwargs.get("log",None),"nbin":kwargs.get("nbin",None)} config = {"binbounds":kwargs.get("binbounds",None),"log":kwargs.get("log",None),"nbin":kwargs.get("nbin",None)}
## power indices ## power indices
about.infos.cflush("INFO: setting power indices ...") about.infos.cflush("INFO: setting power indices ...")
pindex,kindex,rho = gp.get_power_indices(self.para[:(np.size(self.para)-1)//2],self.vol,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=True) pindex,kindex,rho = gp.get_power_indices2(self.para[:(np.size(self.para)-1)//2],self.vol,self.para[-((np.size(self.para)-1)//2):].astype(np.bool),fourier=True)
## bin if ... ## bin if ...
if(config.get("log") is not None)or(config.get("nbin") is not None)or(config.get("binbounds") is not None): if(config.get("log") is not None)or(config.get("nbin") is not None)or(config.get("binbounds") is not None):
pindex,kindex,rho = gp.bin_power_indices(pindex,kindex,rho,**config) pindex,kindex,rho = gp.bin_power_indices(pindex,kindex,rho,**config)
...@@ -3195,6 +3195,9 @@ class rg_space(space): ...@@ -3195,6 +3195,9 @@ class rg_space(space):
about.infos.cprint("INFO: absolute values and phases are plotted.") about.infos.cprint("INFO: absolute values and phases are plotted.")
if(title): if(title):
title += " " title += " "
if(bool(kwargs.get("save",False))):
save_ = os.path.splitext(os.path.basename(str(kwargs.get("save"))))
kwargs.update(save=save_[0]+"_absolute"+save_[1])
self.get_plot(np.absolute(x),title=title+"(absolute)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) self.get_plot(np.absolute(x),title=title+"(absolute)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) # self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) # self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
...@@ -3202,6 +3205,8 @@ class rg_space(space): ...@@ -3202,6 +3205,8 @@ class rg_space(space):
unit = "rad" unit = "rad"
if(cmap is None): if(cmap is None):
cmap = pl.cm.hsv_r cmap = pl.cm.hsv_r
if(bool(kwargs.get("save",False))):
kwargs.update(save=save_[0]+"_phase"+save_[1])
self.get_plot(np.angle(x,deg=False),title=title+"(phase)",vmin=-3.1416,vmax=3.1416,power=False,unit=unit,norm=None,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) ## values in [-pi,pi] self.get_plot(np.angle(x,deg=False),title=title+"(phase)",vmin=-3.1416,vmax=3.1416,power=False,unit=unit,norm=None,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) ## values in [-pi,pi]
return None ## leave method return None ## leave method
else: else:
...@@ -4011,11 +4016,16 @@ class lm_space(space): ...@@ -4011,11 +4016,16 @@ class lm_space(space):
if(np.iscomplexobj(x)): if(np.iscomplexobj(x)):
if(title): if(title):
title += " " title += " "
if(bool(kwargs.get("save",False))):
save_ = os.path.splitext(os.path.basename(str(kwargs.get("save"))))
kwargs.update(save=save_[0]+"_absolute"+save_[1])
self.get_plot(np.absolute(x),title=title+"(absolute)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) self.get_plot(np.absolute(x),title=title+"(absolute)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) # self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) # self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
if(cmap is None): if(cmap is None):
cmap = pl.cm.hsv_r cmap = pl.cm.hsv_r
if(bool(kwargs.get("save",False))):
kwargs.update(save=save_[0]+"_phase"+save_[1])
self.get_plot(np.angle(x,deg=False),title=title+"(phase)",vmin=-3.1416,vmax=3.1416,power=False,norm=None,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) ## values in [-pi,pi] self.get_plot(np.angle(x,deg=False),title=title+"(phase)",vmin=-3.1416,vmax=3.1416,power=False,norm=None,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) ## values in [-pi,pi]
return None ## leave method return None ## leave method
else: else:
...@@ -5413,7 +5423,7 @@ class nested_space(space): ...@@ -5413,7 +5423,7 @@ class nested_space(space):
self.datatype = self.nest[-1].datatype self.datatype = self.nest[-1].datatype
self.discrete = np.prod([nn.discrete for nn in self.nest],axis=0,dtype=np.bool,out=None) self.discrete = np.prod([nn.discrete for nn in self.nest],axis=0,dtype=np.bool,out=None)
self.vol = None ## not needed self.vol = np.prod([nn.get_meta_volume(total=True) for nn in self.nest],axis=0,dtype=None,out=None) ## total volume
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -5738,7 +5748,7 @@ class nested_space(space): ...@@ -5738,7 +5748,7 @@ class nested_space(space):
""" """
if(total): if(total):
## product ## product
return np.prod([nn.get_meta_volume(total=True) for nn in self.nest],axis=0,dtype=None,out=None) return self.vol ## == np.prod([nn.get_meta_volume(total=True) for nn in self.nest],axis=0,dtype=None,out=None)
else: else:
mol = self.nest[0].get_meta_volume(total=False) mol = self.nest[0].get_meta_volume(total=False)
## tensor product ## tensor product
...@@ -8280,26 +8290,52 @@ class operator(object): ...@@ -8280,26 +8290,52 @@ class operator(object):
def det(self): def det(self):
""" """
Computes the determinant of the operator Computes the determinant of the operator.
Returns Returns
------- -------
det : float det : float
The determinant The determinant
""" """
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'det'.")) raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'det'."))
def inverse_det(self): def inverse_det(self):
""" """
Computes the determinant of the inverse operator Computes the determinant of the inverse operator.
Returns Returns
------- -------
det : float det : float
The determinant The determinant
""" """
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'inverse_det'.")) raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'inverse_det'."))
def log_det(self):
"""
Computes the logarithm of the determinant of the operator (if applicable).
Returns
-------
logdet : float
The logarithm of the determinant
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'log_det'."))
def tr_log(self):
"""
Computes the trace of the logarithm of the operator (if applicable).
Returns
-------
logdet : float
The trace of the logarithm
"""
return self.log_det()
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def hat(self,bare=False,domain=None,target=None,**kwargs): def hat(self,bare=False,domain=None,target=None,**kwargs):
...@@ -9119,6 +9155,23 @@ class diagonal_operator(operator): ...@@ -9119,6 +9155,23 @@ class diagonal_operator(operator):
else: else:
raise ValueError(about._errors.cstring("ERROR: singular operator.")) raise ValueError(about._errors.cstring("ERROR: singular operator."))
def log_det(self):
"""
Computes the logarithm of the determinant of the operator.
Returns
-------
logdet : float
The logarithm of the determinant
"""
if(self.uni): ## identity
return 0
elif(self.domain.dim(split=False)<self.domain.dof()): ## hidden degrees of freedom
return self.domain.calc_dot(np.ones(self.domain.dim(split=True),dtype=self.domain.datatype,order='C'),np.log(self.val))
else:
return np.sum(np.log(self.val),axis=None,dtype=None,out=None)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_random_field(self,domain=None,target=None,**kwargs): def get_random_field(self,domain=None,target=None,**kwargs):
...@@ -10201,7 +10254,8 @@ class vecvec_operator(operator): ...@@ -10201,7 +10254,8 @@ class vecvec_operator(operator):
def inverse_diag(self): def inverse_diag(self):
""" """
Inverse is ill-defined for this operator. Inverse is ill-defined for this operator.
""" """
raise AttributeError(about._errors.cstring("ERROR: singular operator.")) raise AttributeError(about._errors.cstring("ERROR: singular operator."))
...@@ -10209,18 +10263,27 @@ class vecvec_operator(operator): ...@@ -10209,18 +10263,27 @@ class vecvec_operator(operator):
def det(self): def det(self):
""" """
Computes the determinant of the operator Computes the determinant of the operator.
Returns Returns
------- -------
det : 0 det : 0
The determinant The determinant
""" """
return 0 return 0
def inverse_det(self): def inverse_det(self):
""" """
Inverse is ill-defined for this operator. Inverse is ill-defined for this operator.
"""
raise AttributeError(about._errors.cstring("ERROR: singular operator."))
def log_det(self):
"""
Logarithm of the determinant is ill-defined for this singular operator.
""" """
raise AttributeError(about._errors.cstring("ERROR: singular operator.")) raise AttributeError(about._errors.cstring("ERROR: singular operator."))
......
...@@ -1077,7 +1077,11 @@ class explicit_operator(operator): ...@@ -1077,7 +1077,11 @@ class explicit_operator(operator):
if(np.any(self._hidden)): if(np.any(self._hidden)):
about.warnings.cprint("WARNING: inappropriate determinant calculation.") about.warnings.cprint("WARNING: inappropriate determinant calculation.")
return np.linalg.det(self.weight(rowpower=0.5,colpower=0.5,overwrite=False)) det = np.linalg.det(self.weight(rowpower=0.5,colpower=0.5,overwrite=False))
if(np.isreal(det)):
return np.asscalar(np.real(det))
else:
return det
def inverse_det(self): def inverse_det(self):
""" """
...@@ -1104,6 +1108,44 @@ class explicit_operator(operator): ...@@ -1104,6 +1108,44 @@ class explicit_operator(operator):
else: else:
raise ValueError(about._errors.cstring("ERROR: singular matrix.")) raise ValueError(about._errors.cstring("ERROR: singular matrix."))
def log_det(self):
"""
Computes the logarithm of the determinant of the operator
(if applicable).
Returns
-------
logdet : float
The logarithm of the determinant
See Also
--------
numpy.linalg.slogdet
Raises
------
ValueError
If `domain` and `target` are unequal or it is non-positive
definite matrix.
"""
if(self.domain!=self.target):
raise ValueError(about._errors.cstring("ERROR: determinant ill-defined."))
if(np.any(self._hidden)):
about.warnings.cprint("WARNING: inappropriate determinant calculation.")
sign,logdet = np.linalg.slogdet(self.weight(rowpower=0.5,colpower=0.5,overwrite=False))
if(abs(sign)<0.1): ## abs(sign) << 1
raise ValueError(about._errors.cstring("ERROR: singular matrix."))
if(sign==-1):
raise ValueError(about._errors.cstring("ERROR: non-positive definite matrix."))
else:
logdet += np.log(sign)
if(np.isreal(logdet)):
return np.asscalar(np.real(logdet))
else:
return logdet
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def __len__(self): def __len__(self):
...@@ -2324,7 +2366,11 @@ def explicify(op,newdomain=None,newtarget=None,ncpu=2,nper=1,loop=False,**quargs ...@@ -2324,7 +2366,11 @@ def explicify(op,newdomain=None,newtarget=None,ncpu=2,nper=1,loop=False,**quargs
raise TypeError(about._errors.cstring("ERROR: invalid input.")) raise TypeError(about._errors.cstring("ERROR: invalid input."))
elif(newdomain is not None)and(newtarget is None)and(op.domain==op.target): elif(newdomain is not None)and(newtarget is None)and(op.domain==op.target):
newtarget = newdomain newtarget = newdomain
return explicit_probing(op=op,function=op.times,domain=newdomain,codomain=newtarget,target=op.domain,ncpu=ncpu,nper=nper,**quargs)(loop=loop) if(newdomain is None)or(newdomain==op.domain):
target_ = None
else:
target_ = op.domain
return explicit_probing(op=op,function=op.times,domain=newdomain,codomain=newtarget,target=target_,ncpu=ncpu,nper=nper,**quargs)(loop=loop)
##----------------------------------------------------------------------------- ##-----------------------------------------------------------------------------
...@@ -575,7 +575,7 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None ...@@ -575,7 +575,7 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
##----------------------------------------------------------------------------- ##-----------------------------------------------------------------------------
def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None,**kwargs): def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None,loglog=False,**kwargs):
""" """
Interpolates a given power spectrum at new k(-indices). Interpolates a given power spectrum at new k(-indices).
...@@ -604,6 +604,9 @@ def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None, ...@@ -604,6 +604,9 @@ def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None,
Scales corresponding to each band in the new power spectrum; Scales corresponding to each band in the new power spectrum;
can be retrieved from `domain` if `kindex` is given can be retrieved from `domain` if `kindex` is given
(default: None). (default: None).
loglog : bool, *optional*
Flag specifying whether the interpolation is done on log-log-scale
(ignoring the monopole) or not (default: False)
Returns Returns
------- -------
...@@ -719,7 +722,13 @@ def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None, ...@@ -719,7 +722,13 @@ def interpolate_power(spec,mode="linear",domain=None,kindex=None,newkindex=None,
spec = np.r_[spec,np.exp(2*np.log(spec[-1])-np.log(spec[-nmirror:-1][::-1]))] spec = np.r_[spec,np.exp(2*np.log(spec[-1])-np.log(spec[-nmirror:-1][::-1]))]
kindex = np.r_[kindex,(2*kindex[-1]-kindex[-nmirror:-1][::-1])] kindex = np.r_[kindex,(2*kindex[-1]-kindex[-nmirror:-1][::-1])]
## interpolation ## interpolation
newspec = ip(kindex,spec,kind=mode,axis=0,copy=True,bounds_error=True,fill_value=np.NAN)(newkindex) if(loglog):
## map to log-log-scale ignoring monopole
logspec = np.log(spec[1:])
logspec = ip(np.log(kindex[1:]),logspec,kind=mode,axis=0,copy=True,bounds_error=True,fill_value=np.NAN)(np.log(newkindex[1:]))
newspec = np.concatenate([[spec[0]],np.exp(logspec)])
else:
newspec = ip(kindex,spec,kind=mode,axis=0,copy=True,bounds_error=True,fill_value=np.NAN)(newkindex)
## check new power spectrum ## check new power spectrum
if(not np.all(np.isfinite(newspec))): if(not np.all(np.isfinite(newspec))):
raise ValueError(about._errors.cstring("ERROR: infinite value(s).")) raise ValueError(about._errors.cstring("ERROR: infinite value(s)."))
......
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2014 Max-Planck-Society
##
## Author: Maksim Greiner
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
## See the GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
from nifty import * from nifty import *
def power_backward_conversion_rg(k_space,p,mean=None,bare=True): def power_backward_conversion_rg(k_space,p,mean=None,bare=True):
...@@ -16,12 +37,12 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True): ...@@ -16,12 +37,12 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True):
Needs to have the same number of entries as Needs to have the same number of entries as
`k_space.get_power_indices()[0]` `k_space.get_power_indices()[0]`
mean : float, *optional* mean : float, *optional*
specifies the mean of the log-normal field. If `mean` is not specifies the mean of the log-normal field. If `mean` is not
specified the function will use the monopole of the power spectrum. specified the function will use the monopole of the power spectrum.
If it is specified the function will NOT use the monopole of the If it is specified the function will NOT use the monopole of the
spectrum (default: None). spectrum (default: None).
WARNING: a mean that is too low can violate positive definiteness WARNING: a mean that is too low can violate positive definiteness
of the log-normal field. In this case the function produces an of the log-normal field. In this case the function produces an
error. error.
bare : bool, *optional* bare : bool, *optional*
whether `p` is the bare power spectrum or not (default: True). whether `p` is the bare power spectrum or not (default: True).
...@@ -34,7 +55,7 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True): ...@@ -34,7 +55,7 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True):
the power spectrum of the underlying Gaussian field, where the the power spectrum of the underlying Gaussian field, where the
monopole has been set to zero. Eventual monopole power has been monopole has been set to zero. Eventual monopole power has been
shifted to the mean. shifted to the mean.
References References
---------- ----------
.. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum"; .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
...@@ -42,41 +63,41 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True): ...@@ -42,41 +63,41 @@ def power_backward_conversion_rg(k_space,p,mean=None,bare=True):
""" """
pindex = k_space.get_power_indices()[2] pindex = k_space.get_power_indices()[2]
V = k_space.vol.prod()**(-1) V = k_space.vol.prod()**(-1)
mono_ind = np.where(pindex==0) mono_ind = np.where(pindex==0)
spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False) spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
if(mean is None): if(mean is None):
mean = 0. mean = 0.
else: else:
spec[0] = 0. spec[0] = 0.
pf = field(k_space,val=spec[pindex]).transform()+mean**2 pf = field(k_space,val=spec[pindex]).transform()+mean**2
if(np.any(pf.val<0.)): if(np.any(pf.val<0.)):
raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean.")) raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
return None return None
p1 = sqrt(log(pf).power()) p1 = sqrt(log(pf).power())
p1[0] = (log(pf)).transform()[mono_ind][0] p1[0] = (log(pf)).transform()[mono_ind][0]
p2 = 0.5*V*log(k_space.calc_weight(spec[pindex],1).sum()+mean**2) p2 = 0.5*V*log(k_space.calc_weight(spec[pindex],1).sum()+mean**2)
logmean = 1/V * (p1[0]-p2) logmean = 1/V * (p1[0]-p2)
p1[0] = 0. p1[0] = 0.
if(np.any(p1<0.)): if(np.any(p1<0.)):
raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean.")) raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
return None return None
if(bare==True): if(bare==True):
return logmean.real,power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real return logmean.real,power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
else: else:
return logmean.real,p1.real return logmean.real,p1.real
def power_forward_conversion_rg(k_space,p,mean=0,bare=True): def power_forward_conversion_rg(k_space,p,mean=0,bare=True):
""" """
This function is designed to convert a theoretical/statistical power This function is designed to convert a theoretical/statistical power
...@@ -101,31 +122,31 @@ def power_forward_conversion_rg(k_space,p,mean=0,bare=True): ...@@ -101,31 +122,31 @@ def power_forward_conversion_rg(k_space,p,mean=0,bare=True):
------- -------
p1 : np.array, p1 : np.array,
the power spectrum of the exponentiated Gaussian field. the power spectrum of the exponentiated Gaussian field.
References References
---------- ----------
.. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum"; .. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
`arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_ `arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
""" """
pindex = k_space.get_power_indices()[2] pindex = k_space.get_power_indices()[2]
spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False) spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
S_x = field(k_space,val=spec[pindex]).transform() S_x = field(k_space,val=spec[pindex]).transform()
S_0 = k_space.calc_weight(spec[pindex],1).sum() S_0 = k_space.calc_weight(spec[pindex],1).sum()
pf = exp(S_x+S_0+2*mean) pf = exp(S_x+S_0+2*mean)
p1 = sqrt(pf.power()) p1 = sqrt(pf.power())
if(bare==True): if(bare==True):
return power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real return power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
else: else:
return p1.real return p1.real
def power_backward_conversion_lm(k_space,p,mean=None): def power_backward_conversion_lm(k_space,p,mean=None):
""" """
...@@ -143,12 +164,12 @@ def power_backward_conversion_lm(k_space,p,mean=None): ...@@ -143,12 +164,12 @@ def power_backward_conversion_lm(k_space,p,mean=None):
Needs to have the same number of entries as Needs to have the same number of entries as
`k_space.get_power_indices()[0]` `k_space.get_power_indices()[0]`
mean : float, *optional* mean : float, *optional*
specifies the mean of the log-normal field. If `mean` is not specifies the mean of the log-normal field. If `mean` is not
specified the function will use the monopole of the power spectrum. specified the function will use the monopole of the power spectrum.
If it is specified the function will NOT use the monopole of the If it is specified the function will NOT use the monopole of the
spectrum. (default: None) spectrum. (default: None)
WARNING: a mean that is too low can violate positive definiteness WARNING: a mean that is too low can violate positive definiteness
of the log-normal field. In this case the function produces an of the log-normal field. In this case the function produces an
error. error.