Skip to content

Commit

Permalink
version check
Browse files Browse the repository at this point in the history 
for ccl
  • Loading branch information
@eunseongleee
eunseongleee committed Sep 27, 2024
1 parent 080c7a2 commit cd3d4dd
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 32 deletions.
12 changes: 7 additions & 5 deletions soliket/clusters/ccl_th.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ class CCL(Theory):
# CCL options
transfer_function: str = 'boltzmann_camb'
matter_pk: str = 'halofit'
baryons_pk: str = 'nobaryons'
baryons_pk: None # CCL v3
md_hmf: str = '200m'
# Params it can accept
params = {'Omega_c': None,
Expand Down Expand Up @@ -126,6 +126,8 @@ def get_can_support_params(self):

def calculate(self, state, want_derived=True, **params_values_dict):
# Generate the CCL cosmology object which can then be used downstream
if self.baryons_pk == 'nobaryons': # For CCL v3, in case people don't want to update SOLikeT .yml configs
self.baryons_pk=None
cosmo = ccl.Cosmology(Omega_c=self.provider.get_param('Omega_c'),
Omega_b=self.provider.get_param('Omega_b'),
h=self.provider.get_param('h'),
Expand All @@ -135,7 +137,7 @@ def calculate(self, state, want_derived=True, **params_values_dict):
m_nu=self.provider.get_param('m_nu'),
transfer_function=self.transfer_function,
matter_power_spectrum=self.matter_pk,
baryons_power_spectrum=self.baryons_pk)
baryonic_effects=self.baryons_pk) # baryonic_effects is CCL v3

state['derived'] = {'H0': cosmo.cosmo.params.H0}
for req_res, attrs in self._required_results.items():
Expand All @@ -149,14 +151,14 @@ def calculate(self, state, want_derived=True, **params_values_dict):
state[req_res] = None
elif req_res == 'nc_data':
if self.md_hmf == '200m':
md = ccl.halos.MassDef200m(c_m='Bhattacharya13')
md = ccl.halos.MassDef(200, 'matter')
elif self.md_hmf == '200c':
md = ccl.halos.MassDef200c(c_m='Bhattacharya13')
md = ccl.halos.MassDef(200, 'critical')
elif self.md_hmf == '500c':
md = ccl.halos.MassDef(500, 'critical')
else:
raise NotImplementedError('Only md_hmf = 200m, 200c and 500c currently supported.')
mf = ccl.halos.MassFuncTinker08(cosmo, mass_def=md)
mf = ccl.halos.MassFuncTinker08(mass_def=md)
state[req_res] = {'HMF': mf,
'md': md}
elif req_res == 'CCL':
Expand Down
110 changes: 87 additions & 23 deletions soliket/clusters/clusters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,7 @@ def _get_integrated(self, pk_intp, **kwargs):
delN2D = np.zeros((len(zarr), Nq))

# integrate over mass
dndzz = integrate.simpson(intgr[None, :,:]*cc, x=self.lnmarr, axis=2)
dndzz = integrate.simpson(intgr[None, ...]*cc, x=self.lnmarr, axis=2)

# integrate over fine z bins and sum over in larger z bins
for i in range(len(zarr)):
Expand Down Expand Up @@ -314,16 +314,16 @@ def _get_completeness(self, marr, zarr, y0, qbin, marr_500c=None, **params):
arg = get_erf_diff(yy0, noise[i], qmin, qmax, qcut, dof=self.selfunc['debiasDOF'])

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]
arg0 = np.float32((lnyy[:, None, None] - mu[i])/(np.sqrt(2.)*scatter))
args = fac * np.exp(np.float32(-arg0**2.)) * cc[:, None, None]

# truncation for debiasing:
qtab_filter = yy0 / noise[i]
qtab_filter[qtab_filter < self.selfunc['debias_cutoff']] = 0.
qtab_filter[qtab_filter >= self.selfunc['debias_cutoff']] = 1.
# end truncation for debiasing.

comp += integrate.simpson(np.float32(args)*qtab_filter[:,None,None], x=lnyy, axis=0)
comp += integrate.simpson(np.float32(args)*qtab_filter[:, None, None], x=lnyy, axis=0)

comp[comp < 0.] = 0.
comp[comp > 1.] = 1.
Expand Down Expand Up @@ -454,7 +454,7 @@ def _get_n_expected(self, pk_intp, **kwargs):

