Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
212 changes: 106 additions & 106 deletions montepython/likelihoods/ska1_lensing/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,14 +37,14 @@ def __init__(self, path, data, command_line):
if (self.use_halofit):
self.need_cosmo_arguments(data, {'non linear':'halofit'})

# Warn if theoretical error and linear cutoff are requested
if (self.use_lmax_lincut and self.theoretical_error!=0):
warnings.warn("A lmax cutoff infered from kmax and a theoretical error are requested. This combination is not implemented and in most cases not necessary as the theoretical error should induce a cutoff by itself.")
# Warn if theoretical error and linear cutoff are requested
if (self.use_lmax_lincut and self.theoretical_error!=0):
warnings.warn("A lmax cutoff infered from kmax and a theoretical error are requested. This combination is not implemented and in most cases not necessary as the theoretical error should induce a cutoff by itself.")

# Define array of l values, and initialize them
# It is a logspace
# find nlmax in order to reach lmax with logarithmic steps dlnl
self.nlmax = np.int(np.log(self.lmax/self.lmin)/self.dlnl)+1
self.nlmax = int(np.log(self.lmax/self.lmin)/self.dlnl)+1
# redefine slightly dlnl so that the last point is always exactly lmax
self.dlnl = np.log(self.lmax/self.lmin)/(self.nlmax-1)
self.l = self.lmin*np.exp(self.dlnl*np.arange(self.nlmax))
Expand All @@ -69,7 +69,7 @@ def __init__(self, path, data, command_line):
# first value is already correctly set to 0.
self.z_bin_edge = np.zeros(self.nbin+1, 'float64')

total_count = 0.
total_count = 0.
for Bin in xrange(self.nbin-1):
bin_count = 0.
z = self.z_bin_edge[Bin]
Expand All @@ -79,23 +79,23 @@ def __init__(self, path, data, command_line):
bin_count += 0.5*(gd_1+gd_2)*self.dz
z += self.dz
self.z_bin_edge[Bin+1] = z
total_count += bin_count
total_count += bin_count
self.z_bin_edge[self.nbin] = self.zmax

# Fill array of discrete z values
self.z = np.linspace(0, self.zmax, num=self.nzmax)

# check normalizations:
# check normalizations:
#print("n_tot = %s" % (n_tot))
#int_nz = self.nzmax
#int_z = np.linspace(0, self.zmax, num=int_nz)
#DDD = np.zeros(int_nz, 'float64')
#for nz in xrange(self.nzmax):
# #for nz2 in xrange(int_nz):
# # DDD[nz2] = self.photo_z_distribution(self.z[nz], int_z[nz2],False)
#int_nz = self.nzmax
#int_z = np.linspace(0, self.zmax, num=int_nz)
#DDD = np.zeros(int_nz, 'float64')
#for nz in xrange(self.nzmax):
# #for nz2 in xrange(int_nz):
# # DDD[nz2] = self.photo_z_distribution(self.z[nz], int_z[nz2],False)
# DDD = self.photo_z_distribution(self.z[nz], int_z)
# D_tot = scipy.integrate.trapz(DDD,int_z)
# print("%s\t%s\t%s" % (self.z[nz], self.galaxy_distribution(self.z[nz]), D_tot ))
# print("%s\t%s\t%s" % (self.z[nz], self.galaxy_distribution(self.z[nz]), D_tot ))

# Fill distribution for each bin (convolving with photo_z distribution)
self.eta_z = np.zeros((self.nzmax, self.nbin), 'float64')
Expand Down Expand Up @@ -157,7 +157,7 @@ def galaxy_distribution(self, z, array=False):
If the array flag is set to True, z is then interpretated as an array,
and not as a single value.

Replaced by Tim Sprenger due to switch from euclid to ska
Replaced by Tim Sprenger due to switch from euclid to ska
"""

z0 = self.par_zm/self.par_a
Expand All @@ -170,58 +170,58 @@ def galaxy_distribution(self, z, array=False):
return galaxy_dist

def photo_z_distribution(self, z, zph, array=True):
"""
Photo z distribution

If the array flag is set to True, zph is then interpretated as an array,
and not as a single value.

Replaced by Tim Sprenger due to switch from euclid to ska
"""

# Note: you must normalize it yourself to one if you want to get nice
# plots of the galaxy distribution function in each bin (otherwise, the
# spectra will remain correct, but each D_i(x) will loot strangely
# normalized when compared to the original D(z)
if not array:
if (z<=self.z_spec):
photo_z_dist = (1.-self.f_spec)*math.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
if (z==zph):
photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
if (z==0.):
photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
elif (z<=self.z_phot):
photo_z_dist = math.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
else:
photo_z_dist = math.exp(
-0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
1.+z)/math.sqrt(2.*math.pi)
else:
if (z<=self.z_spec):
photo_z_dist = (1.-self.f_spec)*np.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
for index_z in xrange(len(zph)):
if (z==zph[index_z]):
photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
if (z==0.):
photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
break
elif (z<=self.z_phot):
photo_z_dist = np.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
else:
photo_z_dist = np.exp(
-0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
1.+z)/math.sqrt(2.*math.pi)


return photo_z_dist
"""
Photo z distribution

