diff --git a/soliket/clusters/clusters.py b/soliket/clusters/clusters.py index 99f8bbad..fa162b5a 100644 --- a/soliket/clusters/clusters.py +++ b/soliket/clusters/clusters.py @@ -233,10 +233,10 @@ def _get_completeness(self, marr, zarr, y0, qbin, marr_500c=None, **params): for i in range(Npatches): if compl_mode == 'erf_prod': - arg.append(get_erf_prod(y0[i], noise[i], qmin, qmax, qcut, kk, Nq)) + arg.append(get_erf_prod(y0[i], noise[i], qmin, qmax, qcut, kk, Nq, dof = self.debiasDOF)) #arg.append(get_stf_prod(y0[i], noise[i], qmin, qmax, qcut, kk, Nq)) elif compl_mode == 'erf_diff': - arg.append(get_erf_diff(y0[i], noise[i], qmin, qmax, qcut)) + arg.append(get_erf_diff(y0[i], noise[i], qmin, qmax, qcut, dof = self.debiasDOF)) #arg.append(get_stf_diff(y0[i], noise[i], qmin, qmax, qcut)) comp = np.einsum('ijk,i->jk', np.nan_to_num(arg), skyfracs) @@ -251,16 +251,23 @@ def _get_completeness(self, marr, zarr, y0, qbin, marr_500c=None, **params): comp = 0. for i in range(Npatches): if compl_mode == 'erf_prod': - arg = get_erf_prod(yy0, noise[i], qmin, qmax, qcut, kk, Nq) + arg = get_erf_prod(yy0, noise[i], qmin, qmax, qcut, kk, Nq, dof = self.debiasDOF) #arg = get_stf_prod(yy0, noise[i], qmin, qmax, qcut, kk, Nq) elif compl_mode == 'erf_diff': - arg = get_erf_diff(yy0, noise[i], qmin, qmax, qcut) + arg = get_erf_diff(yy0, noise[i], qmin, qmax, qcut, dof = self.debiasDOF) #arg = get_stf_diff(yy0, noise[i], qmin, qmax, qcut) cc = np.float32(arg * skyfracs[i]) arg0 = np.float32((lnyy[:, None,None] - mu[i])/(np.sqrt(2.)*scatter)) args = fac * np.exp(np.float32(-arg0**2.)) * cc[:, None,None] - comp += np.trapz(np.float32(args), x=lnyy, axis=0) + + # truncation for debiasing: + qtab_filter = yy0/noise[i] + qtab_filter[qtab_filter=self.debias_cutoff] = 1. + # end truncation for debiasing. + + comp += np.trapz(np.float32(args)*qtab_filter[:,None,None], x=lnyy, axis=0) comp[comp < 0.] = 0. comp[comp > 1.] = 1. @@ -675,8 +682,9 @@ def initialize_common(self): # to print all columns: print(catf[1].columns) # SPT-style SNR bias correction, this should be applied before applying SNR cut - debiasDOF = self.selfunc['debiasDOF'] - qcat = np.sqrt(np.power(qcat, 2) - debiasDOF) + self.debiasDOF = self.selfunc['debiasDOF'] + self.debias_cutoff = self.selfunc['debias_cutoff'] + #qcat = np.sqrt(np.power(qcat, 2) - debiasDOF) # dont debias catalog q! # only above given SNR cut ind = np.where(qcat >= self.qcut)[0] @@ -1195,57 +1203,57 @@ def tinker(sgm, z): # return self.get_dndlnM_at_z_and_M(z,marr) -def get_erf(y, noise, cut): +def get_erf(y, noise, cut, dof = 0): #arg = (y - cut*noise)/np.sqrt(2.)/noise - arg = (y/noise - cut)/np.sqrt(2.) + arg = (np.sqrt((y/noise)**2.+dof) - cut)/np.sqrt(2.) erfc = (special.erf(arg) + 1.)/2. return erfc -def get_stf(y, noise, qcut): - ans = y * 0.0 - ans[y - qcut*noise > 0] = 1.0 - return ans +# def get_stf(y, noise, qcut): +# ans = y * 0.0 +# ans[y - qcut*noise > 0] = 1.0 +# return ans -def get_erf_diff(y, noise, qmin, qmax, qcut): - arg1 = (y/noise - qmax)/np.sqrt(2.) +def get_erf_diff(y, noise, qmin, qmax, qcut, dof = 0): + arg1 = (np.sqrt((y/noise)**2.+dof) - qmax)/np.sqrt(2.) if qmin > qcut: qlim = qmin else: qlim = qcut - arg2 = (y/noise - qlim)/np.sqrt(2.) + arg2 = (np.sqrt((y/noise)**2.+dof) - qlim)/np.sqrt(2.) erf_compl = (special.erf(arg2) - special.erf(arg1)) / 2. return erf_compl -def get_stf_diff(y, noise, qmin, qmax, qcut): - if qmin > qcut: - qmin = qmin - else: - qmin = qcut - ans = y * 0.0 - ans[(y - qmin*noise > 0) & (y - qmax*noise < 0)] = 1.0 - return ans - -def get_erf_prod(y, noise, qmin, qmax, qcut, k, Nq): - - arg0 = get_erf(y, noise, qcut) - arg1 = get_erf(y, noise, qmin) - arg2 = 1. - get_erf(y, noise, qmax) - - if k == 0: arg1 = 1 - if k == Nq-1: arg2 = 1 - - return arg0 * arg1 * arg2 +# def get_stf_diff(y, noise, qmin, qmax, qcut): +# if qmin > qcut: +# qmin = qmin +# else: +# qmin = qcut +# ans = y * 0.0 +# ans[(y - qmin*noise > 0) & (y - qmax*noise < 0)] = 1.0 +# return ans -def get_stf_prod(y, noise, qmin, qmax, qcut, k, Nq): +def get_erf_prod(y, noise, qmin, qmax, qcut, k, Nq, dof = 0): - arg0 = get_stf(y, noise, qcut) - arg1 = get_stf(y, noise, qmin) - arg2 = 1. - get_stf(y, noise, qmax) + arg0 = get_erf(y, noise, qcut, dof) + arg1 = get_erf(y, noise, qmin, dof) + arg2 = 1. - get_erf(y, noise, qmax, dof) if k == 0: arg1 = 1 if k == Nq-1: arg2 = 1 return arg0 * arg1 * arg2 +# +# def get_stf_prod(y, noise, qmin, qmax, qcut, k, Nq): +# +# arg0 = get_stf(y, noise, qcut) +# arg1 = get_stf(y, noise, qmin) +# arg2 = 1. - get_stf(y, noise, qmax) +# +# if k == 0: arg1 = 1 +# if k == Nq-1: arg2 = 1 +# +# return arg0 * arg1 * arg2 def gaussian(xx, mu, sig, noNorm=False): if noNorm: diff --git a/soliket/clusters/input_files/test_binned_lkl_ccl_evaluate.yaml b/soliket/clusters/input_files/test_binned_lkl_ccl_evaluate.yaml index 7c7b1ae2..8da45801 100644 --- a/soliket/clusters/input_files/test_binned_lkl_ccl_evaluate.yaml +++ b/soliket/clusters/input_files/test_binned_lkl_ccl_evaluate.yaml @@ -55,7 +55,8 @@ likelihood: dwnsmpl_bins : 30 # If resolution=downsample, number of bins to use save_dwsmpld : False # Save downsampled Q and rms to npz file and once it exists read those - debiasDOF : 0. + debiasDOF : 3. + debias_cutoff : 2. binning: # redshift bins for number counts z: diff --git a/soliket/clusters/input_files/test_binned_lkl_classy_sz_evaluate.yaml b/soliket/clusters/input_files/test_binned_lkl_classy_sz_evaluate.yaml index fb15736d..c03dc842 100644 --- a/soliket/clusters/input_files/test_binned_lkl_classy_sz_evaluate.yaml +++ b/soliket/clusters/input_files/test_binned_lkl_classy_sz_evaluate.yaml @@ -57,6 +57,7 @@ likelihood: save_dwsmpld : False # Save downsampled Q and rms to npz file and once it exists read those debiasDOF : 0. + debias_cutoff: 0. binning: # redshift bins for number counts z: @@ -191,8 +192,8 @@ theory: use_planck_binned_proba : 0 #use diff of erfs use_skyaveraged_noise: 0 # this will speed-up everything - szcc_dof : 0. - szcc_qtrunc : 0. + szcc_dof : 3. + szcc_qtrunc : 2. stop_at_error: True