Skip to content

Commit

Permalink
Merge d83fc3a into 8ae4c81
Browse files Browse the repository at this point in the history
  • Loading branch information
@pennalima
pennalima committed Sep 28, 2021
2 parents 8ae4c81 + d83fc3a commit 9ebb56c
Show file tree
Hide file tree
Showing 4 changed files with 138 additions and 78 deletions.
79 changes: 51 additions & 28 deletions numcosmo/lss/nc_halo_density_profile.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -350,7 +350,7 @@ nc_halo_density_profile_class_init (NcHaloDensityProfileClass *klass)
*
*/
ncm_model_class_set_sparam (model_class, NC_HALO_DENSITY_PROFILE_C_DELTA, "c_{\\Delta}", "cDelta",
0.5, 10.0, 1.0e-1,
2.5, 10.0, 1.0e-1,
NC_HALO_DENSITY_PROFILE_DEFAULT_PARAMS_ABSTOL, NC_HALO_DENSITY_PROFILE_DEFAULT_C_DELTA,
NCM_PARAM_TYPE_FIXED);

Expand Down Expand Up @@ -505,11 +505,12 @@ _nc_halo_density_profile_prepare_dl_spher_mass (NcHaloDensityProfile *dp)
gsl_integration_workspace **w = ncm_integral_get_workspace ();
gsl_function F;
gdouble err;
gint key = 6;

F.function = &_nc_halo_density_profile_prepare_dl_spher_mass_int;
F.params = dp;

gsl_integration_qag (&F, 0.0, C_DELTA, 0.0, NCM_DEFAULT_PRECISION, NCM_INTEGRAL_PARTITION, 6, *w, &self->dl_spher_mass, &err);
gsl_integration_qag (&F, 0.0, C_DELTA, 0.0, self->reltol, NCM_INTEGRAL_PARTITION, key, *w, &self->dl_spher_mass, &err);

ncm_memory_pool_return (w);
ncm_model_lstate_set_update (NCM_MODEL (dp), PREPARE_DL_SPHER_MASS);
Expand All @@ -526,23 +527,40 @@ typedef struct _NcHaloDensityProfile2D
} NcHaloDensityProfile2D;

static gdouble
_nc_halo_density_profile_prepare_dl_2d_density_X_u (gdouble u, gpointer userdata)
_nc_halo_density_profile_prepare_dl_2d_density_X_u (gdouble lnu, gpointer userdata)
{
NcHaloDensityProfile2D *dp2D = (NcHaloDensityProfile2D *) userdata;
const gdouble u = exp (lnu);

return nc_halo_density_profile_eval_dl_density (dp2D->dp, hypot (dp2D->X, u));
return nc_halo_density_profile_eval_dl_density (dp2D->dp, hypot (dp2D->X, u)) * u;
}

static gdouble
_nc_halo_density_profile_prepare_dl_2d_density_X (gdouble lnX, gpointer userdata)
{
NcHaloDensityProfile2D *dp2D = (NcHaloDensityProfile2D *) userdata;
NcHaloDensityProfilePrivate * const self = dp2D->dp->priv;
gdouble err, dl_2d_density_X;
gdouble err, dl_2d_density_dX, dl_2d_density_X, lnx0, lnx1;
const gdouble dlnx = 2.0 * M_LN10;
gdouble abstol = 0.0;
gint key = 6;

dp2D->X = exp (lnX);
gsl_integration_qagiu (dp2D->F, 0.0, abstol, self->reltol, NCM_INTEGRAL_PARTITION, dp2D->w, &dl_2d_density_X, &err);

dl_2d_density_X = 0.0;
lnx0 = lnX + GSL_LOG_DBL_EPSILON;
lnx1 = lnx0 + dlnx;
do
{
gsl_integration_qag (dp2D->F, lnx0, lnx1, abstol, self->reltol, NCM_INTEGRAL_PARTITION, key, dp2D->w, &dl_2d_density_dX, &err);

dl_2d_density_X += dl_2d_density_dX;
abstol = dl_2d_density_X * self->reltol;

lnx0 = lnx1;
lnx1 = lnx0 + dlnx;

} while (fabs (dl_2d_density_dX / dl_2d_density_X) > self->reltol);

return log (2.0 * dl_2d_density_X);
}
Expand Down Expand Up @@ -576,32 +594,23 @@ _nc_halo_density_profile_prepare_dl_2d_density (NcHaloDensityProfile *dp)
}

