From 79cfe36f298d35498632001969a6254ea2002d51 Mon Sep 17 00:00:00 2001 From: Zhi-ChaoZhao <73989564+Zhi-ChaoZhao@users.noreply.github.com> Date: Thu, 7 Aug 2025 13:18:05 +0800 Subject: [PATCH 1/2] Fix inconsistent use of tabs and spaces in indentation --- .../likelihoods/ska1_lensing/__init__.py | 212 +++++++++--------- 1 file changed, 106 insertions(+), 106 deletions(-) diff --git a/montepython/likelihoods/ska1_lensing/__init__.py b/montepython/likelihoods/ska1_lensing/__init__.py index ef1cf924..6639c940 100644 --- a/montepython/likelihoods/ska1_lensing/__init__.py +++ b/montepython/likelihoods/ska1_lensing/__init__.py @@ -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)) @@ -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] @@ -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') @@ -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 @@ -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): @@ -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') @@ -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() @@ -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 @@ -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]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]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 @@ -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 Date: Thu, 7 Aug 2025 14:39:41 +0800 Subject: [PATCH 2/2] Fix inconsistent use of tabs and spaces in indentation --- montepython/likelihoods/ska1_pk/__init__.py | 279 ++++++++++---------- 1 file changed, 140 insertions(+), 139 deletions(-) diff --git a/montepython/likelihoods/ska1_pk/__init__.py b/montepython/likelihoods/ska1_pk/__init__.py index c7d5082f..15c204ff 100644 --- a/montepython/likelihoods/ska1_pk/__init__.py +++ b/montepython/likelihoods/ska1_pk/__init__.py @@ -4,6 +4,7 @@ import warnings from numpy import newaxis as na from math import exp, log, pi, log10 +from montepython import io_mp import scipy.integrate # Created by Tim Sprenger in 2017 based on euclid_pk @@ -23,17 +24,17 @@ def __init__(self, path, data, command_line): self.need_cosmo_arguments(data, {'P_k_max_1/Mpc': 1.5*self.k_cut(self.zmax)}) # Compute non-linear power spectrum if requested - if self.use_halofit: - self.need_cosmo_arguments(data, {'non linear':'halofit'}) - print("Using halofit") + if self.use_halofit: + self.need_cosmo_arguments(data, {'non linear':'halofit'}) + print("Using halofit") # Deduce the dz step from the number of bins and the edge values of z self.dz = (self.zmax-self.zmin)/self.nbin - # Compute new zmin and zmax which are bin centers - # Need to be defined as edges if zmin can be close to z=0 - self.zmin += self.dz/2. - self.zmax -= self.dz/2. + # Compute new zmin and zmax which are bin centers + # Need to be defined as edges if zmin can be close to z=0 + self.zmin += self.dz/2. + self.zmax -= self.dz/2. # self.z_mean will contain the central values self.z_mean = np.linspace(self.zmin, self.zmax, num=self.nbin) @@ -45,21 +46,21 @@ def __init__(self, path, data, command_line): # Compute the number of galaxies for each \bar z # N_g(\bar z) = int_{\bar z - dz/2}^{\bar z + dz/2} dn/dz dz - # dn/dz fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z - # dn/dz = number of galaxies per redshift and deg^2 - # old approx.: self.N_g = pow(10.,self.c1)*pow(self.z_mean,self.c2)*np.exp(-self.c3*self.z_mean) * self.dz * self.skycov - self.N_g = np.zeros(self.nbin) - N_tot = 0.0 - for index_z in xrange(self.nbin): - self.N_g[index_z], error = scipy.integrate.quad(self.dndz, self.z_mean[index_z]-self.dz/2., self.z_mean[index_z]+self.dz/2.) - assert error/self.N_g[index_z] <= 0.001, ("dndz integration error is bigger than 0.1%") - N_tot += self.N_g[index_z] - - # Ntot output - #print("\nSKA1: Number of detected galaxies and bias in each redshift bin:") - #for index_z in xrange(self.nbin): - # print("z-bin[" + str(self.z_mean[index_z]-self.dz/2.) + "," + str(self.z_mean[index_z]+self.dz/2.) + "]: \tN = %.4g" % (self.N_g[index_z]) + " ,\t b = %.4g" % (b[index_z])) - #print("Total number of detected galaxies: N = %.4g\n" % (N_tot)) + # dn/dz fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z + # dn/dz = number of galaxies per redshift and deg^2 + # old approx.: self.N_g = pow(10.,self.c1)*pow(self.z_mean,self.c2)*np.exp(-self.c3*self.z_mean) * self.dz * self.skycov + self.N_g = np.zeros(self.nbin) + N_tot = 0.0 + for index_z in xrange(self.nbin): + self.N_g[index_z], error = scipy.integrate.quad(self.dndz, self.z_mean[index_z]-self.dz/2., self.z_mean[index_z]+self.dz/2.) + assert error/self.N_g[index_z] <= 0.001, ("dndz integration error is bigger than 0.1%") + N_tot += self.N_g[index_z] + + # Ntot output + #print("\nSKA1: Number of detected galaxies and bias in each redshift bin:") + #for index_z in xrange(self.nbin): + # print("z-bin[" + str(self.z_mean[index_z]-self.dz/2.) + "," + str(self.z_mean[index_z]+self.dz/2.) + "]: \tN = %.4g" % (self.N_g[index_z]) + " ,\t b = %.4g" % (b[index_z])) + #print("Total number of detected galaxies: N = %.4g\n" % (N_tot)) # Define the k values for the integration (from kmin to kmax), at which # the spectrum will be computed (and stored for the fiducial model) @@ -72,7 +73,7 @@ def __init__(self, path, data, command_line): # If the file exists, initialize the fiducial values, the spectrum will # be read first, with k_size values of k and nbin values of z. Then, # H_fid and D_A fid will be read (each with nbin values). - # Then V_fid, b_fid and the fiducial errors on real space coordinates follow. + # Then V_fid, b_fid and the fiducial errors on real space coordinates follow. self.fid_values_exist = False self.pk_nl_fid = np.zeros((self.k_size, 2*self.nbin+1), 'float64') if self.use_linear_rsd: @@ -81,8 +82,8 @@ def __init__(self, path, data, command_line): self.D_A_fid = np.zeros(2*self.nbin+1, 'float64') self.V_fid = np.zeros(self.nbin, 'float64') self.b_fid = np.zeros(self.nbin, 'float64') - self.sigma_A_fid = np.zeros(self.nbin, 'float64') - self.sigma_B_fid = np.zeros(self.nbin, 'float64') + self.sigma_A_fid = np.zeros(self.nbin, 'float64') + self.sigma_B_fid = np.zeros(self.nbin, 'float64') fid_file_path = os.path.join(self.data_directory, self.fiducial_file) if os.path.exists(fid_file_path): @@ -113,22 +114,22 @@ def __init__(self, path, data, command_line): self.sigma_A_fid[index_z] = float(line.split()[0]) self.sigma_B_fid[index_z] = float(line.split()[1]) line = fid_file.readline() - self.sigma_NL_fid = float(line) + self.sigma_NL_fid = float(line) # Else the file will be created in the loglkl() function. return def dndz(self, z): - # dn/dz/deg^2 fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z * sky coverage in deg^2 - dndz = pow(10.,self.c1)*pow(z,self.c2)*np.exp(-self.c3*z)*self.skycov - return dndz + # dn/dz/deg^2 fit formula from 1412.4700v2: 10^c1*z^c2*e^-c3z * sky coverage in deg^2 + dndz = pow(10.,self.c1)*pow(z,self.c2)*np.exp(-self.c3*z)*self.skycov + return dndz def k_cut(self, z,h=0.6693,n_s=0.9619): - kcut = self.kmax*h - # compute kmax according to highest redshift linear cutoff (1509.07562v2) - if self.use_zscaling: - kcut *= pow(1.+z,2./(2.+n_s)) - return kcut + kcut = self.kmax*h + # compute kmax according to highest redshift linear cutoff (1509.07562v2) + if self.use_zscaling: + kcut *= pow(1.+z,2./(2.+n_s)) + return kcut def loglkl(self, cosmo, data): @@ -145,34 +146,34 @@ def loglkl(self, cosmo, data): # Compute V_survey, for each given redshift bin, which is the volume of # a shell times the sky coverage (only fiducial needed): - if self.fid_values_exist is False: - V_survey = np.zeros(self.nbin, 'float64') - for index_z in xrange(self.nbin): - V_survey[index_z] = 4./3.*pi*self.skycov/41253.*( - r[2*index_z+2]**3-r[2*index_z]**3) + if self.fid_values_exist is False: + V_survey = np.zeros(self.nbin, 'float64') + for index_z in xrange(self.nbin): + V_survey[index_z] = 4./3.*pi*self.skycov/41253.*( + r[2*index_z+2]**3-r[2*index_z]**3) # At the center of each bin, compute the bias function, # using fit fromula from 1412.4700v2: c4*e^c5z - if 'beta_0^SKA1' in self.use_nuisance: - b = self.c4*np.exp(self.c5*self.z_mean*data.mcmc_parameters['beta_1^SKA1']['current']*data.mcmc_parameters['beta_1^SKA1']['scale'])*data.mcmc_parameters['beta_0^SKA1']['current']*data.mcmc_parameters['beta_0^SKA1']['scale'] - else: - b = self.c4*np.exp(self.c5*self.z_mean) - - # Compute freq.res. sigma_r = (1+z)^2/H*sigma_nu/nu_21cm, nu in Mhz - # Compute ang.res. sigma_perp = 1.22(1+z)^2*D_A*lambda_21cm/Baseline, Baseline in km - # combine into exp(-k^2*(mu^2*(sig_r^2-sig_perp^2)+sig_perp^2)) independent of cosmo - # used as exp(-k^2*(mu^2*sigma_A+sigma_B)) all fiducial - if self.fid_values_exist is False: - sigma_A = np.zeros(self.nbin,'float64') - sigma_B = np.zeros(self.nbin,'float64') - sigma_A = ((1.+self.z_mean[:])**2/H[1::2]*self.delta_nu/np.sqrt(8.*np.log(2.))/self.nu0)**2 -( - 1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2 - sigma_B = (1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2 - - # sigma_NL in Mpc = nonlinear dispersion scale of RSD (1405.1452v2) - sigma_NL = 0.0 # fiducial would be 7 but when kept constant that is more constraining than keeping 0 - if 'sigma_NL' in self.use_nuisance: - sigma_NL = data.mcmc_parameters['sigma_NL']['current']*data.mcmc_parameters['sigma_NL']['scale'] + if 'beta_0^SKA1' in self.use_nuisance: + b = self.c4*np.exp(self.c5*self.z_mean*data.mcmc_parameters['beta_1^SKA1']['current']*data.mcmc_parameters['beta_1^SKA1']['scale'])*data.mcmc_parameters['beta_0^SKA1']['current']*data.mcmc_parameters['beta_0^SKA1']['scale'] + else: + b = self.c4*np.exp(self.c5*self.z_mean) + + # Compute freq.res. sigma_r = (1+z)^2/H*sigma_nu/nu_21cm, nu in Mhz + # Compute ang.res. sigma_perp = 1.22(1+z)^2*D_A*lambda_21cm/Baseline, Baseline in km + # combine into exp(-k^2*(mu^2*(sig_r^2-sig_perp^2)+sig_perp^2)) independent of cosmo + # used as exp(-k^2*(mu^2*sigma_A+sigma_B)) all fiducial + if self.fid_values_exist is False: + sigma_A = np.zeros(self.nbin,'float64') + sigma_B = np.zeros(self.nbin,'float64') + sigma_A = ((1.+self.z_mean[:])**2/H[1::2]*self.delta_nu/np.sqrt(8.*np.log(2.))/self.nu0)**2 -( + 1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2 + sigma_B = (1./np.sqrt(8.*np.log(2.))*(1+self.z_mean[:])**2 * D_A[1::2]*2.111e-4/self.Baseline)**2 + + # sigma_NL in Mpc = nonlinear dispersion scale of RSD (1405.1452v2) + sigma_NL = 0.0 # fiducial would be 7 but when kept constant that is more constraining than keeping 0 + if 'sigma_NL' in self.use_nuisance: + sigma_NL = data.mcmc_parameters['sigma_NL']['current']*data.mcmc_parameters['sigma_NL']['scale'] # If the fiducial model does not exist, recover the power spectrum and # store it, then exit. @@ -203,8 +204,8 @@ def loglkl(self, cosmo, data): for index_z in xrange(self.nbin): fid_file.write('%.8g %.8g\n' % (V_survey[index_z],b[index_z])) for index_z in xrange(self.nbin): - fid_file.write('%.8g %.8g\n' % (sigma_A[index_z], sigma_B[index_z])) - fid_file.write('%.8g\n' % sigma_NL) + fid_file.write('%.8g %.8g\n' % (sigma_A[index_z], sigma_B[index_z])) + fid_file.write('%.8g\n' % sigma_NL) print('\n') warnings.warn( "Writing fiducial model in %s, for %s likelihood\n" % ( @@ -220,20 +221,20 @@ def loglkl(self, cosmo, data): # Compute the beta_fid function, for observed spectrum, # beta_fid(k_fid,z) = 1/2b(z) * d log(P_nl_fid(k_fid,z))/d log a # = -1/2b(z)* (1+z) d log(P_nl_fid(k_fid,z))/dz - if self.use_linear_rsd: + if self.use_linear_rsd: beta_fid = -0.5/self.b_fid*(1+self.z_mean)*np.log( - self.pk_lin_fid[:, 2::2]/self.pk_lin_fid[:, :-2:2])/self.dz - else: + self.pk_lin_fid[:, 2::2]/self.pk_lin_fid[:, :-2:2])/self.dz + else: beta_fid = -0.5/self.b_fid*(1+self.z_mean)*np.log( - self.pk_nl_fid[:, 2::2]/self.pk_nl_fid[:, :-2:2])/self.dz + self.pk_nl_fid[:, 2::2]/self.pk_nl_fid[:, :-2:2])/self.dz # Compute the tilde P_fid(k_ref,z,mu) = H_fid(z)/D_A_fid(z)**2 ( 1 + beta_fid(k_fid,z)mu^2)^2 P_nl_fid(k_fid,z)exp(-k_fid^2*(mu_fid^2*sigma_A(z)+sigma_B(z))) self.tilde_P_fid = np.zeros((self.k_size, self.nbin, self.mu_size),'float64') self.tilde_P_fid = self.H_fid[na, 1::2, na]/( - self.D_A_fid[na, 1::2, na])**2*self.b_fid[na,:,na]**2*( + self.D_A_fid[na, 1::2, na])**2*self.b_fid[na,:,na]**2*( 1. + beta_fid[:, :, na] * self.mu_fid[na, na, :]**2)**2 * ( - self.pk_nl_fid[:, 1::2, na]) * np.exp(-self.k_fid[:,na,na]**2 * - (self.mu_fid[na, na, :]**2*(self.sigma_A_fid[na,:,na]+self.sigma_NL_fid**2) + self.sigma_B_fid[na,:,na])) + self.pk_nl_fid[:, 1::2, na]) * np.exp(-self.k_fid[:,na,na]**2 * + (self.mu_fid[na, na, :]**2*(self.sigma_A_fid[na,:,na]+self.sigma_NL_fid**2) + self.sigma_B_fid[na,:,na])) ###################### # TH PART @@ -245,11 +246,11 @@ def loglkl(self, cosmo, data): for index_z in xrange(2*self.nbin+1): self.k[index_k,index_z,:] = np.sqrt((1.-self.mu_fid[:]**2)*self.D_A_fid[index_z]**2/D_A[index_z]**2 + self.mu_fid[:]**2*H[index_z]**2/self.H_fid[index_z]**2 )*self.k_fid[index_k] - # Compute values of mu based on fiducial values: - # mu^2 = mu_fid^2 / (mu_fid^2 + ((H_fid*D_A_fid)/(H*D_A))^2)*(1 - mu_fid^2)) - self.mu = np.zeros((self.nbin,self.mu_size),'float64') - for index_z in xrange(self.nbin): - self.mu[index_z,:] = np.sqrt(self.mu_fid[:]**2/(self.mu_fid[:]**2 + ((self.H_fid[2*index_z+1]*self.D_A_fid[2*index_z+1])/(D_A[2*index_z+1]*H[2*index_z+1]))**2 * (1.-self.mu_fid[:]**2))) + # Compute values of mu based on fiducial values: + # mu^2 = mu_fid^2 / (mu_fid^2 + ((H_fid*D_A_fid)/(H*D_A))^2)*(1 - mu_fid^2)) + self.mu = np.zeros((self.nbin,self.mu_size),'float64') + for index_z in xrange(self.nbin): + self.mu[index_z,:] = np.sqrt(self.mu_fid[:]**2/(self.mu_fid[:]**2 + ((self.H_fid[2*index_z+1]*self.D_A_fid[2*index_z+1])/(D_A[2*index_z+1]*H[2*index_z+1]))**2 * (1.-self.mu_fid[:]**2))) # Recover the non-linear power spectrum from the cosmological module on all # the z_boundaries, to compute afterwards beta. This is pk_nl_th from the @@ -266,31 +267,31 @@ def loglkl(self, cosmo, data): if self.use_linear_rsd: pk_lin_th = cosmo.get_pk_cb_lin(self.k,self.z,self.k_size,2*self.nbin+1,self.mu_size) - if self.UseTheoError : - # Recover the non_linear scale computed by halofit. - #self.k_sigma = np.zeros(2*self.nbin+1, 'float64') - #self.k_sigma = cosmo.nonlinear_scale(self.z,2*self.nbin+1) - - # Define the theoretical error envelope - self.alpha = np.zeros((self.k_size,self.nbin,self.mu_size),'float64') - th_c1 = 0.75056 - th_c2 = 1.5120 - th_a1 = 0.014806 - th_a2 = 0.022047 - for index_z in xrange(self.nbin): - k_z = cosmo.h()*pow(1.+self.z_mean[index_z],2./(2.+cosmo.n_s())) - for index_mu in xrange(self.mu_size): - for index_k in xrange(self.k_size): - if self.k[index_k,2*index_z+1,index_mu]/k_z<0.3: - self.alpha[index_k,index_z,index_mu] = th_a1*np.exp(th_c1*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z)) - else: - self.alpha[index_k,index_z,index_mu] = th_a2*np.exp(th_c2*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z)) - - # Define fractional theoretical error variance R/P^2 - self.R_var = np.zeros((self.k_size,self.nbin,self.mu_size),'float64') - for index_k in xrange(self.k_size): - for index_z in xrange(self.nbin): - self.R_var[index_k,index_z,:] = self.V_fid[index_z]/(2.*np.pi)**2*self.k_CorrLength_hMpc*cosmo.h()/self.z_CorrLength*self.dz*self.k_fid[index_k]**2*self.alpha[index_k,index_z,:]**2 + if self.UseTheoError : + # Recover the non_linear scale computed by halofit. + #self.k_sigma = np.zeros(2*self.nbin+1, 'float64') + #self.k_sigma = cosmo.nonlinear_scale(self.z,2*self.nbin+1) + + # Define the theoretical error envelope + self.alpha = np.zeros((self.k_size,self.nbin,self.mu_size),'float64') + th_c1 = 0.75056 + th_c2 = 1.5120 + th_a1 = 0.014806 + th_a2 = 0.022047 + for index_z in xrange(self.nbin): + k_z = cosmo.h()*pow(1.+self.z_mean[index_z],2./(2.+cosmo.n_s())) + for index_mu in xrange(self.mu_size): + for index_k in xrange(self.k_size): + if self.k[index_k,2*index_z+1,index_mu]/k_z<0.3: + self.alpha[index_k,index_z,index_mu] = th_a1*np.exp(th_c1*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z)) + else: + self.alpha[index_k,index_z,index_mu] = th_a2*np.exp(th_c2*np.log10(self.k[index_k,2*index_z+1,index_mu]/k_z)) + + # Define fractional theoretical error variance R/P^2 + self.R_var = np.zeros((self.k_size,self.nbin,self.mu_size),'float64') + for index_k in xrange(self.k_size): + for index_z in xrange(self.nbin): + self.R_var[index_k,index_z,:] = self.V_fid[index_z]/(2.*np.pi)**2*self.k_CorrLength_hMpc*cosmo.h()/self.z_CorrLength*self.dz*self.k_fid[index_k]**2*self.alpha[index_k,index_z,:]**2 # Compute the beta function for nl, # beta(k,z) = 1/2b(z) * d log(P_nl_th (k,z))/d log a @@ -315,47 +316,47 @@ def loglkl(self, cosmo, data): self.P_shot[index_z] = self.H_fid[2*index_z+1]/(self.D_A_fid[2*index_z+1]**2)*self.V_fid[index_z]/self.N_g[index_z] # finally compute chi2, for each z_mean - if self.use_zscaling==0: - # redshift dependent cutoff makes integration more complicated - chi2 = 0.0 - index_kmax = 0 - delta_mu = self.mu_fid[1] - self.mu_fid[0] # equally spaced - integrand_low = 0.0 - integrand_hi = 0.0 - - for index_z in xrange(self.nbin): - # uncomment printers to get contributions from individual redshift bins - #printer1 = chi2*delta_mu - # uncomment to display max. kmin (used to infer kmin~0.02): - #kmin: #print("z=" + str(self.z_mean[index_z]) + " kmin=" + str(34.56/r[2*index_z+1]) + "\tor " + str(6.283/(r[2*index_z+2]-r[2*index_z]))) - for index_k in xrange(1,self.k_size): - if ((self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[self.k_size-index_k]) > -1.e-6): - index_kmax = self.k_size-index_k - break - integrand_low = self.integrand(0,index_z,0)*.5 - for index_k in xrange(1,index_kmax+1): - integrand_hi = self.integrand(index_k,index_z,0)*.5 - chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1]) - integrand_low = integrand_hi - chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax]) - for index_mu in xrange(1,self.mu_size-1): - integrand_low = self.integrand(0,index_z,index_mu) - for index_k in xrange(1,index_kmax+1): - integrand_hi = self.integrand(index_k,index_z,index_mu) - chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1]) - integrand_low = integrand_hi - chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax]) - integrand_low = self.integrand(0,index_z,self.mu_size-1)*.5 - for index_k in xrange(1,index_kmax+1): - integrand_hi = self.integrand(index_k,index_z,self.mu_size-1)*.5 - chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1]) - integrand_low = integrand_hi - chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax]) - #printer2 = chi2*delta_mu-printer1 - #print("%s\t%s" % (self.z_mean[index_z], printer2)) - chi2 *= delta_mu - - else: + if self.use_zscaling==0: + # redshift dependent cutoff makes integration more complicated + chi2 = 0.0 + index_kmax = 0 + delta_mu = self.mu_fid[1] - self.mu_fid[0] # equally spaced + integrand_low = 0.0 + integrand_hi = 0.0 + + for index_z in xrange(self.nbin): + # uncomment printers to get contributions from individual redshift bins + #printer1 = chi2*delta_mu + # uncomment to display max. kmin (used to infer kmin~0.02): + #kmin: #print("z=" + str(self.z_mean[index_z]) + " kmin=" + str(34.56/r[2*index_z+1]) + "\tor " + str(6.283/(r[2*index_z+2]-r[2*index_z]))) + for index_k in xrange(1,self.k_size): + if ((self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[self.k_size-index_k]) > -1.e-6): + index_kmax = self.k_size-index_k + break + integrand_low = self.integrand(0,index_z,0)*.5 + for index_k in xrange(1,index_kmax+1): + integrand_hi = self.integrand(index_k,index_z,0)*.5 + chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1]) + integrand_low = integrand_hi + chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax]) + for index_mu in xrange(1,self.mu_size-1): + integrand_low = self.integrand(0,index_z,index_mu) + for index_k in xrange(1,index_kmax+1): + integrand_hi = self.integrand(index_k,index_z,index_mu) + chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1]) + integrand_low = integrand_hi + chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax]) + integrand_low = self.integrand(0,index_z,self.mu_size-1)*.5 + for index_k in xrange(1,index_kmax+1): + integrand_hi = self.integrand(index_k,index_z,self.mu_size-1)*.5 + chi2 += (integrand_hi+integrand_low)*.5*(self.k_fid[index_k]-self.k_fid[index_k-1]) + integrand_low = integrand_hi + chi2 += integrand_low*(self.k_cut(self.z_mean[index_z],cosmo.h(),cosmo.n_s())-self.k_fid[index_kmax]) + #printer2 = chi2*delta_mu-printer1 + #print("%s\t%s" % (self.z_mean[index_z], printer2)) + chi2 *= delta_mu + + else: chi2 = 0.0 mu_integrand_lo,mu_integrand_hi = 0.0,0.0 k_integrand = np.zeros(self.k_size,'float64')