Skip to content

Commit

Permalink
New minor version. Several improvements.
Browse files Browse the repository at this point in the history 
Better building tool support for lapack/BLAS.
Improved speed using different interfaces for Lapack/BLAS.
Included levmar internally and removed as optional dependency.
Levmar intreface now supports box constraints.
Added virtual functions for lnNorma2 and resample in NcmDataGaussCov.
Updated interval vector usage in NcmFitNLOpt and NcmFitLevmar.
Added support for estimating true width and colour in NcDataSNIACov.
Improved resampling in NcDataSNIACov, now it resamples all data m_B, w and c.
Several new functions in NcmMatrix (support to writing data in col-major order).
New ncm_mset_catalog_calc_ci for calculating confidence intervals for generic function using a NcmMSetCatalog.
Added option for confidence intervals from catalogs in darkenergy.
Added templates as dependency for NcmFitNLOpt enums files.

NcDataSNIACov now uses the (hopefully) the right normalization (new feature still experimental).

Added some misc statistical functions in NcmUtil (experimental).
  • Loading branch information
Sandro Dias Pinto Vitenti committed Mar 1, 2015
1 parent d173055 commit 0545e3c
Show file tree
Hide file tree
Showing 63 changed files with 10,742 additions and 917 deletions.
95 changes: 40 additions & 55 deletions configure.ac
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@ dnl ***************************************************************************