static gdouble
_nc_halo_density_profile_prepare_dl_cyl_mass_X_x_1 (gdouble x, gpointer userdata)
_nc_halo_density_profile_prepare_dl_cyl_mass_X_x_1 (gdouble lnx, gpointer userdata)
{
NcHaloDensityProfile2D *dp2D = (NcHaloDensityProfile2D *) userdata;

if (x == 0.0)
{
return 0.0;
}
else
{
const gdouble x2 = x * x;

return x2 * nc_halo_density_profile_eval_dl_density (dp2D->dp, x);
}
const gdouble x = exp (lnx);
const gdouble x2 = x * x;

return x2 * x * nc_halo_density_profile_eval_dl_density (dp2D->dp, x);
}

static gdouble
_nc_halo_density_profile_prepare_dl_cyl_mass_X_x_2 (gdouble mu, gpointer userdata)
_nc_halo_density_profile_prepare_dl_cyl_mass_X_x_2 (gdouble lnmu, gpointer userdata)
{
NcHaloDensityProfile2D *dp2D = (NcHaloDensityProfile2D *) userdata;
const gdouble mu2 = mu * mu;
const gdouble fa = sqrt (fabs (expm1 (-2.0 * lnmu)));
const gdouble mu = exp (lnmu);

if (mu == 0.0)
return 0.0;
else
return nc_halo_density_profile_eval_dl_density (dp2D->dp, dp2D->X / mu) / (mu2 * (1.0 + sqrt ((1.0 - mu) * (1.0 + mu))));
return nc_halo_density_profile_eval_dl_density (dp2D->dp, dp2D->X * mu) * mu / (1.0 + fa);
}

static gdouble
Expand All @@ -611,16 +620,30 @@ _nc_halo_density_profile_prepare_dl_cyl_mass_X (gdouble lnX, gpointer userdata)
NcHaloDensityProfilePrivate * const self = dp2D->dp->priv;

gdouble dl_cyl_mass_X = 0.0;
gdouble err, dl_cyl_mass_X_i;
gdouble err, dl_cyl_mass_X_i, X3;
gdouble abstol = 0.0;
gint key = 6;

dp2D->X = exp (lnX);
X3 = gsl_pow_3 (dp2D->X);

gsl_integration_qag (dp2D->F, 0.0, dp2D->X, abstol, self->reltol, NCM_INTEGRAL_PARTITION, 6, dp2D->w, &dl_cyl_mass_X_i, &err);
gsl_integration_qag (dp2D->F, 2.0 * GSL_LOG_DBL_EPSILON, lnX, abstol, self->reltol, NCM_INTEGRAL_PARTITION, key, dp2D->w, &dl_cyl_mass_X_i, &err);
dl_cyl_mass_X += dl_cyl_mass_X_i;

gsl_integration_qag (dp2D->F2, 0.0, 1.0, abstol, self->reltol, NCM_INTEGRAL_PARTITION, key, dp2D->w, &dl_cyl_mass_X_i, &err);
dl_cyl_mass_X += X3 * dl_cyl_mass_X_i;

gsl_integration_qag (dp2D->F2, 0.0, 1.0, abstol, self->reltol, NCM_INTEGRAL_PARTITION, 6, dp2D->w, &dl_cyl_mass_X_i, &err);
dl_cyl_mass_X += gsl_pow_3 (dp2D->X) * dl_cyl_mass_X_i;
abstol += dl_cyl_mass_X_i * self->reltol;

