Skip to content

Commit

Permalink
Merge 598f3d7 into d5408a7
Browse files Browse the repository at this point in the history
  • Loading branch information
@damonge
damonge committed Jul 6, 2022
2 parents d5408a7 + 598f3d7 commit 0a63c44
Show file tree
Hide file tree
Showing 8 changed files with 108 additions and 51 deletions.
14 changes: 7 additions & 7 deletions benchmarks/ccl_test_f2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ CTEST2(f2d,constant) {
2,
2,
ccl_f2d_constantgrowth,
1,0,2,
1,0,0,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand All @@ -108,7 +108,7 @@ CTEST2(f2d,constant) {
2,
2,
ccl_f2d_constantgrowth,
1,0,2,
1,0,0,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand All @@ -127,7 +127,7 @@ CTEST2(f2d,constant) {
2,
2,
ccl_f2d_constantgrowth,
1,0,2,
1,0,0,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand Down Expand Up @@ -156,7 +156,7 @@ CTEST2(f2d,a_overflow) {
2, //extrap_hik
ccl_f2d_constantgrowth, //extrap_growth
1, //is_fka_log
0,2,
0,0,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand Down Expand Up @@ -186,7 +186,7 @@ CTEST2(f2d,sanity) {
2, //extrap_hik
ccl_f2d_constantgrowth, //extrap_growth
1, //is_fka_log
1,2,
0,1,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand Down Expand Up @@ -243,7 +243,7 @@ CTEST2(f2d,factorize) {
2, //extrap_hik
ccl_f2d_constantgrowth, //extrap_growth
1, //is_fka_log
1,2,
0,1,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand Down Expand Up @@ -315,7 +315,7 @@ CTEST2(f2d,pk) {
2, //extrap_hik
ccl_f2d_cclgrowth, //extrap_growth
1, //is_fka_log
0,2,
0,0,2,
ccl_f2d_3,
&status);
ASSERT_TRUE(status==0);
Expand Down
3 changes: 3 additions & 0 deletions include/ccl_f2d.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ typedef struct {
int extrap_order_hik; /**< Order of extrapolating polynomial in log(k) for high k (0, 1 or 2)*/
ccl_f2d_extrap_growth_t extrap_linear_growth; /**< Extrapolation type at high redshifts*/
int is_log; /**< Do I hold the values of log(f(k,a))?*/
int extrap_in_log; /**< Do I hold the values of log(f(k,a))?*/
double growth_factor_0; /**< Constant extrapolating growth factor*/
int growth_exponent; /**< Power to which growth should be exponentiated*/
gsl_spline *fk; /**< Spline holding the values of the k-dependent factor*/
Expand All @@ -56,6 +57,7 @@ typedef struct {
* @param extrap_order_lok Order of the polynomial that extrapolates on wavenumbers smaller than the minimum of lk_arr. Allowed values: 0 (constant), 1 (linear extrapolation) and 2 (quadratic extrapolation). Extrapolation happens in ln(k).
* @param extrap_order_hik Order of the polynomial that extrapolates on wavenumbers larger than the maximum of lk_arr. Allowed values: 0 (constant), 1 (linear extrapolation) and 2 (quadratic extrapolation). Extrapolation happens in ln(k).
* @param extrap_linear_growth: ccl_f2d_extrap_growth_t value defining how the function with scale factors below the interpolation range. Allowed values: ccl_f2d_cclgrowth (scale with the CCL linear growth factor), ccl_f2d_constantgrowth (scale by multiplying the function at the earliest available scale factor by a constant number, defined by `growth_factor_0`), ccl_f2d_no_extrapol (throw an error if the function is ever evaluated outside the interpolation range in a). Note that, above the interpolation range (i.e. for low redshifts), the function will be assumed constant.
* @param extrap_in_log: if not zero, extrapolation will take place in log-log scale.
* @param is_fka_log: if not zero, `fka_arr` contains ln(f(k,a)) instead of f(k,a). If the function is factorizable, then `fk_arr` holds ln(K(k)) and `fa_arr` holds ln(A(a)), where f(k,a)=K(k)*A(a).
* @param growth_factor_0: custom growth function. Irrelevant if extrap_linear_growth!=ccl_f2d_constantgrowth.
* @param growth_exponent: power to which the extrapolating growth factor should be exponentiated when extrapolating (e.g. usually 2 for linear power spectra).
Expand All @@ -72,6 +74,7 @@ ccl_f2d_t *ccl_f2d_t_new(int na,double *a_arr,
int extrap_order_hik,
ccl_f2d_extrap_growth_t extrap_linear_growth,
int is_fka_log,
int extrap_in_log,
double growth_factor_0,
int growth_exponent,
ccl_f2d_interp_t interp_type,
Expand Down
5 changes: 3 additions & 2 deletions pyccl/ccl_pk2d.i
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,13 @@ ccl_f2d_t *set_pk2d_new_from_arrays(double* lkarr,int nk,
double* aarr,int na,
double* pkarr,int npk,
int order_lok,int order_hik,
int is_logp,
int is_logp, int extrap_in_log,
int *status)
{
ccl_f2d_t *psp=ccl_f2d_t_new(na,aarr,nk,lkarr,pkarr,NULL,NULL,0,
order_lok,order_hik,ccl_f2d_cclgrowth,
is_logp,0,2,ccl_f2d_3,status);
is_logp,extrap_in_log,0,2,
ccl_f2d_3,status);
return psp;
}

Expand Down
12 changes: 8 additions & 4 deletions pyccl/pk2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,11 +47,11 @@ class Pk2D(object):
extrap_order_lok (int): extrapolation order to be used on k-values
below the minimum of the splines (use 0, 1 or 2). Note that
the extrapolation will be done in either log(P(k)) or P(k),
depending on the value of `is_logp`.
depending on the value of `is_logp` and `extrap_in_log`.
extrap_order_hik (int): extrapolation order to be used on k-values
above the maximum of the splines (use 0, 1 or 2). Note that
the extrapolation will be done in either log(P(k)) or P(k),
depending on the value of `is_logp`.
depending on the value of `is_logp` and `extrap_in_log`.
is_logp (boolean): if True, pkfunc/pkarr return/hold the natural
logarithm of the power spectrum. Otherwise, the true value
of the power spectrum is expected. Note that arrays will be
Expand All @@ -62,10 +62,12 @@ class Pk2D(object):
wavenumber.
empty (bool): if True, just create an empty object, to be filled
out later
extrap_in_log (bool): if True, extrapolation will be done in
log(k)-log(p) space, regardless of the value of `is_logp`.
"""
def __init__(self, pkfunc=None, a_arr=None, lk_arr=None, pk_arr=None,
is_logp=True, extrap_order_lok=1, extrap_order_hik=2,
cosmo=None, empty=False):
cosmo=None, empty=False, extrap_in_log=False):
if empty:
self.has_psp = False
return
Expand Down Expand Up @@ -117,7 +119,9 @@ def __init__(self, pkfunc=None, a_arr=None, lk_arr=None, pk_arr=None,
self.psp, status = lib.set_pk2d_new_from_arrays(lk_arr, a_arr, pkflat,
int(extrap_order_lok),
int(extrap_order_hik),
int(is_logp), status)
int(is_logp),
int(extrap_in_log),
status)
check(status, cosmo=cosmo)
self.has_psp = True

Expand Down
112 changes: 80 additions & 32 deletions src/ccl_f2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ ccl_f2d_t *ccl_f2d_t_copy(ccl_f2d_t *f2d_o, int *status)
f2d->extrap_linear_growth = f2d_o->extrap_linear_growth;
f2d->extrap_order_lok = f2d_o->extrap_order_lok;
f2d->extrap_order_hik = f2d_o->extrap_order_hik;
f2d->extrap_in_log = f2d_o->extrap_in_log;
f2d->is_log = f2d_o->is_log;
f2d->growth_factor_0 = f2d_o->growth_factor_0;
f2d->growth_exponent = f2d_o->growth_exponent;
Expand Down Expand Up @@ -99,6 +100,7 @@ ccl_f2d_t *ccl_f2d_t_new(int na,double *a_arr,
int extrap_order_hik,
ccl_f2d_extrap_growth_t extrap_linear_growth,
int is_fka_log,
int extrap_in_log,
double growth_factor_0,
int growth_exponent,
ccl_f2d_interp_t interp_type,
Expand All @@ -115,6 +117,7 @@ ccl_f2d_t *ccl_f2d_t_new(int na,double *a_arr,
f2d->is_a_constant = ((a_arr == NULL) || ((fka_arr == NULL) && (fa_arr == NULL)));
f2d->extrap_order_lok = extrap_order_lok;
f2d->extrap_order_hik = extrap_order_hik;
f2d->extrap_in_log = extrap_in_log;
f2d->extrap_linear_growth = extrap_linear_growth;
f2d->is_log = is_fka_log;
f2d->growth_factor_0 = growth_factor_0;
Expand Down Expand Up @@ -285,11 +288,17 @@ double ccl_f2d_t_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, int *status
}

// Now extrapolate in k if needed
// Stored in lin-scale but extrapolated in log
int extrap_linlog = ((!(f2d->is_log)) && (f2d->extrap_in_log));
if (is_hik) {
fka_post = fka_pre;
if (f2d->extrap_order_hik > 0) {
double pd;
double pd, pdd;
double dlk = lk-lk_ev;
if ((extrap_linlog) && (fka_pre==0)) {// Avoid division by zero
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
if (f2d->is_factorizable)
spstatus = gsl_spline_eval_deriv_e(f2d->fk, lk_ev, NULL, &pd);
else
Expand All @@ -298,25 +307,43 @@ double ccl_f2d_t_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, int *status
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
fka_post += pd*dlk;
if (extrap_linlog) {
// dlogP/dlogk = P'/P
pd = pd/fka_pre;
// P(k) = P(k_0) * e^{delta_log_k * dlogP/dlogk}
fka_post *= exp(pd*dlk);
}
else
fka_post += pd*dlk;
if (f2d->extrap_order_hik > 1) {
if (f2d->is_factorizable)
spstatus = gsl_spline_eval_deriv2_e(f2d->fk, lk_ev, NULL, &pd);
spstatus = gsl_spline_eval_deriv2_e(f2d->fk, lk_ev, NULL, &pdd);
else
spstatus = gsl_spline2d_eval_deriv_xx_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pd);
spstatus = gsl_spline2d_eval_deriv_xx_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pdd);
if (spstatus) {
*status=CCL_ERROR_SPLINE_EV;
return NAN;
}
fka_post += pd*dlk*dlk*0.5;
if (extrap_linlog) {
// d^2logP/dlogk^2 = P''/P - (P'/P)^2
pdd = pdd/fka_pre-pd*pd;
// Second order correction
fka_post *= exp(0.5*pdd*dlk*dlk);
}
else
fka_post += pdd*dlk*dlk*0.5;
}
}
}
else if (is_lok) {
fka_post = fka_pre;
if (f2d->extrap_order_lok > 0) {
double pd;
double pd,pdd;
double dlk = lk-lk_ev;
if ((extrap_linlog) && (fka_pre==0)) {// Avoid division by zero
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
if (f2d->is_factorizable)
spstatus = gsl_spline_eval_deriv_e(f2d->fk, lk_ev, NULL, &pd);
else
Expand All @@ -325,18 +352,32 @@ double ccl_f2d_t_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, int *status
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
fka_post += pd*dlk;
if (extrap_linlog) {
// dlogP/dlogk = P'/P
pd = pd/fka_pre;
// P(k) = P(k_0) * e^{delta_log_k * dlogP/dlogk}
fka_post *= exp(pd*dlk);
}
else
fka_post += pd*dlk;

if (f2d->extrap_order_lok > 1) {
if (f2d->is_factorizable)
spstatus = gsl_spline_eval_deriv2_e(f2d->fk, lk_ev, NULL, &pd);
spstatus = gsl_spline_eval_deriv2_e(f2d->fk, lk_ev, NULL, &pdd);
else
spstatus = gsl_spline2d_eval_deriv_xx_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pd);
spstatus = gsl_spline2d_eval_deriv_xx_e(f2d->fka, lk_ev, a_ev, NULL, NULL, &pdd);
if (spstatus) {
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
fka_post += pd*dlk*dlk*0.5;
if (extrap_linlog) {
// d^2logP/dlogk^2 = P''/P - (P'/P)^2
pdd = pdd/fka_pre-pd*pd;
// Second order correction
fka_post *= exp(0.5*pdd*dlk*dlk);
}
else
fka_post += pdd*dlk*dlk*0.5;
}
}
}
Expand Down Expand Up @@ -465,6 +506,17 @@ double ccl_f2d_t_dlogf_dlk_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, i
}

// Now extrapolate in k if needed
//Stored in lin-scale but extrapolated in log
int extrap_linlog = ((!(f2d->is_log)) && (f2d->extrap_in_log));
double fka_0;
if ((is_hik || is_lok) && extrap_linlog) {
// Need to evaluate this to calculate the numerical derivative
fka_0 = ccl_f2d_t_eval(f2d, lk_ev, a_ev, cosmo, status);
if (status || (fka_0==0)) {
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
}
if (is_hik) {
fka_post = fka_pre;
if (f2d->extrap_order_hik > 1) {
Expand All @@ -478,7 +530,15 @@ double ccl_f2d_t_dlogf_dlk_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, i
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
fka_post += pd*dlk;
if (extrap_linlog) {
// fka_0 = P_0
// fka_pre = P'_0
// pd = P''_0
// d^2logP/dlogk^2 = P''/P-(P'/P)^2
fka_post = (fka_pre+(pd/fka_0-fka_pre*fka_pre/(fka_0*fka_0))*dlk)/inv_pk0;
}
else
fka_post += pd*dlk;
}
}
else if (is_lok) {
Expand All @@ -494,7 +554,15 @@ double ccl_f2d_t_dlogf_dlk_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, i
*status = CCL_ERROR_SPLINE_EV;
return NAN;
}
fka_post += pd*dlk;
if (extrap_linlog) {
// fka_0 = P_0
// fka_pre = P'_0
// pd = P''_0
// d^2logP/dlogk^2 = P''/P-(P'/P)^2
fka_post = (fka_pre+(pd/fka_0-fka_pre*fka_pre/(fka_0*fka_0))*dlk)/inv_pk0;
}
else
fka_post += pd*dlk;
}
}
else
Expand All @@ -506,26 +574,6 @@ double ccl_f2d_t_dlogf_dlk_eval(ccl_f2d_t *f2d,double lk,double a,void *cosmo, i

// Extrapolation in a in log-space not needed
// (does not contribute to logarithmic k derivative)
if (is_hiz) {
double gz; // Use CCL's growth function
if (f2d->extrap_linear_growth == ccl_f2d_cclgrowth) {
ccl_cosmology *csm = (ccl_cosmology *)cosmo;
if (!csm->computed_growth) {
*status = CCL_ERROR_GROWTH_INIT;
ccl_cosmology_set_status_message(
csm,
"ccl_f2d.c: ccl_f2d_t_eval(): growth factor splines have not been precomputed!");
return NAN;
}
gz = (
ccl_growth_factor(csm, a, status) /
ccl_growth_factor(csm, a_ev, status));
}
else // Use constant growth factor
gz = f2d->growth_factor_0;

fka_post *= pow(gz, f2d->growth_exponent);
}
}

return fka_post;
Expand Down
4 changes: 2 additions & 2 deletions src/ccl_f3d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,15 +226,15 @@ ccl_f3d_t *ccl_f3d_t_new(int na,double *a_arr,
fka1_arr, NULL, NULL, 0,
extrap_order_lok, extrap_order_hik,
extrap_linear_growth,
is_tkka_log,
is_tkka_log, 0,
growth_factor_0, growth_exponent/2,
interp_type, status);
f3d->fka_2 = ccl_f2d_t_new(na, a_arr,
nk, lk_arr,
fka2_arr, NULL, NULL, 0,
extrap_order_lok, extrap_order_hik,
extrap_linear_growth,
is_tkka_log,
is_tkka_log, 0,
growth_factor_0, growth_exponent/2,
interp_type, status);
}
Expand Down

0 comments on commit 0a63c44

Please sign in to comment.