m4_define([numcosmo_major_version], [0])
m4_define([numcosmo_minor_version], [12])
m4_define([numcosmo_micro_version], [0])
m4_define([numcosmo_interface_age], [0])
m4_define([numcosmo_micro_version], [1])
m4_define([numcosmo_interface_age], [1])
m4_define([numcosmo_binary_age],
[m4_eval(100 * numcosmo_minor_version + numcosmo_micro_version)])
m4_define([numcosmo_version],
Expand Down Expand Up @@ -442,12 +442,17 @@ if test "x$have_blas" != "xno"; then
fi
else
AC_MSG_RESULT([GSL's BLAS])
AC_DEFINE([HAVE_GSL_CBLAS_H], [1], [use lapack functions])
AC_DEFINE([HAVE_GSL_CBLAS_H], [1], [use GSL's BLAS])
fi

dnl ***************************************************************************
dnl Check for LAPACK
dnl ***************************************************************************
AC_ARG_WITH(lapacke,
[AS_HELP_STRING([--with-lapacke=<lib>], [use LAPACKe library <lib>])])
case $with_lapacke in
yes | "") ;;
no) lapacke_ok=disable ;;
-* | */* | *.a | *.so | *.so.* | *.o) LAPACKE_LIBS="$with_lapacke" ;;
*) LAPACKE_LIBS="-l$with_lapacke" ;;
esac

AX_LAPACK([],[have_lapack="no"])
AC_MSG_CHECKING([if lapack is available $LAPACK_LIBS])
Expand All @@ -467,60 +472,45 @@ if test "x$have_lapack" != "xno"; then
else
lapack_support="Other"
fi

LIBS="$LAPACK_LIBS $LIBS"
AC_MSG_RESULT([$lapack_support])

if test "x$lapack_mkl" != "x"; then
AC_CHECK_HEADERS([mkl_lapacke.h],
[AC_DEFINE([HAVE_LAPACKE],[1], [have lapacke support])
AC_DEFINE([HAVE_MKL_LAPACKE_H],[1], [have mkl_lapacke.h header])])
[AC_CHECK_FUNC([LAPACKE_dptsv],
[AC_DEFINE([HAVE_LAPACKE],[1], [have lapacke support])
AC_DEFINE([HAVE_MKL_LAPACKE_H],[1], [have mkl_lapacke.h header])
have_lapacke="yes"][])])
AC_CHECK_HEADERS([mkl_lapack.h])
else
AC_CHECK_HEADERS([lapacke.h],
[AC_DEFINE([HAVE_LAPACKE],[1], [have lapacke support])
AC_DEFINE([HAVE_LAPACKE_H],[1], [have lapacke.h header])])
AC_CHECK_HEADERS([clapack.h],
[AC_DEFINE([HAVE_CLAPACK],[1], [have clapack support])
AC_DEFINE([HAVE_CLAPACK_H],[1], [have clapack.h header])])
fi
LIBS="$LAPACK_LIBS $LIBS"
else
AC_MSG_RESULT([no])
fi

dnl ***************************************************************************
dnl Check for levmar
dnl ***************************************************************************

AC_MSG_CHECKING([--with-levmar-path])
AC_ARG_WITH(levmar-path,
[AS_HELP_STRING([--with-levmar-path=PATH],[Levmar library location])],
[LEVMAR_PATH="-L$withval"], [])
AC_MSG_RESULT($LEVMAR_PATH)

if test "x$LEVMAR_PATH" != x; then
LIBS="$LEVMAR_PATH $LIBS"
fi
if test "x$have_lapacke" != "xyes"; then
AC_CHECK_HEADERS([lapacke.h],
[AC_SEARCH_LIBS([LAPACKE_dptsv],[lapacke],
[AC_DEFINE([HAVE_LAPACKE],[1], [have lapacke support])
AC_DEFINE([HAVE_LAPACKE_H],[1], [have lapacke.h header])
have_lapacke="yes"],[])])
fi

AC_MSG_CHECKING([--with-levmar-extra-libs])
AC_ARG_WITH(levmar-extra-libs,
[AS_HELP_STRING([--with-levmar-extra-libs=LIBS],[Levmar extra libraries, default is none])],
[LEVMAR_EXTRA_LIBS="$withval"], [])
AC_MSG_RESULT($LEVMAR_EXTRA_LIBS)
if test "x$have_lapacke" != "xyes"; then
AC_CHECK_HEADERS([clapack.h],
[AC_SEARCH_LIBS([clapack_dpotri],[clapack],
[AC_DEFINE([HAVE_CLAPACK],[1], [have clapack support])
AC_DEFINE([HAVE_CLAPACK_H],[1], [have clapack.h header])
have_clapack="yes"],[])])
fi

if test "x$have_lapacke" = "xyes"; then
lapack_support="$lapack_support+LAPACKe"
fi

AC_CHECK_LIB([levmar], [dlevmar_der],[LIBS="-llevmar $LEVMAR_EXTRA_LIBS $LIBS"],[levmar_missing="yes"],[$LEVMAR_EXTRA_LIBS])
if test "x$have_clapack" = "xyes"; then
lapack_support="$lapack_support+cLAPACK"
fi

AC_CHECK_HEADER(levmar.h, , [AC_CHECK_HEADER(levmar/levmar.h, [levmar_need_prefix="yes"], [levmar_missing="yes"])])
if test "x$levmar_missing" = "xyes"; then
have_levmar_support=""
else
have_levmar_support="#define NUMCOSMO_HAVE_LEVMAR 1"
fi
AC_SUBST(have_levmar_support)

AM_CONDITIONAL([HAVE_LIBLEVMAR], [test "x$have_levmar_support" != x])

if test "x$levmar_need_prefix" = "xyes"; then
AC_DEFINE([NC_LEVMAR_NEED_PREFIX],[1], [Duh!])
AC_MSG_RESULT([no])
fi

dnl ***************************************************************************
Expand Down Expand Up @@ -679,6 +669,7 @@ AC_CONFIG_FILES([
numcosmo.pc
Makefile
numcosmo/Makefile
numcosmo/levmar/Makefile
tools/Makefile
tests/Makefile
docs/Makefile
Expand Down Expand Up @@ -713,12 +704,6 @@ else
echo "Building with Lapack support: ...................................NO"
echo " Requires Lapack Linear Algebra PACKage"
fi
if [ test "x$have_levmar_support" != x ]; then
echo "Building with levmar support: ...................................YES"
else
echo "Building with levmar support: ...................................NO"
echo " Requires levmar Levenberg-Marquardt least squares library"
fi
if [ test "x$have_cfitsio_support" != x ]; then
echo "Building with cfitsio support: ..................................YES"
else
Expand Down
161 changes: 123 additions & 38 deletions darkenergy/darkenergy.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
#include <math.h>
#include <glib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cdf.h>

gint
main (gint argc, gchar *argv[])
Expand All @@ -61,6 +62,7 @@ main (gint argc, gchar *argv[])
gchar *runconf_cmd_line = NULL;
gboolean is_de = FALSE;
NcmRNG *rng = ncm_rng_pool_get ("darkenergy");
NcmMSetCatalog *mcat = NULL;

ncm_cfg_init ();

Expand Down Expand Up @@ -494,21 +496,13 @@ main (gint argc, gchar *argv[])
}
}


/* All initializations done! */

if (de_fit.resample)
{
NcmRNG *resample_rng;
if (de_fit.mc_seed > -1)
resample_rng = ncm_rng_seeded_new (NULL, de_fit.mc_seed);
else
resample_rng = ncm_rng_new (NULL);

ncm_cfg_msg_sepa ();
ncm_message ("# Resampling from fiducial model.\n");
ncm_dataset_resample (dset, fiduc, resample_rng);
ncm_rng_free (resample_rng);
ncm_dataset_resample (dset, fiduc, rng);
}

de_fit.fisher = (de_fit.fisher || (de_fit.nsigma_fisher != -1) || (de_fit.nsigma != -1) || (de_fit.onedim_cr != NULL));
Expand Down Expand Up @@ -600,6 +594,8 @@ main (gint argc, gchar *argv[])
gdouble p_value = ncm_mset_catalog_param_pdf_pvalue (mc->mcat, m2lnL, FALSE);
ncm_message ("# - pvalue for fitted model [% 20.15g] %04.2f%%.\n#\n", m2lnL, 100.0 * p_value);
}
ncm_mset_catalog_clear (&mcat);
mcat = ncm_fit_mc_get_catalog (mc);
ncm_fit_mc_clear (&mc);
}

Expand Down Expand Up @@ -632,6 +628,9 @@ main (gint argc, gchar *argv[])
gdouble p_value = ncm_mset_catalog_param_pdf_pvalue (mcbs->mcat, m2lnL, FALSE);
ncm_message ("# - pvalue for fitted model [% 20.15g] %04.2f%%.\n#\n", m2lnL, 100.0 * p_value);
}