gsl_integration_qag (dp2D->F2, 1.0, 10.0, abstol, self->reltol, NCM_INTEGRAL_PARTITION, key, dp2D->w, &dl_cyl_mass_X_i, &err);
dl_cyl_mass_X += X3 * dl_cyl_mass_X_i;

abstol += dl_cyl_mass_X_i * self->reltol;

gsl_integration_qag (dp2D->F2, 10.0, GSL_LOG_DBL_MAX / 4.0, abstol, self->reltol, NCM_INTEGRAL_PARTITION, key, dp2D->w, &dl_cyl_mass_X_i, &err);
dl_cyl_mass_X += X3 * dl_cyl_mass_X_i;

//printf ("% 22.15e\n", dl_cyl_mass_X_i);

return log (2.0 * dl_cyl_mass_X);
}
Expand Down
7 changes: 4 additions & 3 deletions numcosmo/lss/nc_halo_density_profile_einasto.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,8 @@
* \end{equation}
* where $c_{\Delta}$ is the concentration parameter.
*
* References: Einasto (1965), arXiv:1712.04512.
* References: Einasto (1965), arXiv:1712.04512. $\alpha$ parameter range values: Gao et al. (2008) https://ui.adsabs.harvard.edu/abs/2008MNRAS.387..536G/abstract,
* Dutton \& Maccio (2014) https://academic.oup.com/mnras/article/441/4/3359/1209689.
*/

#ifdef HAVE_CONFIG_H
Expand Down Expand Up @@ -143,7 +144,7 @@ nc_halo_density_profile_einasto_class_init (NcHaloDensityProfileEinastoClass *kl
* FIXME Set correct values (limits)
*/
ncm_model_class_set_sparam (model_class, NC_HALO_DENSITY_PROFILE_EINASTO_ALPHA, "\\alpha", "alpha",
0.155, 0.5, 0.01,
0.12, 0.35, 0.01,
NC_HALO_DENSITY_PROFILE_EINASTO_DEFAULT_PARAMS_ABSTOL, NC_HALO_DENSITY_PROFILE_EINASTO_DEFAULT_ALPHA,
NCM_PARAM_TYPE_FIXED);

Expand All @@ -166,7 +167,7 @@ _nc_halo_density_profile_einasto_eval_dl_spher_mass (NcHaloDensityProfile *dp, c
const gdouble gamma_3_alpha = gsl_sf_gamma (3.0 / ALPHA);
const gdouble arg_2 = 2.0 * pow (C_DELTA, ALPHA) / ALPHA;
const gdouble gamma_inc_P = gsl_sf_gamma_inc_P (3.0 / ALPHA, arg_2);

return (pow (ALPHA / 2.0, 3.0 / ALPHA) * exp (2.0 / ALPHA) / ALPHA * gamma_3_alpha * gamma_inc_P);
}

Expand Down
122 changes: 79 additions & 43 deletions tests/test_nc_halo_density_profile.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,14 +32,17 @@
#include <glib.h>
#include <glib-object.h>

typedef struct _TestNcHaloDensityProfile
typedef struct _TestNcHaloDensityProfile TestNcHaloDensityProfile;

struct _TestNcHaloDensityProfile
{
NcHaloDensityProfile *dp;
NcHICosmo *cosmo;
gdouble z;
gdouble R1, R2, R3;
gint ntests;
} TestNcHaloDensityProfile;
gboolean (*rng_params) (TestNcHaloDensityProfile *test);
};

