Skip to content

GitLab

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

intermediate update.

parent c2dbf0b4
Hide whitespace changes
Inline Side-by-side
Showing with 231 additions and 71 deletions
__init__.py
View file @ fb3d0c82
......@@ -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 *
##-----------------------------------------------------------------------------
......
nifty_cmaps.py
View file @ fb3d0c82
......@@ -237,4 +237,44 @@ 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)
##-----------------------------------------------------------------------------
nifty_core.py
View file @ fb3d0c82
......@@ -514,7 +514,7 @@ class _about(object): ## nifty support class for global settings
"""
## version
self._version = "0.8.4"
self._version = "0.8.9"
## 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
......
nifty_power.py
View file @ fb3d0c82
......@@ -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)."))
......
nifty_power_conversion.py
View file @ fb3d0c82
## 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.")
......
nifty_tools.py
View file @ fb3d0c82
......@@ -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)
......
powerspectrum.py
View file @ fb3d0c82
......@@ -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())))
......
setup.py
View file @ fb3d0c82
......@@ -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"])])
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