diff --git a/.gitignore b/.gitignore
index 6991ad8a9491a97c929d18d68cbda3343a339b73..ebdda7ed5cb0ca821b845769630e81d438ba9731 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,6 +13,7 @@
build
demos/*
+!demos/__init__.py
!demos/demo_faraday.py
!demos/demo_faraday_map.npy
!demos/demo_excaliwir.py
diff --git a/README.rst b/README.rst
index 9e9f12d275c2f67a435e49e9a80316b9bfe9523c..005dd18c5c6b2def265046ab476b7bd700d5bb4b 100644
--- a/README.rst
+++ b/README.rst
@@ -99,7 +99,7 @@ Requirements
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
obtained by cloning the repository::
diff --git a/__init__.py b/__init__.py
index 564ff2ecefcd9fa0d6030dff5673895eee08bc74..4121b411ac2e59bfc35bd7008e16beea14fcb4b6 100644
--- a/__init__.py
+++ b/__init__.py
@@ -26,7 +26,7 @@ from nifty_power import *
from nifty_tools import *
from nifty_explicit import *
-from nifty_power_conversion import *
+#from nifty_power_conversion import *
##-----------------------------------------------------------------------------
diff --git a/demos/__init__.py b/demos/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e9456ac3f7be17e93bbf5678e3670e4cddb9137d
--- /dev/null
+++ b/demos/__init__.py
@@ -0,0 +1,23 @@
+## 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
+
diff --git a/nifty_cmaps.py b/nifty_cmaps.py
index 34fce2a01a9534a6a0ed37a9f956220e48be7aae..014367951b74e051b11024bf6a1acbc3dbba51ae 100644
--- a/nifty_cmaps.py
+++ b/nifty_cmaps.py
@@ -237,4 +237,43 @@ class ncmap(object):
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)
+
##-----------------------------------------------------------------------------
diff --git a/nifty_core.py b/nifty_core.py
index 338f35c543ff2315f201c91283eeec3e5604694f..4a72131b8704d31455199fcd7dfd38430ed8d59c 100644
--- a/nifty_core.py
+++ b/nifty_core.py
@@ -514,7 +514,7 @@ class _about(object): ## nifty support class for global settings
"""
## version
- self._version = "0.8.4"
+ self._version = "0.9.0"
## switches and notifications
self._errors = notification(default=True,ccode=notification._code)
@@ -2485,7 +2485,7 @@ class rg_space(space):
config = {"binbounds":kwargs.get("binbounds",None),"log":kwargs.get("log",None),"nbin":kwargs.get("nbin",None)}
## 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 ...
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)
@@ -3195,6 +3195,9 @@ class rg_space(space):
about.infos.cprint("INFO: absolute values and phases are plotted.")
if(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.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)
@@ -3202,6 +3205,8 @@ class rg_space(space):
unit = "rad"
if(cmap is None):
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]
return None ## leave method
else:
@@ -4011,11 +4016,16 @@ class lm_space(space):
if(np.iscomplexobj(x)):
if(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.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)
if(cmap is None):
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]
return None ## leave method
else:
@@ -5413,7 +5423,7 @@ class nested_space(space):
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.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):
"""
if(total):
## 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:
mol = self.nest[0].get_meta_volume(total=False)
## tensor product
@@ -8280,26 +8290,52 @@ class operator(object):
def det(self):
"""
- Computes the determinant of the operator
+ Computes the determinant of the operator.
Returns
-------
det : float
The determinant
+
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'det'."))
def inverse_det(self):
"""
- Computes the determinant of the inverse operator
+ Computes the determinant of the inverse operator.
Returns
-------
det : float
The determinant
+
"""
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):
@@ -9119,6 +9155,23 @@ class diagonal_operator(operator):
else:
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):
@@ -10201,7 +10254,8 @@ class vecvec_operator(operator):
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."))
@@ -10209,18 +10263,27 @@ class vecvec_operator(operator):
def det(self):
"""
- Computes the determinant of the operator
+ Computes the determinant of the operator.
Returns
-------
det : 0
The determinant
+
"""
return 0
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."))
diff --git a/nifty_explicit.py b/nifty_explicit.py
index 5e2a11072d4872bea877c1a0b278b22c1f9b0238..29c5dd5ba231deb45e45ec64a468b5b85056012a 100644
--- a/nifty_explicit.py
+++ b/nifty_explicit.py
@@ -1077,7 +1077,11 @@ class explicit_operator(operator):
if(np.any(self._hidden)):
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):
"""
@@ -1104,6 +1108,44 @@ class explicit_operator(operator):
else:
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):
@@ -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."))
elif(newdomain is not None)and(newtarget is None)and(op.domain==op.target):
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)
##-----------------------------------------------------------------------------
diff --git a/nifty_power.py b/nifty_power.py
index 00534780814d571703f88442d3c8393b32217453..9c46ab0e364013afef61f9ff6066bc82cf01bafe 100644
--- a/nifty_power.py
+++ b/nifty_power.py
@@ -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).
@@ -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;
can be retrieved from `domain` if `kindex` is given
(default: None).
+ loglog : bool, *optional*
+ Flag specifying whether the interpolation is done on log-log-scale
+ (ignoring the monopole) or not (default: False)
Returns
-------
@@ -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]))]
kindex = np.r_[kindex,(2*kindex[-1]-kindex[-nmirror:-1][::-1])]
## 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
if(not np.all(np.isfinite(newspec))):
raise ValueError(about._errors.cstring("ERROR: infinite value(s)."))
diff --git a/nifty_power_conversion.py b/nifty_power_conversion.py
index 999e18f32e8090f58fcfa0719fc9af21904e2681..50e820d67c1d87aa6dbd9a5fdf6a673f1ca450bf 100644
--- a/nifty_power_conversion.py
+++ b/nifty_power_conversion.py
@@ -1,3 +1,24 @@
+## 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 *
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
`k_space.get_power_indices()[0]`
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.
- 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).
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.
bare : bool, *optional*
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):
the power spectrum of the underlying Gaussian field, where the
monopole has been set to zero. Eventual monopole power has been
shifted to the mean.
-
+
References
----------
.. [#] 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):
"""
pindex = k_space.get_power_indices()[2]
V = k_space.vol.prod()**(-1)
-
+
mono_ind = np.where(pindex==0)
-
+
spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
-
+
if(mean is None):
mean = 0.
else:
spec[0] = 0.
-
+
pf = field(k_space,val=spec[pindex]).transform()+mean**2
-
+
if(np.any(pf.val<0.)):
raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
return None
-
+
p1 = sqrt(log(pf).power())
-
+
p1[0] = (log(pf)).transform()[mono_ind][0]
-
+
p2 = 0.5*V*log(k_space.calc_weight(spec[pindex],1).sum()+mean**2)
-
+
logmean = 1/V * (p1[0]-p2)
-
+
p1[0] = 0.
-
+
if(np.any(p1<0.)):
raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
return None
-
+
if(bare==True):
return logmean.real,power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
else:
return logmean.real,p1.real
-
+
def power_forward_conversion_rg(k_space,p,mean=0,bare=True):
"""
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):
-------
p1 : np.array,
the power spectrum of the exponentiated Gaussian field.
-
+
References
----------
.. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
`arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
"""
-
+
pindex = k_space.get_power_indices()[2]
-
+
spec = power_operator(k_space,spec=p,bare=bare).get_power(bare=False)
-
+
S_x = field(k_space,val=spec[pindex]).transform()
-
+
S_0 = k_space.calc_weight(spec[pindex],1).sum()
-
+
pf = exp(S_x+S_0+2*mean)
-
+
p1 = sqrt(pf.power())
-
+
if(bare==True):
return power_operator(k_space,spec=p1,bare=False).get_power(bare=True).real
else:
return p1.real
-
-
+
+
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
`k_space.get_power_indices()[0]`
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.
- 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)
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.
Returns
@@ -159,13 +180,13 @@ def power_backward_conversion_lm(k_space,p,mean=None):
the power spectrum of the underlying Gaussian field, where the
monopole has been set to zero. Eventual monopole power has been
shifted to the mean.
-
+
References
----------
.. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
`arXiv:1312.1354 <http://arxiv.org/abs/1312.1354>`_
"""
-
+
p = np.copy(p)
if(mean is not None):
p[0] = 4*pi*mean**2
@@ -174,7 +195,7 @@ def power_backward_conversion_lm(k_space,p,mean=None):
C_0_Omega = field(k_space,val=0)
C_0_Omega.val[:len(klen)] = p*sqrt(2*klen+1)/sqrt(4*pi)
C_0_Omega = C_0_Omega.transform()
-
+
if(np.any(C_0_Omega.val<0.)):
raise ValueError(about._errors.cstring("ERROR: spectrum or mean incompatible with positive definiteness.\n Try increasing the mean."))
return None
@@ -186,20 +207,20 @@ def power_backward_conversion_lm(k_space,p,mean=None):
spec = Z.val[:len(klen)]
mean = (spec[0]-0.5*sqrt(4*pi)*log((p*(2*klen+1)/(4*pi)).sum()))/sqrt(4*pi)
-
+
spec[0] = 0.
spec = spec*sqrt(4*pi)/sqrt(2*klen+1)
-
+
spec = np.real(spec)
-
+
if(np.any(spec<0.)):
spec = spec*(spec>0.)
about.warnings.cprint("WARNING: negative modes set to zero.")
return mean.real,spec
-
-
+
+
def power_forward_conversion_lm(k_space,p,mean=0):
"""
This function is designed to convert a theoretical/statistical power
@@ -217,12 +238,12 @@ def power_forward_conversion_lm(k_space,p,mean=0):
`k_space.get_power_indices()[0]`
m : float, *optional*
specifies the mean of the Gaussian field (default: 0).
-
+
Returns
-------
p1 : np.array,
the power spectrum of the exponentiated Gaussian field.
-
+
References
----------
.. [#] M. Greiner and T.A. Ensslin, "Log-transforming the matter power spectrum";
@@ -233,7 +254,7 @@ def power_forward_conversion_lm(k_space,p,mean=0):
C_0_Omega = field(k_space,val=0)
C_0_Omega.val[:len(klen)] = p*sqrt(2*klen+1)/sqrt(4*pi)
C_0_Omega = C_0_Omega.transform()
-
+
C_0_0 = (p*(2*klen+1)/(4*pi)).sum()
exC = exp(C_0_Omega+C_0_0+2*m)
@@ -243,9 +264,9 @@ def power_forward_conversion_lm(k_space,p,mean=0):
spec = Z.val[:len(klen)]
spec = spec*sqrt(4*pi)/sqrt(2*klen+1)
-
+
spec = np.real(spec)
-
+
if(np.any(spec<0.)):
spec = spec*(spec>0.)
about.warnings.cprint("WARNING: negative modes set to zero.")
diff --git a/nifty_tools.py b/nifty_tools.py
index a91d83881b4d81dffaaa04d0d56f960e3dbb6b23..582fdee054f3422aa4173d1e249f0d8d7527c510 100644
--- a/nifty_tools.py
+++ b/nifty_tools.py
@@ -774,10 +774,7 @@ class conjugate_gradient(object):
delta = delta_*np.absolute(gamma)**0.5
self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,np.real(alpha),np.real(beta),np.real(delta)))
- if(ii==limii):
- self.note.cprint("\n... quit.")
- break
- elif(gamma==0):
+ if(gamma==0):
convergence = clevel+1
self.note.cprint(" convergence level : INF\n... done.")
break
@@ -789,6 +786,9 @@ class conjugate_gradient(object):
break
else:
convergence = max(0,convergence-1)
+ if(ii==limii):
+ self.note.cprint("\n... quit.")
+ break
if(self.spam is not None):
self.spam(self.x,ii)
@@ -843,10 +843,7 @@ class conjugate_gradient(object):
delta = delta_*np.absolute(gamma)**0.5
self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,np.real(alpha),np.real(beta),np.real(delta)))
- if(ii==limii):
- self.note.cprint("\n... quit.")
- break
- elif(gamma==0):
+ if(gamma==0):
convergence = clevel+1
self.note.cprint(" convergence level : INF\n... done.")
break
@@ -858,6 +855,9 @@ class conjugate_gradient(object):
break
else:
convergence = max(0,convergence-1)
+ if(ii==limii):
+ self.note.cprint("\n... quit.")
+ break
if(self.spam is not None):
self.spam(self.x,ii)
@@ -1066,10 +1066,7 @@ class steepest_descent(object):
alpha *= a
norm = g.norm() ## gradient norm
- if(ii==limii):
- self.note.cprint("\n... quit.")
- break
- elif(delta==0):
+ if(delta==0):
convergence = clevel+2
self.note.cprint(" convergence level : %u\n... done."%convergence)
break
@@ -1082,6 +1079,9 @@ class steepest_descent(object):
break
else:
convergence = max(0,convergence-1)
+ if(ii==limii):
+ self.note.cprint("\n... quit.")
+ break
if(self.spam is not None):
self.spam(self.x,ii)
diff --git a/powerspectrum.py b/powerspectrum.py
index de135bb57975b8ee1407224f57a83e146769fdb2..b4938708e02461428cabaddf7de8aaf20f979a5f 100644
--- a/powerspectrum.py
+++ b/powerspectrum.py
@@ -385,6 +385,66 @@ def get_power_indices(axes,dgrid,zerocentered,fourier=True):
return ind,klength,rho
+def get_power_indices2(axes,dgrid,zerocentered,fourier=True):
+ """
+ Returns the index of the Fourier grid points in a numpy
+ array, ordered following the zerocentered flag.
+
+ Parameters
+ ----------
+ axes : ndarray
+ An array with the length of each axis.
+
+ dgrid : ndarray
+ An array with the pixel length of each axis.
+
+ zerocentered : bool
+ 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.
+
+ irred : bool : *optional*
+ If True, the function returns an array of all k-vector lengths and
+ their degeneracy factors. If False, just the power index array is
+ returned.
+
+ fourier : bool : *optional*
+ Whether the output should be in Fourier space or not
+ (default=False).
+
+ Returns
+ -------
+ index, klength, rho : ndarrays
+ Returns the power index array, an array of all k-vector lengths and
+ their degeneracy factors.
+
+ """
+
+ ## kdict, klength
+ if(np.any(zerocentered==False)):
+ kdict = np.fft.fftshift(nkdict_fast2(axes,dgrid,fourier),axes=shiftaxes(zerocentered,st_to_zero_mode=True))
+ else:
+ kdict = nkdict_fast2(axes,dgrid,fourier)
+
+ klength,rho,ind = nkdict_to_indices(kdict)
+
+ return ind,klength,rho
+
+
+def nkdict_to_indices(kdict):
+
+ kindex,pindex = np.unique(kdict,return_inverse=True)
+ pindex = pindex.reshape(kdict.shape)
+
+ rho = pindex.flatten()
+ rho.sort()
+ rho = np.unique(rho,return_index=True,return_inverse=False)[1]
+ rho = np.append(rho[1:]-rho[:-1],[np.prod(pindex.shape)-rho[-1]])
+
+ return kindex,rho,pindex
+
+
def bin_power_indices(pindex,kindex,rho,log=False,nbin=None,binbounds=None):
"""
Returns the (re)binned power indices associated with the Fourier grid.
@@ -622,6 +682,31 @@ def nkdict_fast(axes,dgrid,fourier=True):
return np.sqrt(np.sum((temp_vecs),axis=-1))
+def nkdict_fast2(axes,dgrid,fourier=True):
+ """
+ Calculates an n-dimensional array with its entries being the lengths of
+ the k-vectors from the zero point of the grid.
+
+ """
+ if(fourier):
+ dk = dgrid
+ else:
+ dk = np.array([1/dgrid[i]/axes[i] for i in range(len(axes))])
+
+ inds = []
+ for a in axes:
+ inds += [slice(0,a)]
+
+ cords = np.ogrid[inds]
+
+ dists = ((cords[0]-axes[0]//2)*dk[0])**2
+ for ii in range(1,len(axes)):
+ dists = dists + ((cords[ii]-axes[ii]//2)*dk[ii])**2
+ dists = np.sqrt(dists)
+
+ return dists
+
+
def nklength(kdict):
return np.sort(list(set(kdict.flatten())))
diff --git a/setup.py b/setup.py
index 21175525c87332b48aee9f24c75461e5369eeb56..f2ea738b5272120b4dd1a49a3434a223dd1a9c3c 100644
--- a/setup.py
+++ b/setup.py
@@ -23,17 +23,12 @@ from distutils.core import setup
import os
setup(name="nifty",
- version="0.8.0",
+ version="0.9.0",
description="Numerical Information Field Theory",
author="Marco Selig",
author_email="mselig@mpa-garching.mpg.de",
url="http://www.mpa-garching.mpg.de/ift/nifty/",
- packages=["nifty"],
+ packages=["nifty", "nifty.demos"],
package_dir={"nifty": ""},
- package_data={"nifty": ["demos/demo_excaliwir.py",
- "demos/demo_faraday.py",
- "demos/demo_faraday_map.npy",
- "demos/demo_wf1.py",
- "demos/demo_wf2.py",
- "demos/demo_wf3.py"]},
- data_files=[(os.path.expanduser('~') + "/.nifty", ["nifty_config"])])
\ No newline at end of file
+ data_files=[(os.path.expanduser('~') + "/.nifty", ["nifty_config"])])
+