From a4e6c1c66df19fd82e06f7e37fd214a4fe9a5334 Mon Sep 17 00:00:00 2001 From: Robert Reischke Date: Mon, 11 May 2026 13:17:12 +0200 Subject: [PATCH 1/4] fixed some potential memory problems --- pyccl/halos/pk_2pt.py | 1 - pyccl/halos/profiles/profile_base.py | 1 - src/ccl_cls.c | 461 +++++++++++++++------------ src/ccl_fftlog.c | 30 +- src/ccl_power.c | 8 +- 5 files changed, 279 insertions(+), 222 deletions(-) diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py index 6a3ce9af9..5260584f7 100644 --- a/pyccl/halos/pk_2pt.py +++ b/pyccl/halos/pk_2pt.py @@ -102,7 +102,6 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *, for ia, aa in enumerate(a_use): # normalizations norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) - if prof2 == prof: norm2 = norm1 else: diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py index 2c9cdbc84..bf56fc6a9 100644 --- a/pyccl/halos/profiles/profile_base.py +++ b/pyccl/halos/profiles/profile_base.py @@ -516,6 +516,5 @@ def get_normalization(self, cosmo, a, *, hmc=None): class HaloProfilePressure(HaloProfile): """Base for pressure halo profiles.""" - class HaloProfileCIB(HaloProfile): """Base for CIB halo profiles.""" diff --git a/src/ccl_cls.c b/src/ccl_cls.c index 0b886a229..35fec16f2 100644 --- a/src/ccl_cls.c +++ b/src/ccl_cls.c @@ -8,7 +8,8 @@ #include "ccl.h" -typedef struct{ +typedef struct +{ double l; ccl_cosmology *cosmo; ccl_cl_tracer_collection_t *trc1; @@ -17,8 +18,8 @@ typedef struct{ int *status; } integ_cl_par; - -typedef struct{ +typedef struct +{ int chipow; double l1; double l2; @@ -38,25 +39,28 @@ static void update_chi_limits(ccl_cl_tracer_collection_t *trc, int is_union) { int itr; - double chimin_h=1E15; - double chimax_h=-1E15; - for(itr=0; itr < trc->n_tracers; itr++) { + double chimin_h = 1E15; + double chimax_h = -1E15; + for (itr = 0; itr < trc->n_tracers; itr++) + { if (trc->ts[itr]->chi_min < chimin_h) chimin_h = trc->ts[itr]->chi_min; if (trc->ts[itr]->chi_max > chimax_h) chimax_h = trc->ts[itr]->chi_max; } - if(is_union) { - if(chimin_h < *chimin) + if (is_union) + { + if (chimin_h < *chimin) *chimin = chimin_h; - if(chimax_h > *chimax) + if (chimax_h > *chimax) *chimax = chimax_h; } - else { - if(chimin_h > *chimin) + else + { + if (chimin_h > *chimin) *chimin = chimin_h; - if(chimax_h < *chimax) + if (chimax_h < *chimax) *chimax = chimax_h; } } @@ -64,13 +68,15 @@ static void update_chi_limits(ccl_cl_tracer_collection_t *trc, static void get_k_interval(ccl_cosmology *cosmo, ccl_cl_tracer_collection_t *trc1, ccl_cl_tracer_collection_t *trc2, - double l, double *lkmin, double *lkmax) { + double l, double *lkmin, double *lkmax) +{ int itr; // Loop through all tracers and find distance bounds double chi_min1 = 1E15; double chi_max1 = -1E15; - for (itr=0; itr < trc1->n_tracers; itr++) { + for (itr = 0; itr < trc1->n_tracers; itr++) + { if (trc1->ts[itr]->chi_min < chi_min1) chi_min1 = trc1->ts[itr]->chi_min; if (trc1->ts[itr]->chi_max > chi_max1) @@ -79,7 +85,8 @@ static void get_k_interval(ccl_cosmology *cosmo, double chi_min2 = 1E15; double chi_max2 = -1E15; - for (itr=0; itr < trc2->n_tracers; itr++) { + for (itr = 0; itr < trc2->n_tracers; itr++) + { if (trc2->ts[itr]->chi_min < chi_min2) chi_min2 = trc2->ts[itr]->chi_min; if (trc2->ts[itr]->chi_max > chi_max2) @@ -92,40 +99,45 @@ static void get_k_interval(ccl_cosmology *cosmo, double chi_max = fmin(chi_max1, chi_max2); if (chi_min <= 0) - chi_min = 0.5*(l+0.5)/cosmo->spline_params.K_MAX; + chi_min = 0.5 * (l + 0.5) / cosmo->spline_params.K_MAX; // Don't go beyond kmax - *lkmax = log(fmin(cosmo->spline_params.K_MAX, 2*(l+0.5)/chi_min)); - *lkmin = log(fmax(cosmo->spline_params.K_MIN, (l+0.5)/chi_max)); + *lkmax = log(fmin(cosmo->spline_params.K_MAX, 2 * (l + 0.5) / chi_min)); + *lkmin = log(fmax(cosmo->spline_params.K_MIN, (l + 0.5) / chi_max)); } static double transfer_limber_single(ccl_cl_tracer_t *tr, double l, double lk, double k, double chi_l, double a_l, ccl_cosmology *cosmo, ccl_f2d_t *psp, int ignore_jbes_deriv, - int *status) { + int *status) +{ double dd = 0; // Kernel and transfer evaluated at chi_l double w = ccl_cl_tracer_t_get_kernel(tr, chi_l, status); - double t = ccl_cl_tracer_t_get_transfer(tr, lk,a_l, status); + double t = ccl_cl_tracer_t_get_transfer(tr, lk, a_l, status); double fl = ccl_cl_tracer_t_get_f_ell(tr, l, status); - if (tr->der_bessel < 1) { //We don't need l+1 - dd = w*t; - if (tr->der_bessel == -1) { //If we divide by (chi*k)^2 - double lp1h = l+0.5; - dd /= (lp1h*lp1h); + if (tr->der_bessel < 1) + { // We don't need l+1 + dd = w * t; + if (tr->der_bessel == -1) + { // If we divide by (chi*k)^2 + double lp1h = l + 0.5; + dd /= (lp1h * lp1h); } } - else { // We will need l+1 - if(ignore_jbes_deriv) + else + { // We will need l+1 + if (ignore_jbes_deriv) dd = 0; - else { + else + { // Compute chi_{l+1} and a_{l+1} - double lp1h = l+0.5; - double lp3h = l+1.5; - double chi_lp = lp3h/k; + double lp1h = l + 0.5; + double lp3h = l + 1.5; + double chi_lp = lp3h / k; double a_lp = ccl_scale_factor_of_chi(cosmo, chi_lp, status); // Compute power spectrum ratio there @@ -134,40 +146,43 @@ static double transfer_limber_single(ccl_cl_tracer_t *tr, double l, double lk, // Compute kernel and trasfer at chi_{l+1} double w_p = ccl_cl_tracer_t_get_kernel(tr, chi_lp, status); - double t_p = ccl_cl_tracer_t_get_transfer(tr, lk,a_lp, status); + double t_p = ccl_cl_tracer_t_get_transfer(tr, lk, a_lp, status); // sqrt(2l+1/2l+3) - double sqell = sqrt(lp1h*pk_ratio/lp3h); + double sqell = sqrt(lp1h * pk_ratio / lp3h); if (tr->der_bessel == 1) - dd = l*w*t/lp1h-sqell*w_p*t_p; - else //we assume der_bessel=2 here to avoid extra if clause - dd = sqell*2*w_p*t_p/lp3h - (0.25+2*l)*w*t/(lp1h*lp1h); + dd = l * w * t / lp1h - sqell * w_p * t_p; + else // we assume der_bessel=2 here to avoid extra if clause + dd = sqell * 2 * w_p * t_p / lp3h - (0.25 + 2 * l) * w * t / (lp1h * lp1h); } } - return dd*fl; + return dd * fl; } -static double transfer_limber_wrap(double l,double lk, double k, double chi, +static double transfer_limber_wrap(double l, double lk, double k, double chi, double a, ccl_cl_tracer_collection_t *trc, - ccl_cosmology *cosmo,ccl_f2d_t *psp, - int ignore_jbes_deriv, int *status) { + ccl_cosmology *cosmo, ccl_f2d_t *psp, + int ignore_jbes_deriv, int *status) +{ int itr; double transfer = 0; - for (itr=0; itr < trc->n_tracers; itr++) { + for (itr = 0; itr < trc->n_tracers; itr++) + { transfer += transfer_limber_single( - trc->ts[itr], l, lk, k, chi, a, cosmo, psp, ignore_jbes_deriv, status); + trc->ts[itr], l, lk, k, chi, a, cosmo, psp, ignore_jbes_deriv, status); if (*status != 0) return -1; } return transfer; } - -static double cl_integrand(double lk, void *params) { + +static double cl_integrand(double lk, void *params) +{ double d1, d2; integ_cl_par *p = (integ_cl_par *)params; double k = exp(lk); - double chi = (p->l+0.5)/k; + double chi = (p->l + 0.5) / k; double a = ccl_scale_factor_of_chi(p->cosmo, chi, p->status); d1 = transfer_limber_wrap(p->l, lk, k, chi, a, p->trc1, @@ -183,39 +198,46 @@ static double cl_integrand(double lk, void *params) { double pk = ccl_f2d_t_eval(p->psp, lk, a, p->cosmo, p->status); - return k*pk*d1*d2; + return k * pk * d1 * d2; } static void integ_cls_limber_spline(ccl_cosmology *cosmo, - integ_cl_par *ipar, - double lkmin, double lkmax, - double *result, int *status) { + integ_cl_par *ipar, + double lkmin, double lkmax, + double *result, int *status) +{ int ik; int nk = (int)(fmax((lkmax - lkmin) / cosmo->spline_params.DLOGK_INTEGRATION + 0.5, - 1))+1; + 1)) + + 1; double *fk_arr = NULL; double *lk_arr = NULL; lk_arr = ccl_linear_spacing(lkmin, lkmax, nk); - if(lk_arr == NULL) + if (lk_arr == NULL) *status = CCL_ERROR_LOGSPACE; - if(*status == 0) { + if (*status == 0) + { fk_arr = malloc(nk * sizeof(double)); - if(fk_arr == NULL) + if (fk_arr == NULL) *status = CCL_ERROR_MEMORY; } - if(*status == 0) { - for(ik=0; ikstatus)) { - *status = *(ipar->status); - break; + if (*(ipar->status)) + { + *status = *(ipar->status); + break; } } } - if(*status == 0) { + if (*status == 0) + { ccl_integ_spline(1, nk, lk_arr, &fk_arr, 1, -1, result, gsl_interp_akima, status); @@ -225,63 +247,67 @@ static void integ_cls_limber_spline(ccl_cosmology *cosmo, } static void integ_cls_limber_qag_quad(ccl_cosmology *cosmo, - gsl_function *F, - double lkmin, double lkmax, - gsl_integration_workspace *w, - double *result, double *eresult, - int *status) { + gsl_function *F, + double lkmin, double lkmax, + gsl_integration_workspace *w, + double *result, double *eresult, + int *status) +{ int gslstatus; size_t nevals; gsl_integration_cquad_workspace *w_cquad = NULL; // Integrate gslstatus = gsl_integration_qag(F, lkmin, lkmax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - cosmo->gsl_params.N_ITERATION, - cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, - w, result, eresult); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + cosmo->gsl_params.N_ITERATION, + cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, + w, result, eresult); // Test if a round-off error occured in the evaluation of the integral // If so, try another integration function, more robust but potentially slower - if (gslstatus == GSL_EROUND) { + if (gslstatus == GSL_EROUND) + { ccl_raise_gsl_warning(gslstatus, - "ccl_cls.c: integ_cls_limber_qag_quad(): " - "Default GSL integration failure, attempting backup method."); + "ccl_cls.c: integ_cls_limber_qag_quad(): " + "Default GSL integration failure, attempting backup method."); w_cquad = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION); if (w_cquad == NULL) *status = CCL_ERROR_MEMORY; - if (*status == 0) { + if (*status == 0) + { nevals = 0; gslstatus = gsl_integration_cquad(F, lkmin, lkmax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - w_cquad, result, eresult, &nevals); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + w_cquad, result, eresult, &nevals); } } gsl_integration_cquad_workspace_free(w_cquad); - if(*status == 0) + if (*status == 0) *status = gslstatus; } void ccl_angular_cls_limber(ccl_cosmology *cosmo, - ccl_cl_tracer_collection_t *trc1, - ccl_cl_tracer_collection_t *trc2, - ccl_f2d_t *psp, - int nl_out, double *l_out, double *cl_out, - ccl_integration_t integration_method, - int *status) { + ccl_cl_tracer_collection_t *trc1, + ccl_cl_tracer_collection_t *trc2, + ccl_f2d_t *psp, + int nl_out, double *l_out, double *cl_out, + ccl_integration_t integration_method, + int *status) +{ // make sure to init core things for safety - if (!cosmo->computed_distances) { + if (!cosmo->computed_distances) + { *status = CCL_ERROR_DISTANCES_INIT; ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cls_limber(): distance splines have not been precomputed!"); + cosmo, + "ccl_cls.c: ccl_angular_cls_limber(): distance splines have not been precomputed!"); return; } - #pragma omp parallel shared(cosmo, trc1, trc2, l_out, cl_out, \ - nl_out, status, psp, integration_method) \ - default(none) +#pragma omp parallel shared(cosmo, trc1, trc2, l_out, cl_out, \ + nl_out, status, psp, integration_method) default(none) { int clastatus, lind; integ_cl_par ipar; @@ -290,7 +316,8 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, gsl_function F; double lkmin, lkmax, l, result, eresult; - if (local_status == 0) { + if (local_status == 0) + { // Set up integrating function parameters ipar.cosmo = cosmo; ipar.trc1 = trc1; @@ -299,24 +326,30 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, ipar.status = &clastatus; } - if(integration_method == ccl_integration_qag_quad) { - if (local_status == 0) { - w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); - if (w == NULL) { - local_status = CCL_ERROR_MEMORY; - } + if (integration_method == ccl_integration_qag_quad) + { + if (local_status == 0) + { + w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); + if (w == NULL) + { + local_status = CCL_ERROR_MEMORY; + } } - if (local_status == 0) { - // Set up integrating function - F.function = &cl_integrand; - F.params = &ipar; + if (local_status == 0) + { + // Set up integrating function + F.function = &cl_integrand; + F.params = &ipar; } } - #pragma omp for schedule(dynamic) - for (lind=0; lind < nl_out; ++lind) { - if (local_status == 0) { +#pragma omp for schedule(dynamic) + for (lind = 0; lind < nl_out; ++lind) + { + if (local_status == 0) + { l = l_out[lind]; clastatus = 0; ipar.l = l; @@ -324,22 +357,26 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, // Get integration limits get_k_interval(cosmo, trc1, trc2, l, &lkmin, &lkmax); - // Integrate - if(integration_method == ccl_integration_qag_quad) { - integ_cls_limber_qag_quad(cosmo, &F, lkmin, lkmax, w, - &result, &eresult, &local_status); - } - else if(integration_method == ccl_integration_spline) { - integ_cls_limber_spline(cosmo, &ipar, lkmin, lkmax, - &result, &local_status); - } - else - local_status = CCL_ERROR_NOT_IMPLEMENTED; - - if ((*ipar.status == 0) && (local_status == 0)) { - cl_out[lind] = result / (l+0.5); + // Integrate + if (integration_method == ccl_integration_qag_quad) + { + integ_cls_limber_qag_quad(cosmo, &F, lkmin, lkmax, w, + &result, &eresult, &local_status); + } + else if (integration_method == ccl_integration_spline) + { + integ_cls_limber_spline(cosmo, &ipar, lkmin, lkmax, + &result, &local_status); + } + else + local_status = CCL_ERROR_NOT_IMPLEMENTED; + + if ((*ipar.status == 0) && (local_status == 0)) + { + cl_out[lind] = result / (l + 0.5); } - else { + else + { ccl_raise_gsl_warning(local_status, "ccl_cls.c: ccl_angular_cls_limber():"); cl_out[lind] = NAN; local_status = CCL_ERROR_INTEG; @@ -349,16 +386,18 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, gsl_integration_workspace_free(w); - if (local_status) { - #pragma omp atomic write + if (local_status) + { +#pragma omp atomic write *status = local_status; } } - if (*status) { + if (*status) + { ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cls_limber(); integration error\n"); + cosmo, + "ccl_cls.c: ccl_angular_cls_limber(); integration error\n"); } } @@ -367,21 +406,22 @@ void ccl_angular_cls_nonlimber(ccl_cosmology *cosmo, ccl_cl_tracer_collection_t *trc2, ccl_f2d_t *psp, int nl_out, int *l_out, double *cl_out, - int *status) { + int *status) +{ *status = CCL_ERROR_INCONSISTENT; ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cls_nonlimber(); non-Limber integrator not implemented yet\n"); + cosmo, + "ccl_cls.c: ccl_angular_cls_nonlimber(); non-Limber integrator not implemented yet\n"); } static double cov_integrand(double chi, void *params) { - double d1, d2, d3, d4, tkk, ker=1; + double d1, d2, d3, d4, tkk, ker = 1; integ_cov_par *p = (integ_cov_par *)params; - double k1=(p->l1+0.5)/chi; - double k2=(p->l2+0.5)/chi; - double lk1=log(k1); - double lk2=log(k2); + double k1 = (p->l1 + 0.5) / chi; + double k2 = (p->l2 + 0.5) / chi; + double lk1 = log(k1); + double lk2 = log(k2); double a = ccl_scale_factor_of_chi(p->cosmo, chi, p->status); d1 = transfer_limber_wrap(p->l1, lk1, k1, chi, a, p->trc1, @@ -403,43 +443,49 @@ static double cov_integrand(double chi, void *params) tkk = ccl_f3d_t_eval(p->tsp, lk1, lk2, a, p->finda, p->cosmo, p->status); - if(p->ker_extra!=NULL) + if (p->ker_extra != NULL) ker = ccl_f1d_t_eval(p->ker_extra, a); - return d1*d2*d3*d4*tkk*ker/pow(chi, p->chipow); + return d1 * d2 * d3 * d4 * tkk * ker / pow(chi, p->chipow); } static void integ_cov_limber_spline(ccl_cosmology *cosmo, - integ_cov_par *ipar, - double chimin, double chimax, - double *result, int *status) + integ_cov_par *ipar, + double chimin, double chimax, + double *result, int *status) { int ichi; int nchi = (int)(fmax((chimax - chimin) / cosmo->spline_params.DCHI_INTEGRATION + 0.5, - 1))+1; + 1)) + + 1; double *fchi_arr = NULL; double *chi_arr = NULL; chi_arr = ccl_linear_spacing(chimin, chimax, nchi); - if(chi_arr == NULL) + if (chi_arr == NULL) *status = CCL_ERROR_LOGSPACE; - if(*status == 0) { + if (*status == 0) + { fchi_arr = malloc(nchi * sizeof(double)); - if(fchi_arr == NULL) + if (fchi_arr == NULL) *status = CCL_ERROR_MEMORY; } - if(*status == 0) { - for(ichi=0; ichistatus)) { - *status = *(ipar->status); - break; + if (*(ipar->status)) + { + *status = *(ipar->status); + break; } } } - if(*status == 0) { + if (*status == 0) + { ccl_integ_spline(1, nchi, chi_arr, &fchi_arr, 1, -1, result, gsl_interp_akima, status); @@ -450,40 +496,43 @@ static void integ_cov_limber_spline(ccl_cosmology *cosmo, } static void integ_cov_limber_qag_quad(ccl_cosmology *cosmo, - gsl_function *F, - double chimin, double chimax, - gsl_integration_workspace *w, - double *result, double *eresult, - int *status) { + gsl_function *F, + double chimin, double chimax, + gsl_integration_workspace *w, + double *result, double *eresult, + int *status) +{ int gslstatus; size_t nevals; gsl_integration_cquad_workspace *w_cquad = NULL; // Integrate gslstatus = gsl_integration_qag(F, chimin, chimax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - cosmo->gsl_params.N_ITERATION, - cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, - w, result, eresult); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + cosmo->gsl_params.N_ITERATION, + cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, + w, result, eresult); // Test if a round-off error occured in the evaluation of the integral // If so, try another integration function, more robust but potentially slower - if (gslstatus == GSL_EROUND) { + if (gslstatus == GSL_EROUND) + { ccl_raise_gsl_warning(gslstatus, - "ccl_cls.c: ccl_angular_cov_limber(): " - "Default GSL integration failure, attempting backup method."); + "ccl_cls.c: ccl_angular_cov_limber(): " + "Default GSL integration failure, attempting backup method."); w_cquad = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION); if (w_cquad == NULL) *status = CCL_ERROR_MEMORY; - if (*status == 0) { + if (*status == 0) + { nevals = 0; gslstatus = gsl_integration_cquad(F, chimin, chimax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - w_cquad, result, eresult, &nevals); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + w_cquad, result, eresult, &nevals); } } gsl_integration_cquad_workspace_free(w_cquad); - if(*status == 0) + if (*status == 0) *status = gslstatus; } @@ -500,21 +549,21 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, int chi_exponent, ccl_f1d_t *kernel_extra, double prefactor_extra, int *status) { - if(!cosmo->computed_distances) { + if (!cosmo->computed_distances) + { *status = CCL_ERROR_DISTANCES_INIT; ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cl_limber(): distance splines have not been precomputed!"); + cosmo, + "ccl_cls.c: ccl_angular_cl_limber(): distance splines have not been precomputed!"); return; } - #pragma omp parallel shared(cosmo, trc1, trc2, trc3, trc4, tsp, \ - nl1_out, l1_out, nl2_out, l2_out, cov_out, \ - integration_method, chi_exponent, \ - kernel_extra, prefactor_extra, status) \ - default(none) +#pragma omp parallel shared(cosmo, trc1, trc2, trc3, trc4, tsp, \ + nl1_out, l1_out, nl2_out, l2_out, cov_out, \ + integration_method, chi_exponent, \ + kernel_extra, prefactor_extra, status) default(none) { - int clastatus, lind1,lind2; + int clastatus, lind1, lind2; integ_cov_par ipar; gsl_integration_workspace *w = NULL; int local_status = *status; @@ -522,7 +571,7 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, double chimin, chimax; double l1, l2, result, eresult; ccl_a_finder *finda = ccl_a_finder_new_from_f3d(tsp); - + // Find integration limits chimin = 1E15; chimax = -1E15; @@ -530,8 +579,9 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, update_chi_limits(trc2, &chimin, &chimax, 0); update_chi_limits(trc3, &chimin, &chimax, 0); update_chi_limits(trc4, &chimin, &chimax, 0); - - if (local_status == 0) { + + if (local_status == 0) + { // Set up integrating function parameters ipar.cosmo = cosmo; ipar.trc1 = trc1; @@ -545,50 +595,61 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, ipar.chipow = chi_exponent; } - if(integration_method == ccl_integration_qag_quad) { - if (local_status == 0) { - w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); - if (w == NULL) { - local_status = CCL_ERROR_MEMORY; - } + if (integration_method == ccl_integration_qag_quad) + { + if (local_status == 0) + { + w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); + if (w == NULL) + { + local_status = CCL_ERROR_MEMORY; + } } - if (local_status == 0) { - // Set up integrating function - F.function = &cov_integrand; - F.params = &ipar; + if (local_status == 0) + { + // Set up integrating function + F.function = &cov_integrand; + F.params = &ipar; } } - #pragma omp for schedule(dynamic) - for (lind1=0; lind1 < nl1_out; ++lind1) { +#pragma omp for schedule(dynamic) + for (lind1 = 0; lind1 < nl1_out; ++lind1) + { l1 = l1_out[lind1]; ipar.l1 = l1; - for (lind2=0; lind2 < nl2_out; ++lind2) { - if (local_status == 0) { + for (lind2 = 0; lind2 < nl2_out; ++lind2) + { + if (local_status == 0) + { l2 = l2_out[lind2]; clastatus = 0; ipar.l2 = l2; // Integrate - if(integration_method == ccl_integration_qag_quad) { + if (integration_method == ccl_integration_qag_quad) + { integ_cov_limber_qag_quad(cosmo, &F, chimin, chimax, w, &result, &eresult, &local_status); } - else if(integration_method == ccl_integration_spline) { + else if (integration_method == ccl_integration_spline) + { integ_cov_limber_spline(cosmo, &ipar, chimin, chimax, &result, &local_status); } else local_status = CCL_ERROR_NOT_IMPLEMENTED; - if ((*ipar.status == 0) && (local_status == 0)) { - cov_out[lind1+nl1_out*lind2] = result * prefactor_extra; + if ((*ipar.status == 0) && (local_status == 0)) + { + cov_out[lind1 + nl1_out * lind2] = result * prefactor_extra; } - else { + else + { ccl_raise_gsl_warning(local_status, "ccl_cls.c: ccl_angular_cov_limber():"); - cov_out[lind1+nl1_out*lind2] = NAN; + cov_out[lind1 + nl1_out * lind2] = NAN; local_status = CCL_ERROR_INTEG; } } @@ -596,17 +657,19 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, } gsl_integration_workspace_free(w); - - if (local_status) { - #pragma omp atomic write - *status = local_status; + + if (local_status) + { +#pragma omp atomic write + *status = local_status; } ccl_a_finder_free(finda); } - if (*status) { + if (*status) + { ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cov_limber(); integration error\n"); + cosmo, + "ccl_cls.c: ccl_angular_cov_limber(); integration error\n"); } } diff --git a/src/ccl_fftlog.c b/src/ccl_fftlog.c index 5efd1765a..c717e75ae 100644 --- a/src/ccl_fftlog.c +++ b/src/ccl_fftlog.c @@ -208,7 +208,7 @@ static void fht(int npk, int N, double dim, double mu, double q, double kcrc, int noring, double complex* u, int *status) { - fftw_plan forward_plan, reverse_plan; + fftw_plan forward_plan = NULL, reverse_plan = NULL; double L = log(k[N-1]/k[0]) * N/(N-1.); double complex* ulocal = NULL; if(u == NULL) { @@ -225,8 +225,8 @@ static void fht(int npk, int N, } } - fftw_complex* a_tmp; - fftw_complex* b_tmp; + fftw_complex* a_tmp = NULL; + fftw_complex* b_tmp = NULL; if(*status == 0) { a_tmp = fftw_alloc_complex(N); if(a_tmp==NULL) @@ -336,18 +336,15 @@ static void fht(int npk, int N, } //end omp parallel } - if(*status == 0) { - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - } + if(forward_plan != NULL) fftw_destroy_plan(forward_plan); + if(reverse_plan != NULL) fftw_destroy_plan(reverse_plan); free(ulocal); - //TODO: free this up fftw_free(a_tmp); fftw_free(b_tmp); } -/* Compute the discrete Hankel transform of the function a(r) +/* Compute the discrete Hankel transform of the function a(r) * weighted by a power law and the nth derivative of the * (spherical) bessel function. Explicitly, this function computes * \tilde{a}(k)= \int \frac{dr}{(xk)^plaw} a(r) (J/j)^(n)_\mu(xk) @@ -369,11 +366,11 @@ static void general_fht(int npk, int N, kcrc = mu+1.0; mu = mu+0.5*spherical_bessel; - fftw_plan forward_plan, reverse_plan; + fftw_plan forward_plan = NULL, reverse_plan = NULL; double L = log(k[N-1]/k[0]) * N/(N-1.); double complex* ulocal = NULL; if(u == NULL) { - + kcrc = goodkr_new_deriv(N, mu, q, L, spherical_bessel, bessel_deriv,plaw, kcrc); ulocal = malloc (sizeof(complex double)*N); @@ -386,8 +383,8 @@ static void general_fht(int npk, int N, } } - fftw_complex* a_tmp; - fftw_complex* b_tmp; + fftw_complex* a_tmp = NULL; + fftw_complex* b_tmp = NULL; if(*status == 0) { a_tmp = fftw_alloc_complex(N); if(a_tmp==NULL) @@ -500,13 +497,10 @@ static void general_fht(int npk, int N, } //end omp parallel } - if(*status == 0) { - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - } + if(forward_plan != NULL) fftw_destroy_plan(forward_plan); + if(reverse_plan != NULL) fftw_destroy_plan(reverse_plan); free(ulocal); - //TODO: free this up fftw_free(a_tmp); fftw_free(b_tmp); } diff --git a/src/ccl_power.c b/src/ccl_power.c index 709887b32..7e1c76fde 100644 --- a/src/ccl_power.c +++ b/src/ccl_power.c @@ -129,10 +129,12 @@ static ccl_f2d_t *ccl_compute_linpower_analytic(ccl_cosmology* cosmo, void* par, y2d[j] += log_sigma8; } + // Free the first spline unconditionally; re-allocate with normalization only on success. + // Without this, a ccl_sigma8 failure above would leak the first psp_out. + ccl_f2d_t_free(psp_out); + psp_out = NULL; + if(*status==0) { - // Free the previous P(k,a) spline, and allocate a new one to store the - // properly-normalized P(k,a) - ccl_f2d_t_free(psp_out); psp_out = ccl_f2d_t_new(na,z,nk,x,y2d,NULL,NULL,0, 1,2,ccl_f2d_cclgrowth,1,0,2, ccl_f2d_3,status); From 05b9d79a1509bcdebbee0b0d6bb502ea5662d515 Mon Sep 17 00:00:00 2001 From: Robert Reischke Date: Mon, 11 May 2026 13:24:45 +0200 Subject: [PATCH 2/4] update --- pyccl/halos/pk_2pt.py | 1 + pyccl/halos/profiles/profile_base.py | 1 + src/ccl_cls.c | 463 ++++++++++++--------------- 3 files changed, 202 insertions(+), 263 deletions(-) diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py index 5260584f7..0eb049b2e 100644 --- a/pyccl/halos/pk_2pt.py +++ b/pyccl/halos/pk_2pt.py @@ -102,6 +102,7 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *, for ia, aa in enumerate(a_use): # normalizations norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) + if prof2 == prof: norm2 = norm1 else: diff --git a/pyccl/halos/profiles/profile_base.py b/pyccl/halos/profiles/profile_base.py index bf56fc6a9..2c9cdbc84 100644 --- a/pyccl/halos/profiles/profile_base.py +++ b/pyccl/halos/profiles/profile_base.py @@ -516,5 +516,6 @@ def get_normalization(self, cosmo, a, *, hmc=None): class HaloProfilePressure(HaloProfile): """Base for pressure halo profiles.""" + class HaloProfileCIB(HaloProfile): """Base for CIB halo profiles.""" diff --git a/src/ccl_cls.c b/src/ccl_cls.c index 35fec16f2..e32cfeea0 100644 --- a/src/ccl_cls.c +++ b/src/ccl_cls.c @@ -8,8 +8,7 @@ #include "ccl.h" -typedef struct -{ +typedef struct{ double l; ccl_cosmology *cosmo; ccl_cl_tracer_collection_t *trc1; @@ -18,8 +17,8 @@ typedef struct int *status; } integ_cl_par; -typedef struct -{ + +typedef struct{ int chipow; double l1; double l2; @@ -39,28 +38,25 @@ static void update_chi_limits(ccl_cl_tracer_collection_t *trc, int is_union) { int itr; - double chimin_h = 1E15; - double chimax_h = -1E15; - for (itr = 0; itr < trc->n_tracers; itr++) - { + double chimin_h=1E15; + double chimax_h=-1E15; + for(itr=0; itr < trc->n_tracers; itr++) { if (trc->ts[itr]->chi_min < chimin_h) chimin_h = trc->ts[itr]->chi_min; if (trc->ts[itr]->chi_max > chimax_h) chimax_h = trc->ts[itr]->chi_max; } - if (is_union) - { - if (chimin_h < *chimin) + if(is_union) { + if(chimin_h < *chimin) *chimin = chimin_h; - if (chimax_h > *chimax) + if(chimax_h > *chimax) *chimax = chimax_h; } - else - { - if (chimin_h > *chimin) + else { + if(chimin_h > *chimin) *chimin = chimin_h; - if (chimax_h < *chimax) + if(chimax_h < *chimax) *chimax = chimax_h; } } @@ -68,15 +64,13 @@ static void update_chi_limits(ccl_cl_tracer_collection_t *trc, static void get_k_interval(ccl_cosmology *cosmo, ccl_cl_tracer_collection_t *trc1, ccl_cl_tracer_collection_t *trc2, - double l, double *lkmin, double *lkmax) -{ + double l, double *lkmin, double *lkmax) { int itr; // Loop through all tracers and find distance bounds double chi_min1 = 1E15; double chi_max1 = -1E15; - for (itr = 0; itr < trc1->n_tracers; itr++) - { + for (itr=0; itr < trc1->n_tracers; itr++) { if (trc1->ts[itr]->chi_min < chi_min1) chi_min1 = trc1->ts[itr]->chi_min; if (trc1->ts[itr]->chi_max > chi_max1) @@ -85,8 +79,7 @@ static void get_k_interval(ccl_cosmology *cosmo, double chi_min2 = 1E15; double chi_max2 = -1E15; - for (itr = 0; itr < trc2->n_tracers; itr++) - { + for (itr=0; itr < trc2->n_tracers; itr++) { if (trc2->ts[itr]->chi_min < chi_min2) chi_min2 = trc2->ts[itr]->chi_min; if (trc2->ts[itr]->chi_max > chi_max2) @@ -99,45 +92,40 @@ static void get_k_interval(ccl_cosmology *cosmo, double chi_max = fmin(chi_max1, chi_max2); if (chi_min <= 0) - chi_min = 0.5 * (l + 0.5) / cosmo->spline_params.K_MAX; + chi_min = 0.5*(l+0.5)/cosmo->spline_params.K_MAX; // Don't go beyond kmax - *lkmax = log(fmin(cosmo->spline_params.K_MAX, 2 * (l + 0.5) / chi_min)); - *lkmin = log(fmax(cosmo->spline_params.K_MIN, (l + 0.5) / chi_max)); + *lkmax = log(fmin(cosmo->spline_params.K_MAX, 2*(l+0.5)/chi_min)); + *lkmin = log(fmax(cosmo->spline_params.K_MIN, (l+0.5)/chi_max)); } static double transfer_limber_single(ccl_cl_tracer_t *tr, double l, double lk, double k, double chi_l, double a_l, ccl_cosmology *cosmo, ccl_f2d_t *psp, int ignore_jbes_deriv, - int *status) -{ + int *status) { double dd = 0; // Kernel and transfer evaluated at chi_l double w = ccl_cl_tracer_t_get_kernel(tr, chi_l, status); - double t = ccl_cl_tracer_t_get_transfer(tr, lk, a_l, status); + double t = ccl_cl_tracer_t_get_transfer(tr, lk,a_l, status); double fl = ccl_cl_tracer_t_get_f_ell(tr, l, status); - if (tr->der_bessel < 1) - { // We don't need l+1 - dd = w * t; - if (tr->der_bessel == -1) - { // If we divide by (chi*k)^2 - double lp1h = l + 0.5; - dd /= (lp1h * lp1h); + if (tr->der_bessel < 1) { //We don't need l+1 + dd = w*t; + if (tr->der_bessel == -1) { //If we divide by (chi*k)^2 + double lp1h = l+0.5; + dd /= (lp1h*lp1h); } } - else - { // We will need l+1 - if (ignore_jbes_deriv) + else { // We will need l+1 + if(ignore_jbes_deriv) dd = 0; - else - { + else { // Compute chi_{l+1} and a_{l+1} - double lp1h = l + 0.5; - double lp3h = l + 1.5; - double chi_lp = lp3h / k; + double lp1h = l+0.5; + double lp3h = l+1.5; + double chi_lp = lp3h/k; double a_lp = ccl_scale_factor_of_chi(cosmo, chi_lp, status); // Compute power spectrum ratio there @@ -146,43 +134,40 @@ static double transfer_limber_single(ccl_cl_tracer_t *tr, double l, double lk, // Compute kernel and trasfer at chi_{l+1} double w_p = ccl_cl_tracer_t_get_kernel(tr, chi_lp, status); - double t_p = ccl_cl_tracer_t_get_transfer(tr, lk, a_lp, status); + double t_p = ccl_cl_tracer_t_get_transfer(tr, lk,a_lp, status); // sqrt(2l+1/2l+3) - double sqell = sqrt(lp1h * pk_ratio / lp3h); + double sqell = sqrt(lp1h*pk_ratio/lp3h); if (tr->der_bessel == 1) - dd = l * w * t / lp1h - sqell * w_p * t_p; - else // we assume der_bessel=2 here to avoid extra if clause - dd = sqell * 2 * w_p * t_p / lp3h - (0.25 + 2 * l) * w * t / (lp1h * lp1h); + dd = l*w*t/lp1h-sqell*w_p*t_p; + else //we assume der_bessel=2 here to avoid extra if clause + dd = sqell*2*w_p*t_p/lp3h - (0.25+2*l)*w*t/(lp1h*lp1h); } } - return dd * fl; + return dd*fl; } -static double transfer_limber_wrap(double l, double lk, double k, double chi, +static double transfer_limber_wrap(double l,double lk, double k, double chi, double a, ccl_cl_tracer_collection_t *trc, - ccl_cosmology *cosmo, ccl_f2d_t *psp, - int ignore_jbes_deriv, int *status) -{ + ccl_cosmology *cosmo,ccl_f2d_t *psp, + int ignore_jbes_deriv, int *status) { int itr; double transfer = 0; - for (itr = 0; itr < trc->n_tracers; itr++) - { + for (itr=0; itr < trc->n_tracers; itr++) { transfer += transfer_limber_single( - trc->ts[itr], l, lk, k, chi, a, cosmo, psp, ignore_jbes_deriv, status); + trc->ts[itr], l, lk, k, chi, a, cosmo, psp, ignore_jbes_deriv, status); if (*status != 0) return -1; } return transfer; } - -static double cl_integrand(double lk, void *params) -{ + +static double cl_integrand(double lk, void *params) { double d1, d2; integ_cl_par *p = (integ_cl_par *)params; double k = exp(lk); - double chi = (p->l + 0.5) / k; + double chi = (p->l+0.5)/k; double a = ccl_scale_factor_of_chi(p->cosmo, chi, p->status); d1 = transfer_limber_wrap(p->l, lk, k, chi, a, p->trc1, @@ -198,46 +183,39 @@ static double cl_integrand(double lk, void *params) double pk = ccl_f2d_t_eval(p->psp, lk, a, p->cosmo, p->status); - return k * pk * d1 * d2; + return k*pk*d1*d2; } static void integ_cls_limber_spline(ccl_cosmology *cosmo, - integ_cl_par *ipar, - double lkmin, double lkmax, - double *result, int *status) -{ + integ_cl_par *ipar, + double lkmin, double lkmax, + double *result, int *status) { int ik; int nk = (int)(fmax((lkmax - lkmin) / cosmo->spline_params.DLOGK_INTEGRATION + 0.5, - 1)) + - 1; + 1))+1; double *fk_arr = NULL; double *lk_arr = NULL; lk_arr = ccl_linear_spacing(lkmin, lkmax, nk); - if (lk_arr == NULL) + if(lk_arr == NULL) *status = CCL_ERROR_LOGSPACE; - if (*status == 0) - { + if(*status == 0) { fk_arr = malloc(nk * sizeof(double)); - if (fk_arr == NULL) + if(fk_arr == NULL) *status = CCL_ERROR_MEMORY; } - if (*status == 0) - { - for (ik = 0; ik < nk; ik++) - { + if(*status == 0) { + for(ik=0; ikstatus)) - { - *status = *(ipar->status); - break; + if(*(ipar->status)) { + *status = *(ipar->status); + break; } } } - if (*status == 0) - { + if(*status == 0) { ccl_integ_spline(1, nk, lk_arr, &fk_arr, 1, -1, result, gsl_interp_akima, status); @@ -247,67 +225,63 @@ static void integ_cls_limber_spline(ccl_cosmology *cosmo, } static void integ_cls_limber_qag_quad(ccl_cosmology *cosmo, - gsl_function *F, - double lkmin, double lkmax, - gsl_integration_workspace *w, - double *result, double *eresult, - int *status) -{ + gsl_function *F, + double lkmin, double lkmax, + gsl_integration_workspace *w, + double *result, double *eresult, + int *status) { int gslstatus; size_t nevals; gsl_integration_cquad_workspace *w_cquad = NULL; // Integrate gslstatus = gsl_integration_qag(F, lkmin, lkmax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - cosmo->gsl_params.N_ITERATION, - cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, - w, result, eresult); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + cosmo->gsl_params.N_ITERATION, + cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, + w, result, eresult); // Test if a round-off error occured in the evaluation of the integral // If so, try another integration function, more robust but potentially slower - if (gslstatus == GSL_EROUND) - { + if (gslstatus == GSL_EROUND) { ccl_raise_gsl_warning(gslstatus, - "ccl_cls.c: integ_cls_limber_qag_quad(): " - "Default GSL integration failure, attempting backup method."); + "ccl_cls.c: integ_cls_limber_qag_quad(): " + "Default GSL integration failure, attempting backup method."); w_cquad = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION); if (w_cquad == NULL) *status = CCL_ERROR_MEMORY; - if (*status == 0) - { + if (*status == 0) { nevals = 0; gslstatus = gsl_integration_cquad(F, lkmin, lkmax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - w_cquad, result, eresult, &nevals); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + w_cquad, result, eresult, &nevals); } } gsl_integration_cquad_workspace_free(w_cquad); - if (*status == 0) + if(*status == 0) *status = gslstatus; } void ccl_angular_cls_limber(ccl_cosmology *cosmo, - ccl_cl_tracer_collection_t *trc1, - ccl_cl_tracer_collection_t *trc2, - ccl_f2d_t *psp, - int nl_out, double *l_out, double *cl_out, - ccl_integration_t integration_method, - int *status) -{ + ccl_cl_tracer_collection_t *trc1, + ccl_cl_tracer_collection_t *trc2, + ccl_f2d_t *psp, + int nl_out, double *l_out, double *cl_out, + ccl_integration_t integration_method, + int *status) { // make sure to init core things for safety - if (!cosmo->computed_distances) - { + if (!cosmo->computed_distances) { *status = CCL_ERROR_DISTANCES_INIT; ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cls_limber(): distance splines have not been precomputed!"); + cosmo, + "ccl_cls.c: ccl_angular_cls_limber(): distance splines have not been precomputed!"); return; } -#pragma omp parallel shared(cosmo, trc1, trc2, l_out, cl_out, \ - nl_out, status, psp, integration_method) default(none) + #pragma omp parallel shared(cosmo, trc1, trc2, l_out, cl_out, \ + nl_out, status, psp, integration_method) \ + default(none) { int clastatus, lind; integ_cl_par ipar; @@ -316,8 +290,7 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, gsl_function F; double lkmin, lkmax, l, result, eresult; - if (local_status == 0) - { + if (local_status == 0) { // Set up integrating function parameters ipar.cosmo = cosmo; ipar.trc1 = trc1; @@ -326,30 +299,24 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, ipar.status = &clastatus; } - if (integration_method == ccl_integration_qag_quad) - { - if (local_status == 0) - { - w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); - if (w == NULL) - { - local_status = CCL_ERROR_MEMORY; - } + if(integration_method == ccl_integration_qag_quad) { + if (local_status == 0) { + w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); + if (w == NULL) { + local_status = CCL_ERROR_MEMORY; + } } - if (local_status == 0) - { - // Set up integrating function - F.function = &cl_integrand; - F.params = &ipar; + if (local_status == 0) { + // Set up integrating function + F.function = &cl_integrand; + F.params = &ipar; } } -#pragma omp for schedule(dynamic) - for (lind = 0; lind < nl_out; ++lind) - { - if (local_status == 0) - { + #pragma omp for schedule(dynamic) + for (lind=0; lind < nl_out; ++lind) { + if (local_status == 0) { l = l_out[lind]; clastatus = 0; ipar.l = l; @@ -357,26 +324,22 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, // Get integration limits get_k_interval(cosmo, trc1, trc2, l, &lkmin, &lkmax); - // Integrate - if (integration_method == ccl_integration_qag_quad) - { - integ_cls_limber_qag_quad(cosmo, &F, lkmin, lkmax, w, - &result, &eresult, &local_status); + // Integrate + if(integration_method == ccl_integration_qag_quad) { + integ_cls_limber_qag_quad(cosmo, &F, lkmin, lkmax, w, + &result, &eresult, &local_status); + } + else if(integration_method == ccl_integration_spline) { + integ_cls_limber_spline(cosmo, &ipar, lkmin, lkmax, + &result, &local_status); + } + else + local_status = CCL_ERROR_NOT_IMPLEMENTED; + + if ((*ipar.status == 0) && (local_status == 0)) { + cl_out[lind] = result / (l+0.5); } - else if (integration_method == ccl_integration_spline) - { - integ_cls_limber_spline(cosmo, &ipar, lkmin, lkmax, - &result, &local_status); - } - else - local_status = CCL_ERROR_NOT_IMPLEMENTED; - - if ((*ipar.status == 0) && (local_status == 0)) - { - cl_out[lind] = result / (l + 0.5); - } - else - { + else { ccl_raise_gsl_warning(local_status, "ccl_cls.c: ccl_angular_cls_limber():"); cl_out[lind] = NAN; local_status = CCL_ERROR_INTEG; @@ -386,18 +349,16 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo, gsl_integration_workspace_free(w); - if (local_status) - { -#pragma omp atomic write + if (local_status) { + #pragma omp atomic write *status = local_status; } } - if (*status) - { + if (*status) { ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cls_limber(); integration error\n"); + cosmo, + "ccl_cls.c: ccl_angular_cls_limber(); integration error\n"); } } @@ -406,22 +367,21 @@ void ccl_angular_cls_nonlimber(ccl_cosmology *cosmo, ccl_cl_tracer_collection_t *trc2, ccl_f2d_t *psp, int nl_out, int *l_out, double *cl_out, - int *status) -{ + int *status) { *status = CCL_ERROR_INCONSISTENT; ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cls_nonlimber(); non-Limber integrator not implemented yet\n"); + cosmo, + "ccl_cls.c: ccl_angular_cls_nonlimber(); non-Limber integrator not implemented yet\n"); } static double cov_integrand(double chi, void *params) { - double d1, d2, d3, d4, tkk, ker = 1; + double d1, d2, d3, d4, tkk, ker=1; integ_cov_par *p = (integ_cov_par *)params; - double k1 = (p->l1 + 0.5) / chi; - double k2 = (p->l2 + 0.5) / chi; - double lk1 = log(k1); - double lk2 = log(k2); + double k1=(p->l1+0.5)/chi; + double k2=(p->l2+0.5)/chi; + double lk1=log(k1); + double lk2=log(k2); double a = ccl_scale_factor_of_chi(p->cosmo, chi, p->status); d1 = transfer_limber_wrap(p->l1, lk1, k1, chi, a, p->trc1, @@ -443,49 +403,43 @@ static double cov_integrand(double chi, void *params) tkk = ccl_f3d_t_eval(p->tsp, lk1, lk2, a, p->finda, p->cosmo, p->status); - if (p->ker_extra != NULL) + if(p->ker_extra!=NULL) ker = ccl_f1d_t_eval(p->ker_extra, a); - return d1 * d2 * d3 * d4 * tkk * ker / pow(chi, p->chipow); + return d1*d2*d3*d4*tkk*ker/pow(chi, p->chipow); } static void integ_cov_limber_spline(ccl_cosmology *cosmo, - integ_cov_par *ipar, - double chimin, double chimax, - double *result, int *status) + integ_cov_par *ipar, + double chimin, double chimax, + double *result, int *status) { int ichi; int nchi = (int)(fmax((chimax - chimin) / cosmo->spline_params.DCHI_INTEGRATION + 0.5, - 1)) + - 1; + 1))+1; double *fchi_arr = NULL; double *chi_arr = NULL; chi_arr = ccl_linear_spacing(chimin, chimax, nchi); - if (chi_arr == NULL) + if(chi_arr == NULL) *status = CCL_ERROR_LOGSPACE; - if (*status == 0) - { + if(*status == 0) { fchi_arr = malloc(nchi * sizeof(double)); - if (fchi_arr == NULL) + if(fchi_arr == NULL) *status = CCL_ERROR_MEMORY; } - if (*status == 0) - { - for (ichi = 0; ichi < nchi; ichi++) - { + if(*status == 0) { + for(ichi=0; ichistatus)) - { - *status = *(ipar->status); - break; + if(*(ipar->status)) { + *status = *(ipar->status); + break; } } } - if (*status == 0) - { + if(*status == 0) { ccl_integ_spline(1, nchi, chi_arr, &fchi_arr, 1, -1, result, gsl_interp_akima, status); @@ -496,43 +450,40 @@ static void integ_cov_limber_spline(ccl_cosmology *cosmo, } static void integ_cov_limber_qag_quad(ccl_cosmology *cosmo, - gsl_function *F, - double chimin, double chimax, - gsl_integration_workspace *w, - double *result, double *eresult, - int *status) -{ + gsl_function *F, + double chimin, double chimax, + gsl_integration_workspace *w, + double *result, double *eresult, + int *status) { int gslstatus; size_t nevals; gsl_integration_cquad_workspace *w_cquad = NULL; // Integrate gslstatus = gsl_integration_qag(F, chimin, chimax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - cosmo->gsl_params.N_ITERATION, - cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, - w, result, eresult); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + cosmo->gsl_params.N_ITERATION, + cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS, + w, result, eresult); // Test if a round-off error occured in the evaluation of the integral // If so, try another integration function, more robust but potentially slower - if (gslstatus == GSL_EROUND) - { + if (gslstatus == GSL_EROUND) { ccl_raise_gsl_warning(gslstatus, - "ccl_cls.c: ccl_angular_cov_limber(): " - "Default GSL integration failure, attempting backup method."); + "ccl_cls.c: ccl_angular_cov_limber(): " + "Default GSL integration failure, attempting backup method."); w_cquad = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION); if (w_cquad == NULL) *status = CCL_ERROR_MEMORY; - if (*status == 0) - { + if (*status == 0) { nevals = 0; gslstatus = gsl_integration_cquad(F, chimin, chimax, 0, - cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, - w_cquad, result, eresult, &nevals); + cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL, + w_cquad, result, eresult, &nevals); } } gsl_integration_cquad_workspace_free(w_cquad); - if (*status == 0) + if(*status == 0) *status = gslstatus; } @@ -549,21 +500,21 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, int chi_exponent, ccl_f1d_t *kernel_extra, double prefactor_extra, int *status) { - if (!cosmo->computed_distances) - { + if(!cosmo->computed_distances) { *status = CCL_ERROR_DISTANCES_INIT; ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cl_limber(): distance splines have not been precomputed!"); + cosmo, + "ccl_cls.c: ccl_angular_cl_limber(): distance splines have not been precomputed!"); return; } -#pragma omp parallel shared(cosmo, trc1, trc2, trc3, trc4, tsp, \ - nl1_out, l1_out, nl2_out, l2_out, cov_out, \ - integration_method, chi_exponent, \ - kernel_extra, prefactor_extra, status) default(none) + #pragma omp parallel shared(cosmo, trc1, trc2, trc3, trc4, tsp, \ + nl1_out, l1_out, nl2_out, l2_out, cov_out, \ + integration_method, chi_exponent, \ + kernel_extra, prefactor_extra, status) \ + default(none) { - int clastatus, lind1, lind2; + int clastatus, lind1,lind2; integ_cov_par ipar; gsl_integration_workspace *w = NULL; int local_status = *status; @@ -571,7 +522,7 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, double chimin, chimax; double l1, l2, result, eresult; ccl_a_finder *finda = ccl_a_finder_new_from_f3d(tsp); - + // Find integration limits chimin = 1E15; chimax = -1E15; @@ -579,9 +530,8 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, update_chi_limits(trc2, &chimin, &chimax, 0); update_chi_limits(trc3, &chimin, &chimax, 0); update_chi_limits(trc4, &chimin, &chimax, 0); - - if (local_status == 0) - { + + if (local_status == 0) { // Set up integrating function parameters ipar.cosmo = cosmo; ipar.trc1 = trc1; @@ -595,61 +545,50 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, ipar.chipow = chi_exponent; } - if (integration_method == ccl_integration_qag_quad) - { - if (local_status == 0) - { - w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); - if (w == NULL) - { - local_status = CCL_ERROR_MEMORY; - } + if(integration_method == ccl_integration_qag_quad) { + if (local_status == 0) { + w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION); + if (w == NULL) { + local_status = CCL_ERROR_MEMORY; + } } - if (local_status == 0) - { - // Set up integrating function - F.function = &cov_integrand; - F.params = &ipar; + if (local_status == 0) { + // Set up integrating function + F.function = &cov_integrand; + F.params = &ipar; } } -#pragma omp for schedule(dynamic) - for (lind1 = 0; lind1 < nl1_out; ++lind1) - { + #pragma omp for schedule(dynamic) + for (lind1=0; lind1 < nl1_out; ++lind1) { l1 = l1_out[lind1]; ipar.l1 = l1; - for (lind2 = 0; lind2 < nl2_out; ++lind2) - { - if (local_status == 0) - { + for (lind2=0; lind2 < nl2_out; ++lind2) { + if (local_status == 0) { l2 = l2_out[lind2]; clastatus = 0; ipar.l2 = l2; // Integrate - if (integration_method == ccl_integration_qag_quad) - { + if(integration_method == ccl_integration_qag_quad) { integ_cov_limber_qag_quad(cosmo, &F, chimin, chimax, w, &result, &eresult, &local_status); } - else if (integration_method == ccl_integration_spline) - { + else if(integration_method == ccl_integration_spline) { integ_cov_limber_spline(cosmo, &ipar, chimin, chimax, &result, &local_status); } else local_status = CCL_ERROR_NOT_IMPLEMENTED; - if ((*ipar.status == 0) && (local_status == 0)) - { - cov_out[lind1 + nl1_out * lind2] = result * prefactor_extra; + if ((*ipar.status == 0) && (local_status == 0)) { + cov_out[lind1+nl1_out*lind2] = result * prefactor_extra; } - else - { + else { ccl_raise_gsl_warning(local_status, "ccl_cls.c: ccl_angular_cov_limber():"); - cov_out[lind1 + nl1_out * lind2] = NAN; + cov_out[lind1+nl1_out*lind2] = NAN; local_status = CCL_ERROR_INTEG; } } @@ -657,19 +596,17 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, } gsl_integration_workspace_free(w); - - if (local_status) - { -#pragma omp atomic write - *status = local_status; + + if (local_status) { + #pragma omp atomic write + *status = local_status; } ccl_a_finder_free(finda); } - if (*status) - { + if (*status) { ccl_cosmology_set_status_message( - cosmo, - "ccl_cls.c: ccl_angular_cov_limber(); integration error\n"); + cosmo, + "ccl_cls.c: ccl_angular_cov_limber(); integration error\n"); } -} +} \ No newline at end of file From 79c7f600e8e6def73691cd65280262423aad9016 Mon Sep 17 00:00:00 2001 From: Robert Reischke Date: Mon, 11 May 2026 13:30:58 +0200 Subject: [PATCH 3/4] fix --- pyccl/halos/pk_2pt.py | 4 ++-- src/ccl_fftlog.c | 2 ++ src/ccl_power.c | 4 ++-- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py index 0eb049b2e..4bf719a75 100644 --- a/pyccl/halos/pk_2pt.py +++ b/pyccl/halos/pk_2pt.py @@ -102,7 +102,7 @@ def halomod_power_spectrum(cosmo, hmc, k, a, prof, *, for ia, aa in enumerate(a_use): # normalizations norm1 = prof.get_normalization(cosmo, aa, hmc=hmc) - + if prof2 == prof: norm2 = norm1 else: @@ -229,4 +229,4 @@ def halomod_Pk2D(cosmo, hmc, prof, *, return Pk2D(a_arr=a_arr, lk_arr=lk_arr, pk_arr=pk_arr, extrap_order_lok=extrap_order_lok, extrap_order_hik=extrap_order_hik, - is_logp=False) + is_logp=False) \ No newline at end of file diff --git a/src/ccl_fftlog.c b/src/ccl_fftlog.c index c717e75ae..4c72db055 100644 --- a/src/ccl_fftlog.c +++ b/src/ccl_fftlog.c @@ -340,6 +340,7 @@ static void fht(int npk, int N, if(reverse_plan != NULL) fftw_destroy_plan(reverse_plan); free(ulocal); + //TODO: free this up fftw_free(a_tmp); fftw_free(b_tmp); } @@ -501,6 +502,7 @@ static void general_fht(int npk, int N, if(reverse_plan != NULL) fftw_destroy_plan(reverse_plan); free(ulocal); + //TODO: free this up fftw_free(a_tmp); fftw_free(b_tmp); } diff --git a/src/ccl_power.c b/src/ccl_power.c index 7e1c76fde..970aa8134 100644 --- a/src/ccl_power.c +++ b/src/ccl_power.c @@ -129,8 +129,8 @@ static ccl_f2d_t *ccl_compute_linpower_analytic(ccl_cosmology* cosmo, void* par, y2d[j] += log_sigma8; } - // Free the first spline unconditionally; re-allocate with normalization only on success. - // Without this, a ccl_sigma8 failure above would leak the first psp_out. + // Free the first spline unconditionally; re-allocate with normalisation only on success. + // Without this, ccl_sigma8 failure above would leak the first psp_out. ccl_f2d_t_free(psp_out); psp_out = NULL; From a555e1af38b72f3f897dd4e970c99f0ee709f6c3 Mon Sep 17 00:00:00 2001 From: carlosggarcia Date: Wed, 13 May 2026 16:49:28 -0400 Subject: [PATCH 4/4] flaked --- pyccl/halos/pk_2pt.py | 2 +- src/ccl_cls.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyccl/halos/pk_2pt.py b/pyccl/halos/pk_2pt.py index 4bf719a75..6a3ce9af9 100644 --- a/pyccl/halos/pk_2pt.py +++ b/pyccl/halos/pk_2pt.py @@ -229,4 +229,4 @@ def halomod_Pk2D(cosmo, hmc, prof, *, return Pk2D(a_arr=a_arr, lk_arr=lk_arr, pk_arr=pk_arr, extrap_order_lok=extrap_order_lok, extrap_order_hik=extrap_order_hik, - is_logp=False) \ No newline at end of file + is_logp=False) diff --git a/src/ccl_cls.c b/src/ccl_cls.c index e32cfeea0..0b886a229 100644 --- a/src/ccl_cls.c +++ b/src/ccl_cls.c @@ -609,4 +609,4 @@ void ccl_angular_cl_covariance(ccl_cosmology *cosmo, cosmo, "ccl_cls.c: ccl_angular_cov_limber(); integration error\n"); } -} \ No newline at end of file +}