ncm_mset_catalog_clear (&mcat);
mcat = ncm_fit_mcbs_get_catalog (mcbs);
ncm_fit_mcbs_clear (&mcbs);
}

Expand Down Expand Up @@ -669,6 +668,9 @@ main (gint argc, gchar *argv[])
ncm_fit_mcmc_mean_covar (mcmc);
ncm_mset_catalog_param_pdf (mcmc->mcat, 0);
ncm_fit_log_covar (fit);

ncm_mset_catalog_clear (&mcat);
mcat = ncm_fit_mcmc_get_catalog (mcmc);
ncm_fit_mcmc_clear (&mcmc);
}

Expand Down Expand Up @@ -716,6 +718,10 @@ main (gint argc, gchar *argv[])
ncm_fit_esmcmc_mean_covar (esmcmc);
ncm_mset_catalog_param_pdf (esmcmc->mcat, 0);
ncm_fit_log_covar (fit);

ncm_mset_catalog_clear (&mcat);
mcat = ncm_fit_esmcmc_get_catalog (esmcmc);

ncm_fit_esmcmc_clear (&esmcmc);
}

Expand Down Expand Up @@ -918,35 +924,113 @@ main (gint argc, gchar *argv[])

if (de_fit.kinematics_sigma)
{
NcmMSetFunc *dec_param_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_q);
NcmMSetFunc *E2_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_E2);
NcmMSetFunc *Em2_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_Em2);
NcmMSetFunc *dec_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_dec);
NcmMSetFunc *wec_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_wec);
gint i;
if (mcat == NULL)
{
NcmMSetFunc *dec_param_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_q);
NcmMSetFunc *E2_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_E2);
NcmMSetFunc *Em2_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_Em2);
NcmMSetFunc *dec_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_dec);
NcmMSetFunc *wec_func = nc_hicosmo_create_mset_func1 (nc_hicosmo_wec);
struct {
gchar *name;
NcmMSetFunc *func;
} kinematics[5] = {
{"H^2/H_0^2", E2_func},
{"q(z)", dec_param_func},
{"H_0^2/H^2", Em2_func},
{"DEC", dec_func},
{"WEC", wec_func}
};
gint i, j;