If the array flag is set to True, zph is then interpretated as an array,
and not as a single value.

Replaced by Tim Sprenger due to switch from euclid to ska
"""

# Note: you must normalize it yourself to one if you want to get nice
# plots of the galaxy distribution function in each bin (otherwise, the
# spectra will remain correct, but each D_i(x) will loot strangely
# normalized when compared to the original D(z)
if not array:
if (z<=self.z_spec):
photo_z_dist = (1.-self.f_spec)*math.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
if (z==zph):
photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
if (z==0.):
photo_z_dist += self.f_spec/self.zmax*(self.nzmax-1.)
elif (z<=self.z_phot):
photo_z_dist = math.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
else:
photo_z_dist = math.exp(
-0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
1.+z)/math.sqrt(2.*math.pi)
else:
if (z<=self.z_spec):
photo_z_dist = (1.-self.f_spec)*np.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
for index_z in xrange(len(zph)):
if (z==zph[index_z]):
photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
if (z==0.):
photo_z_dist[index_z] += self.f_spec/self.zmax*(self.nzmax-1.)
break
elif (z<=self.z_phot):
photo_z_dist = np.exp(
-0.5*((z-zph)/self.sigma_phot/(1.+z))**2)/self.sigma_phot/(
1.+z)/math.sqrt(2.*math.pi)
else:
photo_z_dist = np.exp(
-0.5*((z-zph)/self.sigma_noz/(1.+z))**2)/self.sigma_noz/(
1.+z)/math.sqrt(2.*math.pi)


return photo_z_dist

def loglkl(self, cosmo, data):

Expand Down Expand Up @@ -251,9 +251,9 @@ def loglkl(self, cosmo, data):
fun[1:]+fun[:-1])*(self.r[nr+1:]-self.r[nr:-1]))
g[nr, Bin] *= 2.*self.r[nr]*(1.+self.z[nr])

# compute the maximum l where most contributions are linear
# as a function of the lower bin number
if self.use_lmax_lincut:
# compute the maximum l where most contributions are linear
# as a function of the lower bin number
if self.use_lmax_lincut:
lintegrand_lincut_o = np.zeros((self.nzmax, self.nbin, self.nbin), 'float64')
lintegrand_lincut_u = np.zeros((self.nzmax, self.nbin, self.nbin), 'float64')
l_lincut = np.zeros((self.nbin, self.nbin), 'float64')
Expand All @@ -276,26 +276,26 @@ def loglkl(self, cosmo, data):
lintegrand_lincut_u[1:, Bin1, Bin2] +
lintegrand_lincut_u[:-1, Bin1, Bin2])*(
self.r[1:]-self.r[:-1]))
z_peak = np.zeros((self.nbin, self.nbin), 'float64')
z_peak = np.zeros((self.nbin, self.nbin), 'float64')
for Bin1 in xrange(self.nbin):
for Bin2 in xrange(Bin1,self.nbin):
z_peak[Bin1, Bin2] = self.zmax
z_peak[Bin1, Bin2] = self.zmax
for index_z in xrange(self.nzmax):
if (self.r[index_z]>l_lincut[Bin1, Bin2]):
z_peak[Bin1, Bin2] = self.z[index_z]
break
if (self.r[index_z]>l_lincut[Bin1, Bin2]):
z_peak[Bin1, Bin2] = self.z[index_z]
break
if self.use_zscaling:
l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()*pow(1.+z_peak[Bin1, Bin2],2./(2.+cosmo.n_s()))
l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()*pow(1.+z_peak[Bin1, Bin2],2./(2.+cosmo.n_s()))
else:
l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()
l_lincut_mean[Bin1] = np.sum(l_lincut[Bin1, :])/(self.nbin-Bin1)
l_lincut[Bin1, Bin2] *= self.kmax_hMpc*cosmo.h()
l_lincut_mean[Bin1] = np.sum(l_lincut[Bin1, :])/(self.nbin-Bin1)

#for Bin1 in xrange(self.nbin):
#for Bin2 in xrange(Bin1,self.nbin):
#print("%s\t%s\t%s\t%s" % (Bin1, Bin2, l_lincut[Bin1, Bin2], l_lincut_mean[Bin1]))
#for Bin1 in xrange(self.nbin):
#for Bin2 in xrange(Bin1,self.nbin):
#print("%s\t%s\t%s\t%s" % (Bin1, Bin2, l_lincut[Bin1, Bin2], l_lincut_mean[Bin1]))

#for nr in xrange(1, self.nzmax-1):
# print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.z[nr], g[nr, 0], g[nr, 1], g[nr, 2], g[nr, 3], g[nr, 4], g[nr, 5], g[nr, 6], g[nr, 7], g[nr, 8], g[nr, 9]))
#for nr in xrange(1, self.nzmax-1):
# print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (self.z[nr], g[nr, 0], g[nr, 1], g[nr, 2], g[nr, 3], g[nr, 4], g[nr, 5], g[nr, 6], g[nr, 7], g[nr, 8], g[nr, 9]))

# Get power spectrum P(k=l/r,z(r)) from cosmological module
kmin_in_inv_Mpc = self.k_min_h_by_Mpc * cosmo.h()
Expand All @@ -318,7 +318,7 @@ def loglkl(self, cosmo, data):
else:
pk[index_l, index_z] = cosmo.pk(self.l[index_l]/self.r[index_z], self.z[index_z])

#print("%s\t%s\t%s" %(self.l[index_l], self.z[index_z], pk[index_l, index_z]))
#print("%s\t%s\t%s" %(self.l[index_l], self.z[index_z], pk[index_l, index_z]))

# Recover the non_linear scale computed by halofit. If no scale was
# affected, set the scale to one, and make sure that the nuisance
Expand All @@ -329,15 +329,15 @@ def loglkl(self, cosmo, data):
else:
k_sigma = cosmo.nonlinear_scale(self.z, self.nzmax)

if not (cosmo.nonlinear_method == 0):
k_sigma_problem = False
for index_z in xrange(self.nzmax-1):
if (k_sigma[index_z+1]<k_sigma[index_z]) or (k_sigma[index_z+1]>2.5):
k_sigma[index_z+1] = 2.5
k_sigma_problem = True
#print("%s\t%s" % (k_sigma[index_z], self.z[index_z]))
if k_sigma_problem:
warnings.warn("There were unphysical (decreasing in redshift or exploding) values of k_sigma (=cosmo.nonlinear_scale(...)). To proceed they were set to 2.5, the highest scale that seems to be stable.")
if not (cosmo.nonlinear_method == 0):
k_sigma_problem = False
for index_z in xrange(self.nzmax-1):
if (k_sigma[index_z+1]<k_sigma[index_z]) or (k_sigma[index_z+1]>2.5):
k_sigma[index_z+1] = 2.5
k_sigma_problem = True
#print("%s\t%s" % (k_sigma[index_z], self.z[index_z]))
if k_sigma_problem:
warnings.warn("There were unphysical (decreasing in redshift or exploding) values of k_sigma (=cosmo.nonlinear_scale(...)). To proceed they were set to 2.5, the highest scale that seems to be stable.")

# Define the alpha function, that will characterize the theoretical
# uncertainty. Chosen to be 0.001 at low k, raise between 0.1 and 0.2
Expand Down Expand Up @@ -504,7 +504,7 @@ def loglkl(self, cosmo, data):
# TODO parallelize this
for index, ell in enumerate(ells):

if self.use_lmax_lincut:
if self.use_lmax_lincut:
CutBin = -1
for zBin in xrange(self.nbin):
if (ell<l_lincut_mean[zBin]):
Expand All @@ -514,7 +514,7 @@ def loglkl(self, cosmo, data):
break
if (CutBin==-1):
break
else:
else:
det_theory = np.linalg.det(Cov_theory[index, :, :])
det_observ = np.linalg.det(Cov_observ[index, :, :])

Expand Down Expand Up @@ -549,14 +549,14 @@ def loglkl(self, cosmo, data):
newCov[:, i] = Cov_observ[index,:, i]#MArchi#newCov[:, i] = Cov_observ[:, i]
det_theory_plus_error_cross_obs += np.linalg.det(
newCov)
try:
function_vector[k] = (2.*ell+1.)*self.fsky*(det_theory_plus_error_cross_obs/det_theory_plus_error + math.log(det_theory_plus_error/det_observ) - self.nbin ) + dof*vector[k]**2
except ValueError:
warnings.warn("ska1_lensing: Could not evaluate chi2 including theoretical error with the current parameters. The corresponding chi2 is now set to nan!")
break
break
break
chi2 = np.nan
try:
function_vector[k] = (2.*ell+1.)*self.fsky*(det_theory_plus_error_cross_obs/det_theory_plus_error + math.log(det_theory_plus_error/det_observ) - self.nbin ) + dof*vector[k]**2
except ValueError:
warnings.warn("ska1_lensing: Could not evaluate chi2 including theoretical error with the current parameters. The corresponding chi2 is now set to nan!")
break
break
break
chi2 = np.nan


# Computing first
Expand All @@ -582,20 +582,20 @@ def loglkl(self, cosmo, data):


else:
if self.use_lmax_lincut:
if self.use_lmax_lincut:
det_cross = 0.
for i in xrange(0,self.nbin-CutBin):
newCov = np.copy(Cov_theory[index, CutBin:, CutBin:])
newCov[:, i] = Cov_observ[index, CutBin:, CutBin+i]
det_cross += np.linalg.det(newCov)
else:
else:
det_cross = 0.
for i in xrange(self.nbin):
newCov = np.copy(Cov_theory[index, :, :])
newCov[:, i] = Cov_observ[index, :, i]
det_cross += np.linalg.det(newCov)

if self.use_lmax_lincut:
if self.use_lmax_lincut:
chi2 += (2.*ell+1.)*self.fsky*(det_cross/det_theory + math.log(det_theory/det_observ) - self.nbin+CutBin)
else:
chi2 += (2.*ell+1.)*self.fsky*(det_cross/det_theory + math.log(det_theory/det_observ) - self.nbin)
Expand Down
Loading