diff --git a/README.rst b/README.rst
index 68bedc88411bb4438611cf449c007ce768073cd0..09ed93af3f96b204ba8a4790d5b5cf62ba9db807 100644
--- a/README.rst
+++ b/README.rst
@@ -96,7 +96,7 @@ Requirements
Download
........
-The latest release is tagged **v0.4.0** and is available as a source package
+The latest release is tagged **v0.4.2** 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 561e8c65ee21e6ce21f6c87701df1259b1c8f5ab..deb66f085e0fd160477bc862d01de1f94443409f 100644
--- a/__init__.py
+++ b/__init__.py
@@ -21,8 +21,9 @@
from __future__ import division
from nifty_core import *
-#from nifty_power import *
-#from nifty_cmaps import *
+from nifty_cmaps import *
+from nifty_power import *
+
##-----------------------------------------------------------------------------
diff --git a/nifty_core.py b/nifty_core.py
index ff4b72821ec20e4db6e578ea8c109567fe11176d..3efbc4eb22d49e416a65533e788c6e88420769ea 100644
--- a/nifty_core.py
+++ b/nifty_core.py
@@ -2178,6 +2178,9 @@ class rg_space(space):
Notes
-----
Only even numbers of grid points per axis are supported.
+ The basis transformations between position `x` and Fourier mode `k`
+ rely on (inverse) fast Fourier transformations using the
+ :math:`exp(2 \pi i k^\dagger x)`-formulation.
Attributes
----------
diff --git a/nifty_power.py b/nifty_power.py
index 0f967e5a7592bc0ab7b1aafa32d889eca797688c..2c50205b73ba0df87442bcc9357c5a4ca0935fb6 100644
--- a/nifty_power.py
+++ b/nifty_power.py
@@ -291,7 +291,7 @@ def _calc_inverse(tk,var,kindex,rho,b1,Amem): ## > computes the inverse Hessian
## inversion
return np.linalg.inv(T2+np.diag(b2,k=0)),b2,Amem
-def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None,rho=None,q=1E-42,alpha=1,perception=(1,0),smoothness=False,var=100,bare=True,**kwargs):
+def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None,rho=None,q=1E-42,alpha=1,perception=(1,0),smoothness=True,var=10,bare=True,**kwargs):
"""
Infers the power spectrum.
@@ -335,9 +335,9 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
(default: (1,0)).
smoothness : bool, *optional*
Indicates whether the smoothness prior is used or not
- (default: False).
+ (default: True).
var : {scalar, list, array}, *optional*
- Variance of the assumed spectral smoothness prior (default: 100).
+ Variance of the assumed spectral smoothness prior (default: 10).
bare : bool, *optional*
Indicates whether the power spectrum entries returned are "bare"
or not (mandatory for the correct incorporation of volume weights)
@@ -503,30 +503,45 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
if(smoothness):
if(not domain.discrete):
numerator = weight_power(domain,numerator,power=-1,pindex=pindex,pundex=pundex) ## bare(!)
- pk = numerator/denominator1 ## bare(!)
## smoothness prior
- tk = np.log(pk)
- Amemory = None
- var_ = var*1.1 # temporally increasing the variance
- breakinfo = False
- while(var_>=var): # slowly lowering the variance
- absdelta = 1
- while(absdelta>1E-3): # solving with fixed variance
- ## solution of A delta = b1 - b2
- Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory)
- delta = np.dot(Ainverse,numerator/pk-denominator2,out=None)
- if(np.abs(delta).max()>absdelta): # increasing variance when speeding up
- var_ *= 1.1
- absdelta = np.abs(delta).max()
- tk += min(1,0.1/absdelta)*delta # adaptive step width
- pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width
- var_ /= 1.1 # lowering the variance when converged
- if(var_<var):
- if(breakinfo): # making sure there's one iteration with the correct variance
- break
- var_ = var
- breakinfo = True
+ divergence = 1
+ while(divergence):
+ pk = numerator/denominator1 ## bare(!)
+ tk = np.log(pk)
+ Amemory = None
+ var_ = var*1.1 # temporally increasing the variance
+ var_OLD = -1
+ breakinfo = False
+ while(var_>=var): # slowly lowering the variance
+ absdelta = 1
+ while(absdelta>1E-3): # solving with fixed variance
+ ## solution of A delta = b1 - b2
+ Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory)
+ delta = np.dot(Ainverse,numerator/pk-denominator2,out=None)
+ if(np.abs(delta).max()>absdelta): # increasing variance when speeding up
+ var_ *= 1.1
+ absdelta = np.abs(delta).max()
+ tk += min(1,0.1/absdelta)*delta # adaptive step width
+ pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width
+ var_ /= 1.1 # lowering the variance when converged
+ if(var_<var):
+ if(breakinfo): # making sure there's one iteration with the correct variance
+ break
+ var_ = var
+ breakinfo = True
+ elif(var_==var_OLD):
+ if(divergence==3):
+ pot = int(np.log10(var_))
+ var = int(1+var_*10**-pot)*10**pot
+ about.warnings.cprint("WARNING: smoothness variance increased ( var = "+str(var)+" ).")
+ break
+ else:
+ divergence += 1
+ else:
+ var_OLD = var_
+ if(breakinfo):
+ break
## weight if ...
if(not domain.discrete)and(not bare):
diff --git a/powerspectrum.py b/powerspectrum.py
index 73e86b322bca3c209643600be87185bfaae83f7e..1d051b294848e76d71a2361a09ac8e672148374d 100644
--- a/powerspectrum.py
+++ b/powerspectrum.py
@@ -218,11 +218,11 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
"""
## field absolutes
if(not fourier):
- foufield = np.fft.fftn(field)
+ 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
- if(np.any(zerocentered==False))and(pindex is None): ## kdict is zerocentered, but pindex is shifted
- foufield = np.fft.fftshift(foufield, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
fieldabs = np.abs(foufield)**2
if(rho is None):
@@ -242,6 +242,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
ps[position] += fieldabs[ii]
rho[position] += 1
else:
+ ## zerocenter pindex
+ if(np.any(zerocentered==False)):
+ pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
## power spectrum
ps = np.zeros(np.max(pindex)+1)
rho = np.zeros(ps.size)
@@ -262,6 +265,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
position = np.searchsorted(klength,kdict[ii])
ps[position] += fieldabs[ii]
else:
+ ## zerocenter pindex
+ if(np.any(zerocentered==False)):
+ pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
## power spectrum
ps = np.zeros(rho.size)
for ii in np.ndindex(pindex.shape):