ncm_message ("# Kinematics data propagating the covariance matrix:\n");
for (j = 0; j < 5; j++)
{
ncm_message ("# Kinematics data propagating the covariance matrix: %s:\n", kinematics[j].name);
for (i = 0; i < de_fit.kinematics_n; i++)
{
gdouble z = de_fit.kinematics_z / (de_fit.kinematics_n - 1.0) * i;
gdouble kf_z, sigma_kf_z;
ncm_fit_function_error (fit, kinematics[j].func, &z, FALSE, &kf_z, &sigma_kf_z);
ncm_message ("% 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g\n",
z,
kf_z,
kf_z - 1.0 * sigma_kf_z, kf_z + 1.0 * sigma_kf_z,
kf_z - 2.0 * sigma_kf_z, kf_z + 2.0 * sigma_kf_z,
kf_z - 3.0 * sigma_kf_z, kf_z + 3.0 * sigma_kf_z);
}
if (j != 4)
ncm_message ("\n\n");
}
ncm_mset_func_free (dec_param_func);
ncm_mset_func_free (E2_func);
ncm_mset_func_free (Em2_func);
ncm_mset_func_free (dec_func);
ncm_mset_func_free (wec_func);
}
else
{
const guint nz = de_fit.kinematics_n;
NcmMSetFunc *dec_param_func = nc_hicosmo_create_mset_arrayfunc1 (nc_hicosmo_q, nz);
NcmMSetFunc *E2_func = nc_hicosmo_create_mset_arrayfunc1 (nc_hicosmo_E2, nz);
NcmMSetFunc *Em2_func = nc_hicosmo_create_mset_arrayfunc1 (nc_hicosmo_Em2, nz);
NcmMSetFunc *dec_func = nc_hicosmo_create_mset_arrayfunc1 (nc_hicosmo_dec, nz);
NcmMSetFunc *wec_func = nc_hicosmo_create_mset_arrayfunc1 (nc_hicosmo_wec, nz);
NcmVector *z_vec = ncm_vector_new (nz);
GArray *p_val = g_array_new (FALSE, FALSE, sizeof (gdouble));
gint i, j;
struct {
gchar *name;
NcmMSetFunc *func;
} kinematics[5] = {
{"H^2/H_0^2", E2_func},
{"q(z)", dec_param_func},
{"H_0^2/H^2", Em2_func},
{"DEC", dec_func},
{"WEC", wec_func}
};

g_array_set_size (p_val, 3);
g_array_index (p_val, gdouble, 0) = gsl_cdf_chisq_P (1.0, 1.0);
g_array_index (p_val, gdouble, 1) = gsl_cdf_chisq_P (4.0, 1.0);
g_array_index (p_val, gdouble, 2) = gsl_cdf_chisq_P (9.0, 1.0);

ncm_message ("# Kinematics data from catalog:\n");

for (i = 0; i < de_fit.kinematics_n; i++)
{
const gdouble z = de_fit.kinematics_z / (de_fit.kinematics_n - 1.0) * i;
ncm_vector_set (z_vec, i, z);
}

for (j = 0; j < 5; j++)
{
NcmMatrix *res = ncm_mset_catalog_calc_ci (mcat, kinematics[j].func, ncm_vector_ptr (z_vec, 0), p_val);
ncm_message ("# Kinematics data from catalog: %s:\n", kinematics[j].name);

for (i = 0; i < nz; i++)
{

ncm_message ("% 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g\n",
ncm_vector_get (z_vec, i),
ncm_matrix_get (res, i, 0),
ncm_matrix_get (res, i, 1), ncm_matrix_get (res, i, 2),
ncm_matrix_get (res, i, 3), ncm_matrix_get (res, i, 4),
ncm_matrix_get (res, i, 5), ncm_matrix_get (res, i, 6));
}
if (j != 4)
ncm_message ("\n\n");
ncm_matrix_free (res);
}
ncm_vector_free (z_vec);
g_array_unref (p_val);
ncm_mset_func_free (dec_param_func);
ncm_mset_func_free (E2_func);
ncm_mset_func_free (Em2_func);
ncm_mset_func_free (dec_func);
ncm_mset_func_free (wec_func);
}