void test_nc_halo_density_profile_nfw_new (TestNcHaloDensityProfile *test, gconstpointer pdata);
void test_nc_halo_density_profile_nfw_eval_density (TestNcHaloDensityProfile *test, gconstpointer pdata);
Expand Down Expand Up @@ -142,13 +145,14 @@ test_nc_halo_density_profile_nfw_new (TestNcHaloDensityProfile *test, gconstpoin
ncm_model_param_set_by_name (NCM_MODEL (dp), "log10MDelta", 15.0);
ncm_model_param_set_by_name (NCM_MODEL (dp), "cDelta", 4.0);

test->cosmo = cosmo;
test->dp = dp;
test->z = 1.0;
test->R1 = 0.3; /* Mpc */
test->R2 = nc_halo_density_profile_r_s (dp, cosmo, test->z);
test->R3 = 10.0;
test->ntests = 1000;
test->cosmo = cosmo;
test->dp = dp;
test->z = 1.0;
test->R1 = 0.3; /* Mpc */
test->R2 = nc_halo_density_profile_r_s (dp, cosmo, test->z);
test->R3 = 10.0;
test->ntests = 1000;
test->rng_params = NULL;

nc_distance_free (dist);
}
Expand Down Expand Up @@ -177,17 +181,35 @@ test_nc_halo_density_profile_hernquist_new (TestNcHaloDensityProfile *test, gcon
ncm_model_param_set_by_name (NCM_MODEL (dp), "log10MDelta", 15.0);
ncm_model_param_set_by_name (NCM_MODEL (dp), "cDelta", 4.0);

test->cosmo = cosmo;
test->dp = dp;
test->z = 1.0;
test->R1 = 0.3; /* Mpc */
test->R2 = nc_halo_density_profile_r_s (dp, cosmo, test->z);
test->R3 = 10.0;
test->ntests = 1000;
test->cosmo = cosmo;
test->dp = dp;
test->z = 1.0;
test->R1 = 0.3; /* Mpc */
test->R2 = nc_halo_density_profile_r_s (dp, cosmo, test->z);
test->R3 = 10.0;
test->ntests = 1000;
test->rng_params = NULL;

nc_distance_free (dist);
}

static gboolean
_test_nc_halo_density_profile_einasto_rng (TestNcHaloDensityProfile *test)
{
NcmModel *model = NCM_MODEL (test->dp);
const gdouble log10MDelta = g_test_rand_double_range (ncm_model_param_get_lower_bound (model, NC_HALO_DENSITY_PROFILE_LOG10M_DELTA), ncm_model_param_get_upper_bound (model, NC_HALO_DENSITY_PROFILE_LOG10M_DELTA));
const gdouble cDelta = g_test_rand_double_range (ncm_model_param_get_lower_bound (model, NC_HALO_DENSITY_PROFILE_C_DELTA), ncm_model_param_get_upper_bound (model, NC_HALO_DENSITY_PROFILE_C_DELTA));
const gdouble alpha = g_test_rand_double_range (ncm_model_param_get_lower_bound (model, NC_HALO_DENSITY_PROFILE_EINASTO_ALPHA), ncm_model_param_get_upper_bound (model, NC_HALO_DENSITY_PROFILE_EINASTO_ALPHA));

ncm_model_param_set_by_name (model, "log10MDelta", log10MDelta);
ncm_model_param_set_by_name (model, "cDelta", cDelta);
ncm_model_param_set_by_name (model, "alpha", alpha);

/*printf ("% 22.15g % 22.15g % 22.15g\n", log10MDelta, cDelta, alpha);*/

return TRUE;
}

