Skip to content

Commit

Permalink
NumCosmo product file (#143)
Browse files Browse the repository at this point in the history 
* Introducing the product-file options.
* Updated Halofit to return linear Pk when all required scales are linear.
* Support for exploration phase in APES.
* Adding calibrate option to numcosmo app.
* Reorganized numcosmo app.
* Removed calibrate_apes tool.
* Testing power-spectra with/without halofit.
* More tests for linear universe.
* Added tests for APES exploration.
* Testing APES MPI.
* Removed unused code and optimizing tests.
  • Loading branch information
@vitenti
vitenti committed Feb 15, 2024
1 parent 5c3f38f commit 31a590e
Show file tree
Hide file tree
Showing 23 changed files with 2,654 additions and 1,268 deletions.
5 changes: 3 additions & 2 deletions notebooks/extern_lib_comparisons/ccl_background.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -428,9 +428,10 @@
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3"
"pygments_lexer": "ipython3",
"version": "3.11.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}
36 changes: 34 additions & 2 deletions numcosmo/math/ncm_fit_esmcmc_walker_apes.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,7 @@ typedef struct _NcmFitESMCMCWalkerAPESPrivate
gboolean use_interp;
gboolean use_threads;
gboolean constructed;
guint exploration;
} NcmFitESMCMCWalkerAPESPrivate;

struct _NcmFitESMCMCWalkerAPES
Expand Down Expand Up @@ -147,6 +148,7 @@ ncm_fit_esmcmc_walker_apes_init (NcmFitESMCMCWalkerAPES *apes)
self->use_interp = FALSE;
self->use_threads = FALSE;
self->constructed = FALSE;
self->exploration = 0;

g_ptr_array_set_free_func (self->thetastar, (GDestroyNotify) ncm_vector_free);
}
Expand Down Expand Up @@ -557,6 +559,9 @@ _ncm_fit_esmcmc_walker_apes_setup (NcmFitESMCMCWalker *walker, NcmMSet *mset, GP
} while (!ncm_mset_fparam_valid_bounds (mset, thetastar_i));
}
}

if (self->exploration > 0)
self->exploration--;
}

static void
Expand Down Expand Up @@ -614,7 +619,10 @@ _ncm_fit_esmcmc_walker_apes_prob (NcmFitESMCMCWalker *walker, GPtrArray *theta,
* exp (- 0.5 * (m2lnL_star - m2lnL_cur)), exp (- 0.5 * (m2lnp_cur - m2lnp_star)));
*/

return exp (-0.5 * ((m2lnL_star - m2lnp_star) - (m2lnL_cur - m2lnp_cur)));
if (self->exploration)
return exp (-0.5 * (m2lnL_star - m2lnL_cur));
else
return exp (-0.5 * ((m2lnL_star - m2lnp_star) - (m2lnL_cur - m2lnp_cur)));
}

static gdouble
Expand All @@ -625,7 +633,10 @@ _ncm_fit_esmcmc_walker_apes_prob_norm (NcmFitESMCMCWalker *walker, GPtrArray *th
const gdouble m2lnp_star = ncm_vector_get (self->m2lnp_star, k);
const gdouble m2lnp_cur = ncm_vector_get (self->m2lnp_cur, k);

return -0.5 * (m2lnp_cur - m2lnp_star);
if (self->exploration)
return 0.0;
else
return -0.5 * (m2lnp_cur - m2lnp_star);
}

static void
Expand Down Expand Up @@ -1108,3 +1119,24 @@ ncm_fit_esmcmc_walker_apes_set_cov_robust (NcmFitESMCMCWalkerAPES *apes)
ncm_stats_dist_kde_set_cov_type (NCM_STATS_DIST_KDE (self->sd1), NCM_STATS_DIST_KDE_COV_TYPE_ROBUST);
}

