Commit d876195c authored by Marco Selig's avatar Marco Selig
Browse files

version update.

parent 80236a7a
...@@ -96,7 +96,7 @@ Requirements ...@@ -96,7 +96,7 @@ Requirements
Download 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 at `<https://github.com/mselig/nifty/tags>`_. The current version can be
obtained by cloning the repository:: obtained by cloning the repository::
......
...@@ -21,8 +21,9 @@ ...@@ -21,8 +21,9 @@
from __future__ import division from __future__ import division
from nifty_core import * from nifty_core import *
#from nifty_power import * from nifty_cmaps import *
#from nifty_cmaps import * from nifty_power import *
##----------------------------------------------------------------------------- ##-----------------------------------------------------------------------------
......
...@@ -2178,6 +2178,9 @@ class rg_space(space): ...@@ -2178,6 +2178,9 @@ class rg_space(space):
Notes Notes
----- -----
Only even numbers of grid points per axis are supported. 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 Attributes
---------- ----------
......
...@@ -291,7 +291,7 @@ def _calc_inverse(tk,var,kindex,rho,b1,Amem): ## > computes the inverse Hessian ...@@ -291,7 +291,7 @@ def _calc_inverse(tk,var,kindex,rho,b1,Amem): ## > computes the inverse Hessian
## inversion ## inversion
return np.linalg.inv(T2+np.diag(b2,k=0)),b2,Amem 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. Infers the power spectrum.
...@@ -335,9 +335,9 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None ...@@ -335,9 +335,9 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
(default: (1,0)). (default: (1,0)).
smoothness : bool, *optional* smoothness : bool, *optional*
Indicates whether the smoothness prior is used or not Indicates whether the smoothness prior is used or not
(default: False). (default: True).
var : {scalar, list, array}, *optional* 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* bare : bool, *optional*
Indicates whether the power spectrum entries returned are "bare" Indicates whether the power spectrum entries returned are "bare"
or not (mandatory for the correct incorporation of volume weights) 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 ...@@ -503,30 +503,45 @@ def infer_power(m,domain=None,Sk=None,D=None,pindex=None,pundex=None,kindex=None
if(smoothness): if(smoothness):
if(not domain.discrete): if(not domain.discrete):
numerator = weight_power(domain,numerator,power=-1,pindex=pindex,pundex=pundex) ## bare(!) numerator = weight_power(domain,numerator,power=-1,pindex=pindex,pundex=pundex) ## bare(!)
pk = numerator/denominator1 ## bare(!)
## smoothness prior ## smoothness prior
tk = np.log(pk) divergence = 1
Amemory = None while(divergence):
var_ = var*1.1 # temporally increasing the variance pk = numerator/denominator1 ## bare(!)
breakinfo = False tk = np.log(pk)
while(var_>=var): # slowly lowering the variance Amemory = None
absdelta = 1 var_ = var*1.1 # temporally increasing the variance
while(absdelta>1E-3): # solving with fixed variance var_OLD = -1
## solution of A delta = b1 - b2 breakinfo = False
Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory) while(var_>=var): # slowly lowering the variance
delta = np.dot(Ainverse,numerator/pk-denominator2,out=None) absdelta = 1
if(np.abs(delta).max()>absdelta): # increasing variance when speeding up while(absdelta>1E-3): # solving with fixed variance
var_ *= 1.1 ## solution of A delta = b1 - b2
absdelta = np.abs(delta).max() Ainverse,denominator2,Amemory = _calc_inverse(tk,var_,kindex,rho,denominator1,Amemory)
tk += min(1,0.1/absdelta)*delta # adaptive step width delta = np.dot(Ainverse,numerator/pk-denominator2,out=None)
pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width if(np.abs(delta).max()>absdelta): # increasing variance when speeding up
var_ /= 1.1 # lowering the variance when converged var_ *= 1.1
if(var_<var): absdelta = np.abs(delta).max()
if(breakinfo): # making sure there's one iteration with the correct variance tk += min(1,0.1/absdelta)*delta # adaptive step width
break pk *= np.exp(min(1,0.1/absdelta)*delta) # adaptive step width
var_ = var var_ /= 1.1 # lowering the variance when converged
breakinfo = True 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 ... ## weight if ...
if(not domain.discrete)and(not bare): if(not domain.discrete)and(not bare):
......
...@@ -218,11 +218,11 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k ...@@ -218,11 +218,11 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
""" """
## field absolutes ## field absolutes
if(not fourier): 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: else:
foufield = field 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 fieldabs = np.abs(foufield)**2
if(rho is None): if(rho is None):
...@@ -242,6 +242,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k ...@@ -242,6 +242,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
ps[position] += fieldabs[ii] ps[position] += fieldabs[ii]
rho[position] += 1 rho[position] += 1
else: else:
## zerocenter pindex
if(np.any(zerocentered==False)):
pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
## power spectrum ## power spectrum
ps = np.zeros(np.max(pindex)+1) ps = np.zeros(np.max(pindex)+1)
rho = np.zeros(ps.size) rho = np.zeros(ps.size)
...@@ -262,6 +265,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k ...@@ -262,6 +265,9 @@ def calc_ps_fast(field,axes,dgrid,zerocentered=False,fourier=False,pindex=None,k
position = np.searchsorted(klength,kdict[ii]) position = np.searchsorted(klength,kdict[ii])
ps[position] += fieldabs[ii] ps[position] += fieldabs[ii]
else: else:
## zerocenter pindex
if(np.any(zerocentered==False)):
pindex = np.fft.fftshift(pindex, axes=shiftaxes(zerocentered,st_to_zero_mode=True))
## power spectrum ## power spectrum
ps = np.zeros(rho.size) ps = np.zeros(rho.size)
for ii in np.ndindex(pindex.shape): for ii in np.ndindex(pindex.shape):
......
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