ncm_message ("# Kinematics data:\n");
for (i = 0; i < de_fit.kinematics_n; i++)
{
gdouble z = de_fit.kinematics_z / (de_fit.kinematics_n - 1.0) * i;
gdouble q_z, E2_z, Em2_z, dec, wec;
gdouble sigma_q_z, sigma_E2_z, sigma_Em2_z, sigma_dec, sigma_wec;
ncm_fit_function_error (fit, dec_param_func, &z, FALSE, &q_z, &sigma_q_z);
ncm_fit_function_error (fit, E2_func, &z, FALSE, &E2_z, &sigma_E2_z);
ncm_fit_function_error (fit, Em2_func, &z, FALSE, &Em2_z, &sigma_Em2_z);
ncm_fit_function_error (fit, dec_func, &z, FALSE, &dec, &sigma_dec);
ncm_fit_function_error (fit, wec_func, &z, FALSE, &wec, &sigma_wec);
ncm_message ("% 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g % 20.15g\n",
z,
q_z, sigma_q_z,
E2_z, sigma_E2_z,
Em2_z, sigma_Em2_z,
dec, sigma_dec,
wec, sigma_wec);
}
ncm_mset_func_free (dec_param_func);
ncm_mset_func_free (E2_func);
ncm_mset_func_free (Em2_func);
}

if (ca_array != NULL)
Expand All @@ -964,7 +1048,8 @@ main (gint argc, gchar *argv[])

if (fiduc != NULL)
ncm_mset_free (fiduc);


ncm_mset_catalog_clear (&mcat);
ncm_serialize_global_reset ();
ncm_model_free (NCM_MODEL (model));
ncm_mset_free (mset);
Expand Down
2 changes: 0 additions & 2 deletions darkenergy/de_options.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,9 +187,7 @@ _nc_de_print_fit_list (const gchar *option_name, const gchar *value, gpointer da

ncm_cfg_enum_print_all (NCM_TYPE_FIT_GSLMM_ALGOS, "Minimization algorithims [gsl-mm]");
ncm_cfg_enum_print_all (NCM_TYPE_FIT_GSLMMS_ALGOS, "Minimization algorithims [gsl-mms]");
#ifdef NUMCOSMO_HAVE_LEVMAR
ncm_cfg_enum_print_all (NCM_TYPE_FIT_LEVMAR_ALGOS, "Minimization algorithims [levmar]");
#endif
#ifdef NUMCOSMO_HAVE_NLOPT
ncm_cfg_enum_print_all (NCM_TYPE_FIT_NLOPT_ALGORITHM, "Minimization algorithims [nlopt]");
#endif
Expand Down

0 comments on commit 0545e3c

Please sign in to comment.