/**
* ncm_fit_esmcmc_walker_apes_set_exploration:
* @apes: a #NcmFitESMCMCWalkerAPES
* @exploration: a guint
*
* Sets the exploration parameter to be used when building the
* posterior approximation. During the exploration phase, the
* new samples are accepted considering only the posterior.
* This makes the exploration phase to be more efficient to find
* the global maximum of the posterior, but this phase should be
* discarded in the final analysis.
*
*/
void
ncm_fit_esmcmc_walker_apes_set_exploration (NcmFitESMCMCWalkerAPES *apes, guint exploration)
{
NcmFitESMCMCWalkerAPESPrivate * const self = ncm_fit_esmcmc_walker_apes_get_instance_private (apes);

self->exploration = exploration;
}

2 changes: 2 additions & 0 deletions numcosmo/math/ncm_fit_esmcmc_walker_apes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,8 @@ void ncm_fit_esmcmc_walker_apes_set_cov_fixed_from_mset (NcmFitESMCMCWalkerAPES
void ncm_fit_esmcmc_walker_apes_set_cov_robust_diag (NcmFitESMCMCWalkerAPES *apes);
void ncm_fit_esmcmc_walker_apes_set_cov_robust (NcmFitESMCMCWalkerAPES *apes);

void ncm_fit_esmcmc_walker_apes_set_exploration (NcmFitESMCMCWalkerAPES *apes, guint exploration);

G_END_DECLS

#endif /* _NCM_FIT_ESMCMC_WALKER_APES_H_ */
Expand Down
2 changes: 1 addition & 1 deletion numcosmo/math/ncm_powspec.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -739,7 +739,7 @@ ncm_powspec_eval (NcmPowspec *powspec, NcmModel *model, const gdouble z, const g
* @model: a #NcmModel
* @z: time $z$
* @k: a #NcmVector
* @Pk: (out caller-allocates): a #NcmVector
* @Pk: a #NcmVector
*
* Evaluates the power spectrum @powspec at $z$ and in the knots
* contained in @k and puts the result in @Pk.
Expand Down
142 changes: 97 additions & 45 deletions numcosmo/nc_powspec_mnl_halofit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -427,14 +427,35 @@ _nc_powspec_mnl_halofit_linear_scale (NcPowspecMNLHaloFit *pshf, NcHICosmo *cosm
FDF.params = &vps;

/* Check if f(lnR_min) f(lnR_max) are both positive/negative. */
if (_nc_powspec_mnl_halofit_varm1 (lnR_min, &vps) * _nc_powspec_mnl_halofit_varm1 (lnR_max, &vps) > 0.0)
g_error ("_nc_powspec_mnl_halofit_linear_scale: cannot find linear scale. "
"The lnR-interval (% 22.15g, % 22.15g) with lnVar (% 22.15g, % 22.15g) "
"does not seem to include lnVar == 0.0. This interval can be increased by increasing "
"the lnk-interval for the linear power spectrum.",
lnR_min, lnR_max,
_nc_powspec_mnl_halofit_varm1 (lnR_min, &vps),
_nc_powspec_mnl_halofit_varm1 (lnR_max, &vps));
{
const gdouble varm1_min = _nc_powspec_mnl_halofit_varm1 (lnR_min, &vps);
const gdouble varm1_max = _nc_powspec_mnl_halofit_varm1 (lnR_max, &vps);

if (varm1_min * varm1_max > 0.0)
{
if ((varm1_min < 0.0) && (varm1_max < 0.0))
{
const gdouble Rnlin = 0.0;

/*
* All considered scales are linear, variance==1 R is below the lower
* limit, returning 0,0
*/

return Rnlin;
}
else
{
g_error ("_nc_powspec_mnl_halofit_linear_scale: cannot find linear scale. "
"The lnR-interval (% 22.15g, % 22.15g) with lnVar (% 22.15g, % 22.15g) "
"does not seem to include lnVar == 0.0. This interval can be increased by increasing "
"the lnk-interval for the linear power spectrum.",
lnR_min, lnR_max,
_nc_powspec_mnl_halofit_varm1 (lnR_min, &vps),
_nc_powspec_mnl_halofit_varm1 (lnR_max, &vps));
}
}
}

gsl_root_fdfsolver_set (self->linear_scale_solver, &FDF, lnR);

Expand Down Expand Up @@ -555,8 +576,15 @@ _nc_powspec_mnl_halofit_prepare_nl (NcPowspecMNLHaloFit *pshf, NcmModel *model)
gdouble R1 = _nc_powspec_mnl_halofit_linear_scale_z (z1, Fznl.params);
gboolean found = FALSE;

/* If the model is linear for all chosen scales, the linear Pk is enough. */
if (R0 < R_min)
g_error ("_nc_powspec_mnl_halofit_prepare_nl: linear universe or too large R_min, in the latter case increase k_max (R0 == % 21.15g, R_min == % 21.15g).", R0, R_min);
{
self->znl = -1.0;
ncm_spline_clear (&self->Rsigma);
self->Rsigma = ncm_spline_cubic_notaknot_new ();

return;
}

while (R1 > R_min)
{
Expand All @@ -580,13 +608,19 @@ _nc_powspec_mnl_halofit_prepare_nl (NcPowspecMNLHaloFit *pshf, NcmModel *model)
gsl_root_fsolver_set (self->znl_solver, &Fznl, z0, z1);

do {
gdouble znl_try;

iter++;
status = gsl_root_fsolver_iterate (self->znl_solver);

self->znl = gsl_root_fsolver_root (self->znl_solver);
z0 = gsl_root_fsolver_x_lower (self->znl_solver);
z1 = gsl_root_fsolver_x_upper (self->znl_solver);
status = gsl_root_test_interval (z0, z1, 0.0, 1.0e-3);
znl_try = gsl_root_fsolver_root (self->znl_solver);
z0 = gsl_root_fsolver_x_lower (self->znl_solver);
z1 = gsl_root_fsolver_x_upper (self->znl_solver);
status = gsl_root_test_interval (z0, z1, 0.0, 1.0e-3);

/* Due to the non-continuous nature of Fznl we need to check if we are not in a bump. */
if (fabs (Fznl.function (znl_try, Fznl.params)) < fabs (Fznl.function (self->znl, Fznl.params)))
self->znl = znl_try;
} while (status == GSL_CONTINUE && iter < max_iter);

if (iter >= max_iter)
Expand Down Expand Up @@ -812,21 +846,29 @@ _nc_powspec_mnl_halofit_eval (NcmPowspec *powspec, NcmModel *model, const gdoubl
NcHICosmo *cosmo = NC_HICOSMO (model);
NcPowspecMNLHaloFit *pshf = NC_POWSPEC_MNL_HALOFIT (powspec);
NcPowspecMNLHaloFitPrivate * const self = pshf->priv;
const gboolean linscale = (z > self->znl);
const gboolean applysmooth = (z + 1.0 > self->znl);
const gdouble zhf = linscale ? self->znl : z;
const gdouble Pklin = ncm_powspec_eval (NCM_POWSPEC (self->psml), model, z, k);
gdouble Pknln;

if (zhf != self->z)
_nc_powspec_mnl_halofit_preeval (pshf, cosmo, zhf);
if (self->znl < 0.0)
{
return ncm_powspec_eval (NCM_POWSPEC (self->psml), model, z, k);
}
else
{
const gboolean linscale = (z > self->znl);
const gboolean applysmooth = (z + 1.0 > self->znl);
const gdouble zhf = linscale ? self->znl : z;
const gdouble Pklin = ncm_powspec_eval (NCM_POWSPEC (self->psml), model, z, k);
gdouble Pknln;

if (zhf != self->z)
_nc_powspec_mnl_halofit_preeval (pshf, cosmo, zhf);

Pknln = _nc_powspec_mnl_halofit_Pklin2Pknln (pshf, cosmo, k, Pklin);
Pknln = _nc_powspec_mnl_halofit_Pklin2Pknln (pshf, cosmo, k, Pklin);

if (applysmooth)
Pknln = ncm_util_smooth_trans (Pknln, Pklin, self->znl, 1.0, z);
if (applysmooth)
Pknln = ncm_util_smooth_trans (Pknln, Pklin, self->znl, 1.0, z);

return Pknln;
return Pknln;
}
}

static void
Expand All @@ -835,34 +877,42 @@ _nc_powspec_mnl_halofit_eval_vec (NcmPowspec *powspec, NcmModel *model, const gd
NcHICosmo *cosmo = NC_HICOSMO (model);
NcPowspecMNLHaloFit *pshf = NC_POWSPEC_MNL_HALOFIT (powspec);
NcPowspecMNLHaloFitPrivate * const self = pshf->priv;
const gboolean linscale = (z > self->znl);
const gboolean applysmooth = (z + 1.0 > self->znl);
const gdouble zhf = linscale ? self->znl : z;
gdouble theta0, theta1;

ncm_powspec_eval_vec (NCM_POWSPEC (self->psml), model, z, k, Pk);
if (self->znl < 0.0)
{
ncm_powspec_eval_vec (NCM_POWSPEC (self->psml), model, z, k, Pk);
}
else
{
const gboolean linscale = (z > self->znl);
const gboolean applysmooth = (z + 1.0 > self->znl);
const gdouble zhf = linscale ? self->znl : z;
gdouble theta0, theta1;

if (applysmooth)
ncm_util_smooth_trans_get_theta (self->znl, 1.0, z, &theta0, &theta1);
ncm_powspec_eval_vec (NCM_POWSPEC (self->psml), model, z, k, Pk);

if (zhf != self->z)
_nc_powspec_mnl_halofit_preeval (pshf, cosmo, zhf);
if (applysmooth)
ncm_util_smooth_trans_get_theta (self->znl, 1.0, z, &theta0, &theta1);

{
const guint len = ncm_vector_len (k);
guint i;
if (zhf != self->z)
_nc_powspec_mnl_halofit_preeval (pshf, cosmo, zhf);

for (i = 0; i < len; i++)
{
const gdouble ki = ncm_vector_get (k, i);
const gdouble Pklin = ncm_vector_get (Pk, i);
const guint len = ncm_vector_len (k);
guint i;

for (i = 0; i < len; i++)
{
const gdouble ki = ncm_vector_get (k, i);
const gdouble Pklin = ncm_vector_get (Pk, i);

const gdouble Pknln = _nc_powspec_mnl_halofit_Pklin2Pknln (pshf, cosmo, ki, Pklin);
const gdouble Pknln = _nc_powspec_mnl_halofit_Pklin2Pknln (pshf, cosmo, ki, Pklin);

if (applysmooth)
ncm_vector_set (Pk, i, theta0 * Pknln + theta1 * Pklin);
else
ncm_vector_set (Pk, i, Pknln);
if (applysmooth)
ncm_vector_set (Pk, i, theta0 * Pknln + theta1 * Pklin);
else
ncm_vector_set (Pk, i, Pknln);
}
}
}
}
Expand All @@ -874,7 +924,9 @@ _nc_powspec_mnl_halofit_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk)
NcPowspecMNLHaloFitPrivate * const self = pshf->priv;

ncm_powspec_get_nknots (NCM_POWSPEC (self->psml), Nz, Nk);
*Nz = ncm_vector_len (ncm_spline_peek_xv (self->Rsigma));

if (self->znl > 0.0)
*Nz = ncm_vector_len (ncm_spline_peek_xv (self->Rsigma));
}

/**
Expand Down
Loading

0 comments on commit 31a590e

Please sign in to comment.