if self.selfunc['method'] == 'injection':
Pfunc = self.Pfunc_inj(marr, zz, **kwargs)
Pfunc = np.repeat(Pfunc[:,:, np.newaxis], Ynoise.shape[0], axis=2)
Pfunc = np.repeat(Pfunc[..., np.newaxis], Ynoise.shape[0], axis=2)
else:
Pfunc = self.PfuncY(Ynoise, marr, zz, kwargs)

Expand All @@ -465,7 +465,7 @@ def _get_n_expected(self, pk_intp, **kwargs):
zcut_arr = np.arange(zcut)
Ntot = np.sum(np.delete(Nz, zcut_arr, 0))
else:
Nz = integrate.simpson(dndlnm * Pfunc[:,:,index], self.lnmarr[:, None], axis=0)
Nz = integrate.simpson(dndlnm * Pfunc[..., index], self.lnmarr[:, None], axis=0)
Ntot += integrate.simpson(Nz * dVdz, x=zz) * frac

self.log.info("Total predicted N = {}".format(Ntot))
Expand Down Expand Up @@ -531,13 +531,22 @@ def Prob_per_cluster(z, tsz_signal, tsz_signal_err, tile_name):
Pfunc_ind = self.Pfunc_per(tile_index, marr, c_z, c_y, c_yerr, kwargs).T
dn_dlnm = np.squeeze(dndlnm_intp(c_z, self.lnmarr))
dVdz = get_dVdz(self, c_z, dVdz_interp=True)
ans = integrate.simpson(dn_dlnm * Pfunc_ind * dVdz[None,:], x=self.lnmarr[:,None], axis=0)
#ans = integrate.simpson(dn_dlnm * Pfunc_ind * dVdz[None,:] * self.skyfracs.sum(), x=self.lnmarr[:,None], axis=0)

if self.selfunc['bias_handler'] == 'theory':

skyfracs = self.skyfracs / self.skyfracs.sum()
ans = 0
for index, frac in enumerate(skyfracs):
ans += integrate.simpson(dn_dlnm[None, ...] * Pfunc_ind[index,...] * dVdz[None, None, :] * frac, x=self.lnmarr[None, :, None], axis=1)

else:
ans = integrate.simpson(dn_dlnm * Pfunc_ind * dVdz[None, :], x=self.lnmarr[:, None], axis=0)

return ans

return Prob_per_cluster

def P_Yo(self, tile_index, LgY, marr, z, params):
def P_Yo(self, tile_index, LgY, Ynoise, marr, z, params):

if self.theorypred['md_hmf'] != self.theorypred['md_ym']:
marr_ymmd = convert_masses(self, marr, z)
Expand All @@ -550,13 +559,19 @@ def P_Yo(self, tile_index, LgY, marr, z, params):

Ytilde = get_y0(self, marr_ymmd, z, marr_500c, use_Q=True, Ez_interp=True, tile_index=tile_index, **params)

Ytilde = np.repeat(Ytilde[:,:, np.newaxis], LgY.shape[-1], axis=2)
LgY = np.repeat(LgY[np.newaxis, :,:], Ytilde.shape[0], axis=0)
Ytilde = np.repeat(Ytilde[..., np.newaxis], LgY.shape[-1], axis=2)
LgY = np.repeat(LgY[np.newaxis, ...], Ytilde.shape[0], axis=0)

Y = np.exp(LgY)
if self.selfunc['bias_handler'] == 'theory':
Ytilde = np.repeat(Ytilde[..., np.newaxis], Ynoise.shape[0], axis=3)
LgY = np.repeat(LgY[..., np.newaxis], Ynoise.shape[0], axis=3)

numer = -1.0 * (np.log(Y / Ytilde)) ** 2
trueSNR = Ytilde / Ynoise[None, None, None, :]
opt_bias_corr_factor = _opt_bias_func(trueSNR, params['opt_bias_A'], params['opt_bias_B'])
Ytilde = Ytilde * opt_bias_corr_factor