void
test_nc_halo_density_profile_einasto_new (TestNcHaloDensityProfile *test, gconstpointer pdata)
{
Expand All @@ -211,14 +233,16 @@ test_nc_halo_density_profile_einasto_new (TestNcHaloDensityProfile *test, gconst

ncm_model_param_set_by_name (NCM_MODEL (dp), "log10MDelta", 15.0);
ncm_model_param_set_by_name (NCM_MODEL (dp), "cDelta", 4.0);
ncm_model_param_set_by_name (NCM_MODEL (dp), "alpha", 0.25);

test->cosmo = cosmo;
test->dp = dp;
test->z = 1.0;
test->R1 = 0.3; /* Mpc */
test->R2 = nc_halo_density_profile_r_s (dp, cosmo, test->z);
test->R3 = 10.0;
test->ntests = 1000;
test->cosmo = cosmo;
test->dp = dp;
test->z = 1.0;
test->R1 = 0.3; /* Mpc */
test->R2 = nc_halo_density_profile_r_s (dp, cosmo, test->z);
test->R3 = 10.0;
test->ntests = 1000;
test->rng_params = &_test_nc_halo_density_profile_einasto_rng;

nc_distance_free (dist);
}
Expand Down Expand Up @@ -253,47 +277,59 @@ void
test_nc_halo_density_profile_eval_dl_spher_mass (TestNcHaloDensityProfile *test, gconstpointer pdata)
{
NcHaloDensityProfile *dp = test->dp;
gint np = 10;
gint i;
for (i = 0; i < test->ntests; i++)

do
{
const gdouble X = pow (10.0, g_test_rand_double_range (-3.0, 3.0));
const gdouble ISigma = nc_halo_density_profile_eval_dl_spher_mass (dp, X);
const gdouble NISigma = nc_halo_density_profile_eval_numint_dl_spher_mass (dp, X);

for (i = 0; i < test->ntests; i++)
{
const gdouble X = pow (10.0, g_test_rand_double_range (-3.0, 3.0));
const gdouble ISigma = nc_halo_density_profile_eval_dl_spher_mass (dp, X);
const gdouble NISigma = nc_halo_density_profile_eval_numint_dl_spher_mass (dp, X);

ncm_assert_cmpdouble_e (ISigma, ==, NISigma, nc_halo_density_profile_get_reltol (dp) * 1.0e1, 0.0);
}
ncm_assert_cmpdouble_e (ISigma, ==, NISigma, nc_halo_density_profile_get_reltol (dp) * 1.0e1, 0.0);
}
} while ((test->rng_params != NULL) && test->rng_params (test) && np--);
}

void
test_nc_halo_density_profile_eval_dl_2d_density (TestNcHaloDensityProfile *test, gconstpointer pdata)
{
NcHaloDensityProfile *dp = test->dp;
gint np = 10;
gint i;

for (i = 0; i < test->ntests; i++)
do
{
const gdouble X = pow (10.0, g_test_rand_double_range (-3.0, 3.0));
const gdouble ISigma = nc_halo_density_profile_eval_dl_2d_density (dp, X);
const gdouble NISigma = nc_halo_density_profile_eval_numint_dl_2d_density (dp, X);
for (i = 0; i < test->ntests; i++)
{
const gdouble X = pow (10.0, g_test_rand_double_range (-3.0, 3.0));
const gdouble ISigma = nc_halo_density_profile_eval_dl_2d_density (dp, X);
const gdouble NISigma = nc_halo_density_profile_eval_numint_dl_2d_density (dp, X);

ncm_assert_cmpdouble_e (ISigma, ==, NISigma, nc_halo_density_profile_get_reltol (dp) * 1.0e1, 0.0);
}
ncm_assert_cmpdouble_e (ISigma, ==, NISigma, nc_halo_density_profile_get_reltol (dp) * 1.0e1, 0.0);
}
} while ((test->rng_params != NULL) && test->rng_params (test) && np--);
}

