from __future__ import division
import types
import numpy as np
from scipy.integrate import odeint
from astropy.cosmology import w0waCDM
from astropy.units import rad
try:
from ..extern import _nicaea
_nicaea=_nicaea
except ImportError:
_nicaea=None
##########################################################
##########Nicaea built-in redshift distributions##########
##########################################################
_nicaea_builtin = dict()
_nicaea_builtin["ludo"] = 5
_nicaea_builtin["jonben"] = 5
_nicaea_builtin["ymmk"] = 5
_nicaea_builtin["single"] = 2
_nicaea_builtin["hist"] = 0
################################################################################
###########Useful to check validity of redshift distribution inputs#############
################################################################################
def _check_redshift(z,distribution,distribution_parameters,**kwargs):
#Parse redshift distribution from input
if type(z)==np.float:
nzbins = 1
nofz = ["single"]
Nnz = np.array([2],dtype=np.int32)
par_nz = np.array([z,z])
elif type(distribution)==types.FunctionType:
assert type(z)==np.ndarray and z.ndim==1,"distribution is a callable, hence z must be a one dimensional array!"
if distribution_parameters is None:
distribution_parameters="one"
assert distribution_parameters in ["one","all"],"distribution is a callable, hence distribution_parameters must be either 'one' or 'all'"
if distribution_parameters=="one":
nzbins = 1
nofz = ["hist"]
Nnz = np.array([2*len(z)-1],dtype=np.int32)
#Format paramaters accordingly
midpoints = 0.5*(z[1:]+z[:-1])
gal_in_bin = distribution(midpoints,**kwargs)
par_nz = np.concatenate((np.array([z[0],z[-1]]),z[1:-1],gal_in_bin))
assert len(par_nz)==Nnz[0]
else:
nzbins = len(z)-1
nofz = ["hist"]*nzbins
Nnz = np.ones(nzbins,dtype=np.int32)*3
#Format paramaters accordingly
midpoints = 0.5*(z[1:]+z[:-1])
gal_in_bin = distribution(midpoints,**kwargs)
par_nz = np.zeros(3*nzbins)
for n in range(nzbins):
par_nz[3*n] = z[n]
par_nz[3*n+1] = z[n+1]
par_nz[3*n+2] = gal_in_bin[n]
else:
assert z is None,"If you want to specify the distribution parameters manually z must be None"
assert type(distribution)==list and type(distribution_parameters)==list,"If you want to specify the distribution parameters manually both distribution and distribution parameters must be lists"
assert len(distribution)==len(distribution_parameters),"The lists must have the same length"
nzbins = len(distribution)
nofz = distribution
Nnz = list()
for n,distr_in_bin in enumerate(nofz):
assert distr_in_bin in _nicaea_builtin.keys(),"{0} must be one of NICAEA built in types: {1}".format(distr_in_bin,_nicaea_builtin.keys())
if distr_in_bin=="hist":
#Check that number of parameters is odd to conform to standard
assert len(distribution_parameters[n]%2==1)
Nnz.append(len(distribution_parameters[n]))
else:
#Check that number of parameters conforms to standard
assert len(distribution_parameters[n])==_nicaea_builtin[distr_in_bin],"You must provide exactly {0} parameters for distribution {1}".format(_nicaea_builtin[distr_in_bin],distr_in_bin)
Nnz.append(len(distribution_parameters[n]))
#Build the ordered array par_nz
assert len(Nnz)==nzbins
par_nz = np.concatenate(distribution_parameters)
#Return in NICAEA friendly format
return nzbins,nofz,Nnz,par_nz
##################################################################
######Useful for integrating the linear growth factor ODE#########
##################################################################
def _delta_dot(delta,z,cosmology):
#Densities
Ode = cosmology.Ode(z)
w = cosmology.w(z)
#Inhomogeneous term in the ODE
inhomogeneous = np.array([0.0,-1.5 * Ode * (1+3*w) / (1+z)**2])
#Return
return _jacobian_delta_dot(delta,z,cosmology).dot(delta) + inhomogeneous
def _jacobian_delta_dot(delta,z,cosmology):
#Densities
Om = cosmology.Om(z)
Ode = cosmology.Ode(z)
w = cosmology.w(z)
#Linear term in the ODE
linear21 = 1.5*Om / (1+z)**2
linear22 = -0.5 * (Om+Ode*(1+3*w)) / (1+z)
#Return the jacobian
return np.array([[0.0,1.0],[linear21,linear22]])
##########################################
##########NicaeaSettings class############
##########################################
[docs]class NicaeaSettings(dict):
"""
Class handler of the code settings (non linear modeling, tomography, transfer function, etc...)
"""
def __init__(self):
super(NicaeaSettings,self).__init__()
#Default settings
self["snonlinear"]="smith03_revised"
self["stransfer"]="eisenhu"
self["sgrowth"]="growth_de"
self["sde_param"]="linder"
self["normmode"]="norm_s8"
self["stomo"]="tomo_all"
self["sreduced"]="none"
self["q_mag_size"]=1.0
@classmethod
[docs] def default(cls):
"""
Generate default settings
:returns: NicaeaSettings defaults instance
"""
return cls()
@property
def knobs(self):
"""
Lists available settings to tune
"""
return self.keys()
[docs] def available(self,knob):
"""
Given a settings, lists all the possible values
"""
#Available settings
if knob=="snonlinear":
return [ "linear", "pd96", "smith03", "smith03_de", "coyote10", "coyote13", "halodm", "smith03_revised" ]
elif knob=="stransfer":
return [ "bbks", "eisenhu", "eisenhu_osc", "be84" ]
elif knob=="sgrowth":
return [ "heath", "growth_de" ]
elif knob=="sde_param":
return [ "jassal", "linder", "earlyDE", "poly_DE" ]
elif knob=="normmode":
return [ "norm_s8" , "norm_as" ]
elif knob=="stomo":
return [ "tomo_all", "tomo_auto_only", "tomo_cross_only" ]
elif knob=="sreduced":
return ["none", "reduced_K10"]
elif knob=="q_mag_size":
print("Positive float")
return None
else:
raise ValueError("{0} is not a tunable setting!".format(knob))
##########################################
##########Nicaea class####################
##########################################
[docs]class Nicaea(w0waCDM):
"""
Main class handler for the python bindings of the NICAEA cosmological code, written by M. Kilbinger & collaborators
"""
def __init__(self,H0=72.0,Om0=0.26,Ode0=0.74,Ob0=0.046,w0=-1.0,wa=0.0,sigma8=0.800,ns=0.960,name=None):
super(Nicaea,self).__init__(H0,Om0,Ode0,w0=w0,wa=wa,Ob0=Ob0,name=name)
self.sigma8=sigma8
self.ns=ns
def __repr__(self):
astropy_string = super(Nicaea,self).__repr__()
pieces = astropy_string.split(",")
si8_piece = u" sigma8={0}".format(self.sigma8)
ns_piece = u" ns={0}".format(self.ns)
return ",".join(pieces[:3] + [si8_piece,ns_piece] + pieces[3:])
@classmethod
[docs] def fromCosmology(cls,cosmo):
"""
Builds a Nicaea instance from one of astropy.cosmology objects, from which it inherits all the cosmological parameter values
:param cosmo: one of astropy cosmology instances
:type cosmo: astropy FLRW
:returns: Nicaea instance with the cosmological parameters inherited from cosmo
"""
#Get the cosmological parameter values out of the cosmo object
H0 = cosmo.H0.value
Om0 = cosmo.Om0
Ode0 = cosmo.Ode0
#Dark energy
if hasattr(cosmo,"w0"):
w0=cosmo.w0
else:
w0=-1.0
if hasattr(cosmo,"wa"):
wa=cosmo.wa
else:
wa=0.0
#Neutrinos
Neff = cosmo.Neff
Onu0 = cosmo.Onu0
#Set these manually to default
ns = 0.960
sigma8 = 0.800
Ob0 = 0.046
#Instantiate
return cls(H0=H0,Om0=Om0,Ode0=Ode0,Ob0=Ob0,w0=w0,wa=wa,sigma8=sigma8,ns=ns)
################################################################################################################
##################################################################
###########Growth factor and its derivative#######################
##################################################################
def growth_factor(self,z,delta0=np.array([1.0,0.0]),**kwargs):
"""
Computes the linear growth factor at redshift z, given initial conditions on the density contrast
:param z: array of redshifts on which to compute the growth factor
:type z: array.
:param delta0: initial condition (delta0,Ddelta0/dz) for the ODE
:type delta0: array.
:param kwargs: the keyword arguments are passed to odeint
:type kwargs: dict.
:returns: growth factor and its redshift derivative
:rtype: array
"""
return odeint(_delta_dot,y0=delta0,t=z,Dfun=_jacobian_delta_dot,args=(self,),**kwargs) / delta0[0]
################################################################################################################
[docs] def convergencePowerSpectrum(self,ell,z=2.0,distribution=None,distribution_parameters=None,settings=None,**kwargs):
"""
Computes the convergence power spectrum for the given cosmological parameters and redshift distribution using NICAEA
:param ell: multipole moments at which to compute the power spectrum
:type ell: array.
:param z: redshift bins for the sources; if a single float is passed, single redshift is assumed
:type z: float., array or None
:param distribution: redshift distribution of the sources (normalization not necessary); if None a single redshift is expected; if callable, z must be an array, and a single redshift bin is considered, with the galaxy distribution specified by the call of distribution(z); if a list is passed, each element must be a NICAEA type
:type distribution: None,callable or list
:param distribution_parameters: relevant only when distribution is a list or callable. When distribution is callable, distribution_parameters has to be one between "one" and "all" to decide if one or multiple redshift bins have to be considered. If it is a list, each element in it should be the tuple of parameters expected by the correspondent NICAEA distribution type
:type distribution_parameters: str. or list.
:param settings: NICAEA code settings
:type settings: NicaeaSettings instance
:param kwargs: the keyword arguments are passed to the distribution, if callable
:type kwargs: dict.
:returns: ( NlxNz array ) computed power spectrum at the selected multipoles (when computing the cross components these are returned in row major C ordering)
"""
if _nicaea is None:
raise ImportError("You need to install the Nicaea bindings to use this routine! Check your GSL/FFTW3 installations!")
assert isinstance(ell,np.ndarray)
#If no settings provided, use the default ones
if settings is None:
settings=NicaeaSettings.default()
#Check sanity of input
nzbins,nofz,Nnz,par_nz = _check_redshift(z,distribution,distribution_parameters,**kwargs)
#Compute the power spectrum via NICAEA
power_spectrum_nicaea = _nicaea.shearPowerSpectrum(self.Om0,self.Ode0,self.w0,self.wa,self.H0.value/100.0,self.Ob0,self.Onu0,self.Neff,self.sigma8,self.ns,nzbins,ell,Nnz,nofz,par_nz,settings,None)
#Return
if power_spectrum_nicaea.shape[1]==1:
return power_spectrum_nicaea[:,0]
else:
return power_spectrum_nicaea
[docs] def shearTwoPoint(self,theta,z=2.0,distribution=None,distribution_parameters=None,settings=None,kind="+",**kwargs):
"""
Computes the shear two point function for the given cosmological parameters and redshift distribution using NICAEA
:param theta: angles at which to compute the two point function
:type theta: array. with units
:param z: redshift bins for the sources; if a single float is passed, single redshift is assumed
:type z: float., array or None
:param distribution: redshift distribution of the sources (normalization not necessary); if None a single redshift is expected; if callable, z must be an array, and a single redshift bin is considered, with the galaxy distribution specified by the call of distribution(z); if a list is passed, each element must be a NICAEA type
:type distribution: None,callable or list
:param distribution_parameters: relevant only when distribution is a list or callable. When distribution is callable, distribution_parameters has to be one between "one" and "all" to decide if one or multiple redshift bins have to be considered. If it is a list, each element in it should be the tuple of parameters expected by the correspondent NICAEA distribution type
:type distribution_parameters: str. or list.
:param settings: NICAEA code settings
:type settings: NicaeaSettings instance
:param kind: must be "+" or "-"
:type kind: str.
:param kwargs: the keyword arguments are passed to the distribution, if callable
:type kwargs: dict.
:returns: ( NtxNz array ) computed two point function at the selected angles (when computing the cross components these are returned in row major C ordering)
"""
if _nicaea is None:
raise ImportError("You need to install the Nicaea bindings to use this routine! Check your GSL/FFTW3 installations!")
assert isinstance(theta,np.ndarray)
#If no settings provided, use the default ones
if settings is None:
settings=NicaeaSettings.default()
#Convert angles in radians
theta_rad = theta.to(rad).value
#Check sanity of input
nzbins,nofz,Nnz,par_nz = _check_redshift(z,distribution,distribution_parameters,**kwargs)
#Plus or minus?
if kind=="+":
pm = 1
elif kind=="-":
pm = -1
else:
raise ValueError("kind must be either + or -")
#Compute the two point function using NICAEA
two_point_function_nicaea = _nicaea.shear2pt(self.Om0,self.Ode0,self.w0,self.wa,self.H0.value/100.0,self.Ob0,self.Onu0,self.Neff,self.sigma8,self.ns,nzbins,theta_rad,Nnz,nofz,par_nz,settings,pm)
#Return
if two_point_function_nicaea.shape[1]==1:
return two_point_function_nicaea[:,0]
else:
return two_point_function_nicaea