Y = np.exp(LgY)
numer = -1.0 * (np.log(Y / Ytilde)) ** 2
ans = (
1.0 / (params["scatter_sz"] * np.sqrt(2 * np.pi)) *
np.exp(numer / (2.0 * params["scatter_sz"] ** 2))
Expand All @@ -577,8 +592,8 @@ def P_Yo_vec(self, LgY, Ynoise, marr, z, params):
Y = np.exp(LgY).T
Ytilde = get_y0(self, marr_ymmd, z, marr_500c, use_Q=True, Ez_interp=False, **params)

Y = np.repeat(Y[np.newaxis, :,:,:], Ytilde.shape[0], axis=0)
Ytilde = np.repeat(Ytilde[:, np.newaxis, :,:], Y.shape[1], axis=1)
Y = np.repeat(Y[np.newaxis, ...], Ytilde.shape[0], axis=0)
Ytilde = np.repeat(Ytilde[:, np.newaxis, ...], Y.shape[1], axis=1)

if self.selfunc['bias_handler'] == 'theory':
trueSNR = Ytilde / Ynoise[:, None, None, None]
Expand Down Expand Up @@ -611,7 +626,7 @@ def P_of_gt_SN(self, LgY, marr, z, Ynoise, params):
qcut = np.outer(np.ones(np.shape(Ytilde)), self.qcut)
qcut_a = np.reshape(qcut, (Ytilde.shape[0], Ytilde.shape[1], Ytilde.shape[2]))

Ynoise = np.outer(Ynoise, np.ones(np.shape(Ytilde[0,:,:])))
Ynoise = np.outer(Ynoise, np.ones(np.shape(Ytilde[0, ...])))
Ynoise_a = np.reshape(Ynoise, (Ytilde.shape[0], Ytilde.shape[1], Ytilde.shape[2]))

if self.selfunc['bias_handler'] == 'theory':
Expand Down Expand Up @@ -649,7 +664,7 @@ def P_of_gt_SN(self, LgY, marr, z, Ynoise, params):
qtab_filter[qtab_filter >= self.selfunc['debias_cutoff']] = 1.
# end truncation for debiasing.

ans = integrate.simpson(P_Y * sig_thresh * qtab_filter[None,None,:,:], x=LgY, axis=2) #* np.log(10) # why log10?
ans = integrate.simpson(P_Y * sig_thresh * qtab_filter[None, None, ...], x=LgY, axis=2) #* np.log(10) # why log10?

return ans

Expand All @@ -665,6 +680,9 @@ def Y_prob(self, Y_c, LgY, Ynoise):
return ans

def Pfunc_per(self, tile_index, marr, z, Y_c, Y_err, params):

Ynoise = self.noise

if params["scatter_sz"] == 0:
if self.theorypred['md_hmf'] != self.theorypred['md_ym']:
marr_ymmd = convert_masses(self, marr, z)
Expand All @@ -680,8 +698,21 @@ def Pfunc_per(self, tile_index, marr, z, Y_c, Y_err, params):
# P_Y_sig = np.nan_to_num(self.Y_prob(Y_c[:, None], LgYtilde, Y_err[:, None]))
# ans = P_Y_sig

ans = np.nan_to_num(get_erf(Ytilde, Y_err[:, None], self.qcut)) # dN/dz
if self.selfunc['bias_handler'] == 'theory':

Ytilde = np.repeat(Ytilde[..., np.newaxis], Ynoise.shape[-1], axis=2)

Ynoise = np.outer(Ynoise, np.ones(np.shape(Ytilde[..., 0])))
Ynoise_a = np.reshape(Ynoise, (Ytilde.shape[0], Ytilde.shape[1], Ytilde.shape[2]))

trueSNR = Ytilde / Ynoise_a
opt_bias_corr_factor = _opt_bias_func(trueSNR, params['opt_bias_A'], params['opt_bias_B'])
Ytilde = Ytilde * opt_bias_corr_factor

ans = np.nan_to_num(get_erf(Ytilde, Y_err[:, None, None], self.qcut))

else:
ans = np.nan_to_num(get_erf(Ytilde, Y_err[:, None], self.qcut))
else:
LgY = self.lny
Y = np.exp(LgY)
Expand All @@ -690,10 +721,22 @@ def Pfunc_per(self, tile_index, marr, z, Y_c, Y_err, params):
P_Y_sig = np.nan_to_num(get_erf(Y, Y_err[:, None], self.qcut))
#P_Y_sig = self.Y_prob(Y_c[:, None], LgY, Y_err[:, None])

P_Y = np.nan_to_num(self.P_Yo(tile_index, LgYa, marr, z, params))
P_Y_sig = np.repeat(P_Y_sig[:, np.newaxis, :], P_Y.shape[1], axis=1)
LgY = LgY[None, None, :]
ans = integrate.simpson(P_Y * P_Y_sig, LgY, np.diff(LgY), axis=2)
if self.selfunc['bias_handler'] == 'theory':

P_Y = np.nan_to_num(self.P_Yo(tile_index, LgYa, Ynoise, marr, z, params))

P_Y_sig = np.repeat(P_Y_sig[:, np.newaxis, :], P_Y.shape[1], axis=1)
LgY = LgY[None, None, :]

P_Y_sig = np.repeat(P_Y_sig[..., np.newaxis], Ynoise.shape[0], axis=3)
LgY = np.repeat(LgY[..., np.newaxis], Ynoise.shape[0], axis=3)
ans = integrate.simpson(P_Y * P_Y_sig, LgY, np.diff(LgY), axis=2)

else:
P_Y = np.nan_to_num(self.P_Yo(tile_index, LgYa, Ynoise, marr, z, params))
P_Y_sig = np.repeat(P_Y_sig[:, np.newaxis, :], P_Y.shape[1], axis=1)
LgY = LgY[None, None, :]
ans = integrate.simpson(P_Y * P_Y_sig, LgY, np.diff(LgY), axis=2)

return ans

Expand Down Expand Up @@ -761,6 +804,27 @@ def initialize_common(self):
cat_tsz_signal_err = cat_tab['fixed_err_y_c'].data.astype(float)
cat_tile_name = np.array(cat_tab['tileName'].data, dtype = str)

# #----------------------- zcut
keep = zcat < self.binning['z']['zmax']
qcat = qcat[keep]
cat_tsz_signal = cat_tsz_signal[keep]
cat_tsz_signal_err = cat_tsz_signal_err[keep]
cat_tile_name = cat_tile_name[keep]
zcat = zcat[keep]
# #----------------------- zcut


# #----------------------- split
# if self.selfunc['split'] == 'east':
# RAcat = cat_tab['RADeg'].data.astype(float)
# keep = RAcat/15. > 12.
# zcat = zcat[keep]
# qcat = qcat[keep]
# cat_tsz_signal = cat_tsz_signal[keep]
# cat_tsz_signal_err = cat_tsz_signal_err[keep]
# cat_tile_name = cat_tile_name[keep]
# #----------------------- split

# Optimization bias handler
if self.selfunc['bias_handler'] not in ['theory', 'catalog']:
raise NotImplementedError('bias_handler should be either "theory" or "catalog"')
Expand Down Expand Up @@ -1458,7 +1522,7 @@ def rel(m):
else:
splQ = 1.

y0 = A0 * (Ez ** 2.) * (mb / Mpivot) ** (1. + B0) * splQ
y0 = A0 * (Ez ** C0) * (mb / Mpivot) ** (1. + B0) * splQ
y0[y0 <= 0] = 1e-9

return y0
Expand Down
5 changes: 2 additions & 3 deletions soliket/clusters/nemo_mocks.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,10 +93,9 @@ def make_nemo_mock(configdict):
def get_nemo_pred(configdict, zbins):

selFn=completeness.SelFn(configdict['path2selFn'], SNRCut = configdict['predSNRCut'], zStep = configdict['selFnZStep'],
enableDrawSample = configdict['makeMock'], massFunction = configdict['massFunc'],
applyRelativisticCorrection = configdict['relativisticCorrection'],
massFunction = configdict['massFunc'], applyRelativisticCorrection = configdict['relativisticCorrection'],
rhoType = configdict['rhoType'], delta = configdict['delta'], method=configdict['method'],
QSource=configdict['QSource'])
QSource=configdict['QSource']) #enableDrawSample = configdict['makeMock'],

predMz=selFn.compMz*selFn.mockSurvey.clusterCount
#TODO: Ask Matt where the minimal mass comes from
Expand Down
2 changes: 1 addition & 1 deletion soliket/poisson.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ def _get_n_expected(self, **kwargs):
raise NotImplementedError

def logp(self, **kwargs):
pk_intp = self.theory.get_Pk_interpolator()
pk_intp = self.provider.get_Pk_interpolator()
rate_densities = self._get_rate_fn(pk_intp, **kwargs)
n_expected = self._get_n_expected(pk_intp, **kwargs)
return self.data.loglike(rate_densities, n_expected)
8 changes: 8 additions & 0 deletions version.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
# Note that we need to fall back to the hard-coded version if either
# setuptools_scm can't be imported or setuptools_scm can't determine the
# version, so we catch the generic 'Exception'.
try:
from setuptools_scm import get_version
version = get_version(root='..', relative_to=__file__)
except Exception:
version = '0.1rc1.dev101+g7bc3f87'

0 comments on commit cd3d4dd

Please sign in to comment.