void
test_nc_halo_density_profile_eval_dl_cyl_mass (TestNcHaloDensityProfile *test, gconstpointer pdata)
{
NcHaloDensityProfile *dp = test->dp;
gint np = 10;
gint i;

for (i = 0; i < test->ntests; i++)
do
{
const gdouble X = pow (10.0, g_test_rand_double_range (-3.0, 3.0));
const gdouble ISigma = nc_halo_density_profile_eval_dl_cyl_mass (dp, X);
const gdouble NISigma = nc_halo_density_profile_eval_numint_dl_cyl_mass (dp, X);
for (i = 0; i < test->ntests; i++)
{
const gdouble X = pow (10.0, g_test_rand_double_range (-3.0, 3.0));
const gdouble ISigma = nc_halo_density_profile_eval_dl_cyl_mass (dp, X);
const gdouble NISigma = nc_halo_density_profile_eval_numint_dl_cyl_mass (dp, X);

ncm_assert_cmpdouble_e (ISigma, ==, NISigma, nc_halo_density_profile_get_reltol (dp) * 1.0e1, 0.0);
}
ncm_assert_cmpdouble_e (ISigma, ==, NISigma, nc_halo_density_profile_get_reltol (dp) * 1.0e1, 0.0);
}
} while ((test->rng_params != NULL) && test->rng_params (test) && np--);
}

8 changes: 4 additions & 4 deletions tests/test_ncm_fit_esmcmc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -377,7 +377,7 @@ test_ncm_fit_esmcmc_run (TestNcmFitESMCMC *test, gconstpointer pdata)
/*fprintf (stderr, "retrying[%2d] %d %d\n", max_tries, run, (gint)(run * 2));*/
run = run * 2;

/*ncm_fit_esmcmc_set_mtype (test->esmcmc, NCM_FIT_RUN_MSGS_SIMPLE);*/
/*ncm_fit_esmcmc_set_mtype (test->esmcmc, NCM_FIT_RUN_MSGS_NONE);*/

ncm_fit_esmcmc_start_run (test->esmcmc);
ncm_fit_esmcmc_run (test->esmcmc, run);
Expand All @@ -393,7 +393,7 @@ test_ncm_fit_esmcmc_run (TestNcmFitESMCMC *test, gconstpointer pdata)
void
test_ncm_fit_esmcmc_run_lre (TestNcmFitESMCMC *test, gconstpointer pdata)
{
const gint run = test->dim * g_test_rand_int_range (1500, 2000) / test->nrun_div;
const gint run = MAX (test->dim * g_test_rand_int_range (1500, 2000) / test->nrun_div, 100);
NcmMatrix *data_cov = ncm_matrix_dup (NCM_DATA_GAUSS_COV (test->data_mvnd)->cov);
gboolean ok = FALSE;
gint max_tries = 15;
Expand Down Expand Up @@ -592,7 +592,7 @@ test_ncm_fit_esmcmc_run_restart_from_cat (TestNcmFitESMCMC *test, gconstpointer
void
test_ncm_fit_esmcmc_run_lre_auto_trim (TestNcmFitESMCMC *test, gconstpointer pdata)
{
const gint run = test->dim * g_test_rand_int_range (6000, 12000) / test->nrun_div;
const gint run = MAX (test->dim * g_test_rand_int_range (6000, 12000) / test->nrun_div, 100);
gdouble prec = 1.0e-2;
NcmMatrix *data_cov = ncm_matrix_dup (NCM_DATA_GAUSS_COV (test->data_mvnd)->cov);
NcmMatrix *data_cor = ncm_matrix_dup (NCM_DATA_GAUSS_COV (test->data_mvnd)->cov);
Expand Down Expand Up @@ -660,7 +660,7 @@ test_ncm_fit_esmcmc_run_lre_auto_trim (TestNcmFitESMCMC *test, gconstpointer pda
void
test_ncm_fit_esmcmc_run_lre_auto_trim_vol (TestNcmFitESMCMC *test, gconstpointer pdata)
{
const gint run = test->dim * g_test_rand_int_range (6000, 12000) / test->nrun_div;
const gint run = MAX (test->dim * g_test_rand_int_range (6000, 12000) / test->nrun_div, 100);
gdouble prec = 1.0e-2;
gdouble lnnorm_sd, lnnorm;

Expand Down

0 comments on commit 9ebb56c

Please sign in to comment.