diff --git a/.gitignore b/.gitignore index ba4be7dd5..5f67726c8 100644 --- a/.gitignore +++ b/.gitignore @@ -214,3 +214,15 @@ test/test_SharedPtr test/output_test_serialEnv test/test_Environment/input_test_serialEnv test/test_serialEnv +test/test_jeffreys_pdf +test/test_gaussian_pdf_gradient +test/test_log_normal_pdf_gradient +test/test_beta_pdf_gradient +test/output_test_intercomm0_gravity_1 +test/output_test_intercomm0_gravity_2 +test/test_intercomm0_gravity +test/test_intercomm0/gravity_1proc.txt +test/test_intercomm0/gravity_2proc.txt +test/test_intercomm0/test_intercomm0_gravity_run.sh +test/test_seq_of_vec_hdf5_write +test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh diff --git a/CHANGES b/CHANGES index be1920e46..bac940af3 100644 --- a/CHANGES +++ b/CHANGES @@ -3,7 +3,13 @@ QUESO: Quantification of Uncertainty for Estimation, Simulation, and Optimization. ----------------------------------------------------- -Version 0.54.0 +Version 0.55.0 (in progress) + * Fix bug in JeffreysJointPdf + * Adding gradient support for Beta, Log Normal, and Gaussian PDFs + * Fix MLE/MAP bug when subEnvironments have more than one MPI process + * Add support for HDF5 output + +Version 0.54.0 (Oct 8, 2015) * Fix memory leak in test_GslVector * Make MPI optional * Optimise GSLMatrix::operator() and GslVector::operator[] @@ -13,7 +19,7 @@ Version 0.54.0 * Fix some compiler warnings * Fix HDF5 #include issues -Version 0.53.0 +Version 0.53.0 (Jun 11, 2015) * Add linear interpolation surrogates * Refactor input options processing * Refactor existing queso errors and asserts @@ -21,42 +27,37 @@ Version 0.53.0 * Add basic scoped pointer wrappers * Better error message reporting for bad sample covariance matrices -Version 0.52.0 +Version 0.52.0 (Apr 9, 2015) * Add canned Gaussian likelihoods * Fix bug in logit transform logic -Version 0.51.0 +Version 0.51.0 (Nov 29, 2014) * Add canned likelihood for scalar GPMSA use-case a la Higdon et al * Add a logit-transformed transition kernel for more efficient proposals * Adding jeffreys distribution as an available prior distribution * Adding likelihood value caching to ML sampler * Adding Eric Wright to contributors list -Version 0.50.0 +Version 0.50.0 (Nov 16, 2013) * Implemented MCMC sampling on function spaces * Added function space MCMC unit tests * Added serial and parallel function space MCMC examples * Replace exit() calls with a macro exception handler -Version 0.47.1 +Version 0.47.1 (Sep 23, 2013) * created two simple statistical examples, simpleStatisticalForwardProblem and simpleStatisticalInverseProblem * included description of the new examples on QUESO user's manual, including how to run and to process the generated data * updated QUESO user's manual to include description of validationCycle example * updated links for how to obtain QUESO and contact information (github/git) + * fixed broken Makefile -Version 0.47.0 +Version 0.47.0 (Sep 4, 2013) * removed example statisticalInverseProblem as it has a bug - it fails to recreate the problem described in the Normal test from Laine's MCMC tool: http://helios.fmi.fi/~lainema/mcmc/ex/normalex.html. -Version 0.47.1 (23 Sep 2013) - - * Update manual and fix broken Makefile - Version 0.47.0 (4 Sep 2013) - * Migrated library to github Version 0.46 (13 May 2013) - * changed libGRVY linkage to be optional * added additional regression tests to 'make check' * added rng (random number generator) class to accommodate either GSL or BOOST rng's diff --git a/Makefile.am b/Makefile.am index ec9036893..543fc37d3 100644 --- a/Makefile.am +++ b/Makefile.am @@ -17,6 +17,10 @@ AUX_DIST = build-aux/install-sh \ pkgconfigdir = $(libdir)/pkgconfig pkgconfig_DATA = queso.pc queso-deps.pc +# Support for queso-config in $(exec_prefix)/bin +queso_configdir = $(exec_prefix)/bin +queso_config_SCRIPTS = src/apps/queso-config + if CODE_COVERAGE_ENABLED lcov_dir=$(top_builddir)/docs/html/lcov diff --git a/configure.ac b/configure.ac index 8e2d95948..49a1d70c3 100644 --- a/configure.ac +++ b/configure.ac @@ -2,7 +2,7 @@ # Process this file with autoconf to produce a configure script. AC_PREREQ(2.65) -AC_INIT([queso], [0.54.0], [queso-users@googlegroups.com]) +AC_INIT([queso], [0.55.0], [queso-users@googlegroups.com]) PACKAGE_DESCRIPTION="The parallel C++ statistical library QUESO: Quantification of uncertainty for estimation, simulation and optimization" AC_SUBST([PACKAGE_DESCRIPTION]) PACKAGE_URL="https://github.com/libqueso/queso" @@ -68,8 +68,8 @@ AC_ARG_ENABLE(debug, [ --enable-debug=[no/yes] turn on debugging [default=$debug_default]],, enable_debug=$debug_default) if test "x$enable_debug" = "xyes"; then - CFLAGS="$CFLAGS -DDEBUG" - CXXFLAGS="$CXXFLAGS -DDEBUG" + QUESO_CPPFLAGS="$CPPFLAGS -DDEBUG" + AC_SUBST(QUESO_CPPFLAGS) AC_MSG_RESULT(yes) else AC_MSG_RESULT(no) @@ -241,6 +241,8 @@ AC_CONFIG_FILES([ test/test_InputOptionsParser/test_options_default.txt test/test_optimizer/input_test_optimizer_options test/test_Environment/input_test_serialEnv + test/test_intercomm0/gravity_1proc.txt + test/test_intercomm0/gravity_2proc.txt ]) AC_CONFIG_FILES([test/test_Regression/test_cobra_samples_diff.sh @@ -270,6 +272,22 @@ AC_CONFIG_FILES([ chmod +x test/test_SequenceOfVectors/test_unifiedPositionsOfMaximum.sh ]) +AC_CONFIG_FILES([ + test/test_intercomm0/test_intercomm0_gravity_run.sh +], +[ + chmod +x test/test_intercomm0/test_intercomm0_gravity_run.sh +]) + +AC_CONFIG_FILES([ + test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh +], +[ + chmod +x test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh +]) + +AC_CONFIG_FILES(src/apps/queso-config, [chmod +x src/apps/queso-config]) + dnl ---------------------------------------------- dnl Collect files for licence header stamping here dnl ---------------------------------------------- diff --git a/examples/Makefile.am b/examples/Makefile.am index f742beb3d..0e3adf9fa 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,31 +1,32 @@ ## Process this file with automake to produce Makefile.in -QUESO_CPPFLAGS = -I$(top_builddir)/src/core/inc -QUESO_CPPFLAGS += -I$(top_builddir)/inc -QUESO_CPPFLAGS += $(BOOST_CPPFLAGS) -QUESO_CPPFLAGS += $(GSL_CFLAGS) -QUESO_CPPFLAGS += $(GRVY_CFLAGS) -QUESO_CPPFLAGS += $(ANN_CFLAGS) +AM_CPPFLAGS = $(QUESO_CPPFLAGS) +AM_CPPFLAGS += -I$(top_builddir)/src/core/inc +AM_CPPFLAGS += -I$(top_builddir)/inc +AM_CPPFLAGS += $(BOOST_CPPFLAGS) +AM_CPPFLAGS += $(GSL_CFLAGS) +AM_CPPFLAGS += $(GRVY_CFLAGS) +AM_CPPFLAGS += $(ANN_CFLAGS) LIBS += $(GSL_LIBS) if TRILINOS_ENABLED - QUESO_CPPFLAGS += -I$(TRILINOS_INCLUDE) + AM_CPPFLAGS += -I$(TRILINOS_INCLUDE) LIBS += -lteuchoscore -lteuchoscomm -lteuchosnumerics -lteuchosparameterlist -lteuchosremainder -lepetra endif if GLPK_ENABLED - QUESO_CPPFLAGS += $(GLPK_CFLAGS) + AM_CPPFLAGS += $(GLPK_CFLAGS) LIBS += $(GLPK_LIBS) endif if HDF5_ENABLED - QUESO_CPPFLAGS += $(HDF5_CFLAGS) + AM_CPPFLAGS += $(HDF5_CFLAGS) LIBS += $(HDF5_LIBS) endif if LIBMESH_ENABLED - QUESO_CPPFLAGS += $(LIBMESH_CPPFLAGS) + AM_CPPFLAGS += $(LIBMESH_CPPFLAGS) LIBS += $(LIBMESH_LIBS) endif @@ -38,7 +39,7 @@ operator_DATA = $(top_srcdir)/examples/infinite_dim/operator/mesh.e operator_PROGRAMS = operator operator_SOURCES = infinite_dim/operator/operator.C operator_LDADD = $(top_builddir)/src/libqueso.la -operator_CPPFLAGS = $(QUESO_CPPFLAGS) +operator_CPPFLAGS = $(AM_CPPFLAGS) dist_operator_DATA = $(operator_SOURCES) $(operator_DATA) gaussian_random_fielddir = $(prefix)/examples/infinite_dim/gaussian_fields @@ -46,7 +47,7 @@ gaussian_random_field_DATA = $(top_srcdir)/examples/infinite_dim/gaussian_fields gaussian_random_field_PROGRAMS = gaussian_random_field gaussian_random_field_SOURCES = infinite_dim/gaussian_fields/gaussian_random_field.C gaussian_random_field_LDADD = $(top_builddir)/src/libqueso.la -gaussian_random_field_CPPFLAGS = $(QUESO_CPPFLAGS) +gaussian_random_field_CPPFLAGS = $(AM_CPPFLAGS) dist_gaussian_random_field_DATA = $(gaussian_random_field_SOURCES) $(gaussian_random_field_DATA) gaussian_random_field_genmeshdir = $(prefix)/examples/infinite_dim/gaussian_fields @@ -54,7 +55,7 @@ gaussian_random_field_genmesh_DATA = $(top_srcdir)/examples/infinite_dim/gaussia gaussian_random_field_genmesh_PROGRAMS = gaussian_random_field_genmesh gaussian_random_field_genmesh_SOURCES = infinite_dim/gaussian_fields/gaussian_random_field_genmesh.C gaussian_random_field_genmesh_LDADD = $(top_builddir)/src/libqueso.la -gaussian_random_field_genmesh_CPPFLAGS = $(QUESO_CPPFLAGS) +gaussian_random_field_genmesh_CPPFLAGS = $(AM_CPPFLAGS) dist_gaussian_random_field_genmesh_DATA = $(gaussian_random_field_genmesh_SOURCES) $(gaussian_random_field_genmesh_DATA) invprobdir = $(prefix)/examples/infinite_dim @@ -62,7 +63,7 @@ invprob_DATA = $(top_srcdir)/examples/infinite_dim/inverse_options.in invprob_PROGRAMS = inverse_problem inverse_problem_SOURCES = infinite_dim/inverse_problem.C inverse_problem_LDADD = $(top_builddir)/src/libqueso.la -inverse_problem_CPPFLAGS = $(QUESO_CPPFLAGS) +inverse_problem_CPPFLAGS = $(AM_CPPFLAGS) dist_invprob_DATA = $(inverse_problem_SOURCES) $(invprob_DATA) # Parallel inverse problem example @@ -71,7 +72,7 @@ parinvprob_DATA = $(top_srcdir)/examples/infinite_dim/parallel_inverse_options.i parinvprob_PROGRAMS = parallel_inverse_problem parallel_inverse_problem_SOURCES = infinite_dim/parallel_inverse_problem.C parallel_inverse_problem_LDADD = $(top_builddir)/src/libqueso.la -parallel_inverse_problem_CPPFLAGS = $(QUESO_CPPFLAGS) +parallel_inverse_problem_CPPFLAGS = $(AM_CPPFLAGS) dist_parinvprob_DATA = $(parallel_inverse_problem_SOURCES) $(parinvprob_DATA) endif @@ -95,7 +96,7 @@ exSimpleStatisticalInverseProblem_gsl_HEADERS += simpleStatisticalInverseProblem exSimpleStatisticalInverseProblem_gsl_LDADD = $(top_builddir)/src/libqueso.la -exSimpleStatisticalInverseProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/simpleStatisticalInverseProblem/src $(QUESO_CPPFLAGS) +exSimpleStatisticalInverseProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/simpleStatisticalInverseProblem/src $(AM_CPPFLAGS) # Make sure everything relevant gets added when doing 'make dist' dist_exSimpleStatisticalInverseProblem_gsl_DATA = @@ -122,7 +123,7 @@ exSimpleStatisticalForwardProblem_gsl_HEADERS += simpleStatisticalForwardProblem exSimpleStatisticalForwardProblem_gsl_LDADD = $(top_builddir)/src/libqueso.la -exSimpleStatisticalForwardProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/simpleStatisticalForwardProblem/src $(QUESO_CPPFLAGS) +exSimpleStatisticalForwardProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/simpleStatisticalForwardProblem/src $(AM_CPPFLAGS) dist_exSimpleStatisticalForwardProblem_gsl_DATA = ${exSimpleStatisticalForwardProblem_gsl_DATA} dist_exSimpleStatisticalForwardProblem_gsl_DATA += ${exSimpleStatisticalForwardProblem_gsl_HEADERS} @@ -146,7 +147,7 @@ dist_exSimpleStatisticalForwardProblem_gsl_DATA += ${exSimpleStatisticalForwardP # # exStatisticalForwardProblem_gsl_LDADD = $(top_builddir)/src/libqueso.la # -# exStatisticalForwardProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/statisticalForwardProblem/src $(QUESO_CPPFLAGS) +# exStatisticalForwardProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/statisticalForwardProblem/src $(AM_CPPFLAGS) # #dist_exStatisticalForwardProblem_gsl_DATA = ${exStatisticalForwardProblem_gsl_DATA} # @@ -171,7 +172,7 @@ dist_exSimpleStatisticalForwardProblem_gsl_DATA += ${exSimpleStatisticalForwardP # # exStatisticalInverseProblem_gsl_LDADD = $(top_builddir)/src/libqueso.la # -# exStatisticalInverseProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/statisticalInverseProblem/src $(QUESO_CPPFLAGS) +# exStatisticalInverseProblem_gsl_CPPFLAGS = -I$(top_srcdir)/examples/statisticalInverseProblem/src $(AM_CPPFLAGS) # #dist_exStatisticalInverseProblem_gsl_DATA = ${exStatisticalInverseProblem_gsl_DATA} \ # #${exStatisticalInverseProblem_gsl_HEADERS} \ @@ -200,7 +201,7 @@ exTgaValidationCycle_gsl_HEADERS += validationCycle/src/exTgaValidationCycle_qoi exTgaValidationCycle_gsl_LDADD = $(top_builddir)/src/libqueso.la -exTgaValidationCycle_gsl_CPPFLAGS = -I$(top_srcdir)/examples/validationCycle/src $(QUESO_CPPFLAGS) +exTgaValidationCycle_gsl_CPPFLAGS = -I$(top_srcdir)/examples/validationCycle/src $(AM_CPPFLAGS) dist_exTgaValidationCycle_gsl_DATA = dist_exTgaValidationCycle_gsl_DATA += $(exTgaValidationCycle_gsl_DATA) @@ -236,7 +237,7 @@ tga2_gsl_HEADERS += validationCycle2/src/tga2_func.h tga2_gsl_LDADD = $(top_builddir)/src/libqueso.la -tga2_gsl_CPPFLAGS = -I$(top_srcdir)/examples/validationCycle2/src $(QUESO_CPPFLAGS) +tga2_gsl_CPPFLAGS = -I$(top_srcdir)/examples/validationCycle2/src $(AM_CPPFLAGS) dist_tga2_gsl_DATA = dist_tga2_gsl_DATA += $(tga2_gsl_DATA) @@ -259,7 +260,7 @@ exInfoTheory_gsl_HEADERS = infoTheoryProblem/src/exInfoTheory_appl.h exInfoTheory_gsl_LDADD = $(top_builddir)/src/libqueso.la -exInfoTheory_gsl_CPPFLAGS = -I$(top_srcdir)/examples/infoTheoryProblem/src $(QUESO_CPPFLAGS) +exInfoTheory_gsl_CPPFLAGS = -I$(top_srcdir)/examples/infoTheoryProblem/src $(AM_CPPFLAGS) dist_exInfoTheory_gsl_DATA = dist_exInfoTheory_gsl_DATA += ${exInfoTheory_gsl_DATA} @@ -290,7 +291,7 @@ gravity_gsl_HEADERS += gravity/src/gravity_qoi.h gravity_gsl_LDADD = $(top_builddir)/src/libqueso.la -gravity_gsl_CPPFLAGS = -I$(top_srcdir)/examples/gravity/src $(QUESO_CPPFLAGS) +gravity_gsl_CPPFLAGS = -I$(top_srcdir)/examples/gravity/src $(AM_CPPFLAGS) dist_gravity_gsl_DATA = dist_gravity_gsl_DATA += $(gravity_gsl_DATA) @@ -319,7 +320,7 @@ bimodal_gsl_HEADERS += bimodal/src/bimodal_likelihood.h bimodal_gsl_LDADD = $(top_builddir)/src/libqueso.la -bimodal_gsl_CPPFLAGS = -I$(top_srcdir)/examples/bimodal/src $(QUESO_CPPFLAGS) +bimodal_gsl_CPPFLAGS = -I$(top_srcdir)/examples/bimodal/src $(AM_CPPFLAGS) dist_bimodal_gsl_DATA = dist_bimodal_gsl_DATA += $(bimodal_gsl_DATA) @@ -353,7 +354,7 @@ hysteretic_gsl_HEADERS += hysteretic/src/example_hyst.h hysteretic_gsl_LDADD = $(top_builddir)/src/libqueso.la -hysteretic_gsl_CPPFLAGS = -I$(top_srcdir)/examples/hysteretic/src $(QUESO_CPPFLAGS) +hysteretic_gsl_CPPFLAGS = -I$(top_srcdir)/examples/hysteretic/src $(AM_CPPFLAGS) dist_hysteretic_gsl_DATA = dist_hysteretic_gsl_DATA += $(hysteretic_gsl_DATA) @@ -368,7 +369,7 @@ gpmsa_scalardir = $(prefix)/examples/gp/scalar gpmsa_scalar_PROGRAMS = gpmsa_scalar gpmsa_scalar_SOURCES = gp/scalar/gpmsa_scalar.C gpmsa_scalar_LDADD = $(top_builddir)/src/libqueso.la -gpmsa_scalar_CPPFLAGS = -I$(top_srcdir)/examples/gp/scalar $(QUESO_CPPFLAGS) +gpmsa_scalar_CPPFLAGS = -I$(top_srcdir)/examples/gp/scalar $(AM_CPPFLAGS) dist_gpmsa_scalar_DATA = dist_gpmsa_scalar_DATA += ${gpmsa_scalar_SOURCES} @@ -383,37 +384,37 @@ scalarCovariancedir = $(prefix)/examples/gaussian_likelihoods scalarCovariance_PROGRAMS = scalarCovariance scalarCovariance_SOURCES = gaussian_likelihoods/scalarCovariance.C scalarCovariance_LDADD = $(top_builddir)/src/libqueso.la -scalarCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/scalarCovariance $(QUESO_CPPFLAGS) +scalarCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/scalarCovariance $(AM_CPPFLAGS) diagonalCovariancedir = $(prefix)/examples/gaussian_likelihoods diagonalCovariance_PROGRAMS = diagonalCovariance diagonalCovariance_SOURCES = gaussian_likelihoods/diagonalCovariance.C diagonalCovariance_LDADD = $(top_builddir)/src/libqueso.la -diagonalCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/diagonalCovariance $(QUESO_CPPFLAGS) +diagonalCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/diagonalCovariance $(AM_CPPFLAGS) fullCovariancedir = $(prefix)/examples/gaussian_likelihoods fullCovariance_PROGRAMS = fullCovariance fullCovariance_SOURCES = gaussian_likelihoods/fullCovariance.C fullCovariance_LDADD = $(top_builddir)/src/libqueso.la -fullCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/fullCovariance $(QUESO_CPPFLAGS) +fullCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/fullCovariance $(AM_CPPFLAGS) blockDiagonalCovariancedir = $(prefix)/examples/gaussian_likelihoods blockDiagonalCovariance_PROGRAMS = blockDiagonalCovariance blockDiagonalCovariance_SOURCES = gaussian_likelihoods/blockDiagonalCovariance.C blockDiagonalCovariance_LDADD = $(top_builddir)/src/libqueso.la -blockDiagonalCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/blockDiagonalCovariance $(QUESO_CPPFLAGS) +blockDiagonalCovariance_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/blockDiagonalCovariance $(AM_CPPFLAGS) fullCovarianceRandomCoefficientdir = $(prefix)/examples/gaussian_likelihoods fullCovarianceRandomCoefficient_PROGRAMS = fullCovarianceRandomCoefficient fullCovarianceRandomCoefficient_SOURCES = gaussian_likelihoods/fullCovarianceRandomCoefficient.C fullCovarianceRandomCoefficient_LDADD = $(top_builddir)/src/libqueso.la -fullCovarianceRandomCoefficient_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/fullCovarianceRandomCoefficient $(QUESO_CPPFLAGS) +fullCovarianceRandomCoefficient_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/fullCovarianceRandomCoefficient $(AM_CPPFLAGS) blockDiagonalCovarianceRandomCoefficientsdir = $(prefix)/examples/gaussian_likelihoods blockDiagonalCovarianceRandomCoefficients_PROGRAMS = blockDiagonalCovarianceRandomCoefficients blockDiagonalCovarianceRandomCoefficients_SOURCES = gaussian_likelihoods/blockDiagonalCovarianceRandomCoefficients.C blockDiagonalCovarianceRandomCoefficients_LDADD = $(top_builddir)/src/libqueso.la -blockDiagonalCovarianceRandomCoefficients_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/blockDiagonalCovarianceRandomCoefficients $(QUESO_CPPFLAGS) +blockDiagonalCovarianceRandomCoefficients_CPPFLAGS = -I$(top_srcdir)/examples/gaussian_likelihoods/blockDiagonalCovarianceRandomCoefficients $(AM_CPPFLAGS) #################################### # Interpolation Surrogate examples # @@ -424,11 +425,11 @@ interpsurrg_PROGRAMS = 4d_interp 4d_interp_read 4d_interp_SOURCES = interpolation_surrogate/4d_interp.C 4d_interp_LDADD = $(top_builddir)/src/libqueso.la -4d_interp_CPPFLAGS = $(QUESO_CPPFLAGS) +4d_interp_CPPFLAGS = $(AM_CPPFLAGS) 4d_interp_read_SOURCES = interpolation_surrogate/4d_interp_read.C 4d_interp_read_LDADD = $(top_builddir)/src/libqueso.la -4d_interp_read_CPPFLAGS = $(QUESO_CPPFLAGS) +4d_interp_read_CPPFLAGS = $(AM_CPPFLAGS) dist_interpsurrg_DATA = dist_interpsurrg_DATA += ${4d_interp_SOURCES} @@ -445,7 +446,7 @@ pointers_PROGRAMS = example_SharedPtr example_SharedPtr_SOURCES = pointers/example_SharedPtr.C example_SharedPtr_LDADD = $(top_builddir)/src/libqueso.la -example_SharedPtr_CPPFLAGS = $(QUESO_CPPFLAGS) +example_SharedPtr_CPPFLAGS = $(AM_CPPFLAGS) dist_pointers_DATA = dist_pointers_DATA += ${example_SharedPtr_SOURCES} diff --git a/m4/common/boost.m4 b/m4/common/boost.m4 index 3d54bfd2d..7a9298151 100644 --- a/m4/common/boost.m4 +++ b/m4/common/boost.m4 @@ -1,5 +1,5 @@ # boost.m4: Locate Boost headers and libraries for autoconf-based projects. -# Copyright (C) 2007, 2008, 2009, 2010, 2011 Benoit Sigoure +# Copyright (C) 2007-2011, 2014 Benoit Sigoure # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -22,7 +22,7 @@ # along with this program. If not, see . m4_define([_BOOST_SERIAL], [m4_translit([ -# serial 18 +# serial 25 ], [# ], [])]) @@ -59,7 +59,8 @@ m4_pattern_forbid([^_?(BOOST|Boost)_]) # It could be useful to turn this into a macro which extracts the # value of any macro. m4_define([_BOOST_SED_CPP], -[AC_LANG_PREPROC_REQUIRE()dnl +[AC_LANG_PUSH([C++])dnl +AC_LANG_PREPROC_REQUIRE()dnl AC_REQUIRE([AC_PROG_SED])dnl AC_LANG_CONFTEST([AC_LANG_SOURCE([[$2]])]) AS_IF([dnl eval is necessary to expand ac_cpp. @@ -71,13 +72,31 @@ dnl strip `\n' with backquotes, not the `\r'. This results in dnl boost_cv_lib_version='1_37\r' for instance, which breaks dnl everything else. dnl Cannot use 'dnl' after [$4] because a trailing dnl may break AC_CACHE_CHECK +dnl +dnl Beware that GCC 5, when expanding macros, may embed # line directives +dnl a within single line: +dnl +dnl # 1 "conftest.cc" +dnl # 1 "" +dnl # 1 "" +dnl # 1 "conftest.cc" +dnl # 1 "/opt/local/include/boost/version.hpp" 1 3 +dnl # 2 "conftest.cc" 2 +dnl boost-lib-version = +dnl # 2 "conftest.cc" 3 +dnl "1_56" +dnl +dnl So get rid of the # lines, and glue the remaining ones together. (eval "$ac_cpp conftest.$ac_ext") 2>&AS_MESSAGE_LOG_FD | + grep -v '#' | tr -d '\r' | + tr -s '\n' ' ' | $SED -n -e "$1" >conftest.i 2>&1], [$3], [$4]) rm -rf conftest* -])# AC_EGREP_CPP +AC_LANG_POP([C++])dnl +])# _BOOST_SED_CPP @@ -91,7 +110,7 @@ rm -rf conftest* # On # success, defines HAVE_BOOST. On failure, calls the optional # ACTION-IF-NOT-FOUND action if one was supplied. # Otherwise aborts with an error message. -AC_DEFUN([BOOST_REQUIRE], +AC_DEFUN_ONCE([BOOST_REQUIRE], [AC_REQUIRE([AC_PROG_CXX])dnl AC_REQUIRE([AC_PROG_GREP])dnl echo "$as_me: this is boost.m4[]_BOOST_SERIAL" >&AS_MESSAGE_LOG_FD @@ -107,13 +126,6 @@ AC_ARG_WITH([boost], [AS_HELP_STRING([--with-boost=DIR], [prefix of Boost $1 @<:@guess@:>@])])dnl AC_ARG_VAR([BOOST_ROOT],[Location of Boost installation])dnl -# If BOOST_DIR is set and BOOST_ROOT is not, use BOOST_DIR. -m4_pattern_allow([^BOOST_DIR$])dnl -if test x"$BOOST_DIR" != x; then - if test x"$BOOST_ROOT" = x; then - BOOST_ROOT="$BOOST_DIR" - fi -fi # If BOOST_ROOT is set and the user has not provided a value to # --with-boost, then treat BOOST_ROOT as if it the user supplied it. if test x"$BOOST_ROOT" != x; then @@ -213,7 +225,7 @@ AC_LANG_POP([C++])dnl AC_CACHE_CHECK([for Boost's header version], [boost_cv_lib_version], [m4_pattern_allow([^BOOST_LIB_VERSION$])dnl - _BOOST_SED_CPP([/^boost-lib-version = /{s///;s/\"//g;p;q;}], + _BOOST_SED_CPP([[/^boost-lib-version = /{s///;s/[\" ]//g;p;q;}]], [#include boost-lib-version = BOOST_LIB_VERSION], [boost_cv_lib_version=`cat conftest.i`])]) @@ -221,13 +233,14 @@ boost-lib-version = BOOST_LIB_VERSION], boost_major_version=`echo "$boost_cv_lib_version" | sed 's/_//;s/_.*//'` case $boost_major_version in #( '' | *[[!0-9]]*) - AC_MSG_ERROR([invalid value: boost_major_version=$boost_major_version]) + AC_MSG_ERROR([invalid value: boost_major_version='$boost_major_version']) ;; esac fi CPPFLAGS=$boost_save_CPPFLAGS ])# BOOST_REQUIRE + # BOOST_STATIC() # -------------- # Add the "--enable-static-boost" configure argument. If this argument is given @@ -239,6 +252,7 @@ AC_DEFUN([BOOST_STATIC], [enable_static_boost=yes], [enable_static_boost=no])])# BOOST_STATIC + # BOOST_FIND_HEADER([HEADER-NAME], [ACTION-IF-NOT-FOUND], [ACTION-IF-FOUND]) # -------------------------------------------------------------------------- # Wrapper around AC_CHECK_HEADER for Boost headers. Useful to check for @@ -271,14 +285,16 @@ fi ])# BOOST_FIND_HEADER -# BOOST_FIND_LIB([LIB-NAME], [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST], -# [CXX-PROLOGUE]) -# ------------------------------------------------------------------------- -# Look for the Boost library LIB-NAME (e.g., LIB-NAME = `thread', for -# libboost_thread). Check that HEADER-NAME works and check that -# libboost_LIB-NAME can link with the code CXX-TEST. The optional argument -# CXX-PROLOGUE can be used to include some C++ code before the `main' -# function. +# BOOST_FIND_LIBS([COMPONENT-NAME], [CANDIDATE-LIB-NAMES], +# [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST], +# [CXX-PROLOGUE]) +# -------------------------------------------------------------- +# Look for the Boost library COMPONENT-NAME (e.g., `thread', for +# libboost_thread) under the possible CANDIDATE-LIB-NAMES (e.g., +# "thread_win32 thread"). Check that HEADER-NAME works and check that +# libboost_LIB-NAME can link with the code CXX-TEST. The optional +# argument CXX-PROLOGUE can be used to include some C++ code before +# the `main' function. # # Invokes BOOST_FIND_HEADER([HEADER-NAME]) (see above). # @@ -292,7 +308,7 @@ fi # builds. Some sample values for PREFERRED-RT-OPT: (nothing), mt, d, mt-d, gdp # ... If you want to make sure you have a specific version of Boost # (eg, >= 1.33) you *must* invoke BOOST_REQUIRE before this macro. -AC_DEFUN([BOOST_FIND_LIB], +AC_DEFUN([BOOST_FIND_LIBS], [AC_REQUIRE([BOOST_REQUIRE])dnl AC_REQUIRE([_BOOST_FIND_COMPILER_TAG])dnl AC_REQUIRE([BOOST_STATIC])dnl @@ -306,32 +322,69 @@ AS_VAR_PUSHDEF([Boost_lib], [boost_cv_lib_$1])dnl AS_VAR_PUSHDEF([Boost_lib_LDFLAGS], [boost_cv_lib_$1_LDFLAGS])dnl AS_VAR_PUSHDEF([Boost_lib_LDPATH], [boost_cv_lib_$1_LDPATH])dnl AS_VAR_PUSHDEF([Boost_lib_LIBS], [boost_cv_lib_$1_LIBS])dnl -BOOST_FIND_HEADER([$3]) +BOOST_FIND_HEADER([$4]) boost_save_CPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS $BOOST_CPPFLAGS" -# Now let's try to find the library. The algorithm is as follows: first look -# for a given library name according to the user's PREFERRED-RT-OPT. For each -# library name, we prefer to use the ones that carry the tag (toolset name). -# Each library is searched through the various standard paths were Boost is -# usually installed. If we can't find the standard variants, we try to -# enforce -mt (for instance on MacOSX, libboost_threads.dylib doesn't exist -# but there's -obviously- libboost_threads-mt.dylib). AC_CACHE_CHECK([for the Boost $1 library], [Boost_lib], - [Boost_lib=no - case "$2" in #( - mt | mt-) boost_mt=-mt; boost_rtopt=;; #( - mt* | mt-*) boost_mt=-mt; boost_rtopt=`expr "X$2" : 'Xmt-*\(.*\)'`;; #( - *) boost_mt=; boost_rtopt=$2;; + [_BOOST_FIND_LIBS($@)]) +case $Boost_lib in #( + (no) _AC_MSG_LOG_CONFTEST + AC_MSG_ERROR([cannot find the flags to link with Boost $1]) + ;; +esac +AC_SUBST(AS_TR_CPP([BOOST_$1_LDFLAGS]), [$Boost_lib_LDFLAGS])dnl +AC_SUBST(AS_TR_CPP([BOOST_$1_LDPATH]), [$Boost_lib_LDPATH])dnl +AC_SUBST([BOOST_LDPATH], [$Boost_lib_LDPATH])dnl +AC_SUBST(AS_TR_CPP([BOOST_$1_LIBS]), [$Boost_lib_LIBS])dnl +CPPFLAGS=$boost_save_CPPFLAGS +AS_VAR_POPDEF([Boost_lib])dnl +AS_VAR_POPDEF([Boost_lib_LDFLAGS])dnl +AS_VAR_POPDEF([Boost_lib_LDPATH])dnl +AS_VAR_POPDEF([Boost_lib_LIBS])dnl +AC_LANG_POP([C++])dnl +fi +]) + + +# BOOST_FIND_LIB([LIB-NAME], +# [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST], +# [CXX-PROLOGUE]) +# -------------------------------------------------------------- +# Backward compatibility wrapper for BOOST_FIND_LIBS. +AC_DEFUN([BOOST_FIND_LIB], +[BOOST_FIND_LIBS([$1], $@)]) + + +# _BOOST_FIND_LIBS([LIB-NAME], [CANDIDATE-LIB-NAMES], +# [PREFERRED-RT-OPT], [HEADER-NAME], [CXX-TEST], +# [CXX-PROLOGUE]) +# -------------------------------------------------------------- +# Real implementation of BOOST_FIND_LIBS: rely on these local macros: +# Boost_lib, Boost_lib_LDFLAGS, Boost_lib_LDPATH, Boost_lib_LIBS +# +# The algorithm is as follows: first look for a given library name +# according to the user's PREFERRED-RT-OPT. For each library name, we +# prefer to use the ones that carry the tag (toolset name). Each +# library is searched through the various standard paths were Boost is +# usually installed. If we can't find the standard variants, we try +# to enforce -mt (for instance on MacOSX, libboost_thread.dylib +# doesn't exist but there's -obviously- libboost_thread-mt.dylib). +AC_DEFUN([_BOOST_FIND_LIBS], +[Boost_lib=no + case "$3" in #( + (mt | mt-) boost_mt=-mt; boost_rtopt=;; #( + (mt* | mt-*) boost_mt=-mt; boost_rtopt=`expr "X$3" : 'Xmt-*\(.*\)'`;; #( + (*) boost_mt=; boost_rtopt=$3;; esac if test $enable_static_boost = yes; then boost_rtopt="s$boost_rtopt" fi # Find the proper debug variant depending on what we've been asked to find. case $boost_rtopt in #( - *d*) boost_rt_d=$boost_rtopt;; #( - *[[sgpn]]*) # Insert the `d' at the right place (in between `sg' and `pn') + (*d*) boost_rt_d=$boost_rtopt;; #( + (*[[sgpn]]*) # Insert the `d' at the right place (in between `sg' and `pn') boost_rt_d=`echo "$boost_rtopt" | sed 's/\(s*g*\)\(p*n*\)/\1\2/'`;; #( - *) boost_rt_d='-d';; + (*) boost_rt_d='-d';; esac # If the PREFERRED-RT-OPT are not empty, prepend a `-'. test -n "$boost_rtopt" && boost_rtopt="-$boost_rtopt" @@ -342,8 +395,8 @@ AC_CACHE_CHECK([for the Boost $1 library], [Boost_lib], AC_MSG_ERROR([the libext variable is empty, did you invoke Libtool?]) boost_save_ac_objext=$ac_objext # Generate the test file. - AC_LANG_CONFTEST([AC_LANG_PROGRAM([#include <$3> -$5], [$4])]) + AC_LANG_CONFTEST([AC_LANG_PROGRAM([#include <$4> +$6], [$5])]) dnl Optimization hacks: compiling C++ is slow, especially with Boost. What dnl we're trying to do here is guess the right combination of link flags dnl (LIBS / LDFLAGS) to use a given library. This can take several @@ -365,21 +418,22 @@ dnl start the for loops). [AC_MSG_ERROR([cannot compile a test that uses Boost $1])]) ac_objext=$boost_save_ac_objext boost_failed_libs= -# Don't bother to ident the 6 nested for loops, only the 2 innermost ones -# matter. +# Don't bother to ident the following nested for loops, only the 2 +# innermost ones matter. +for boost_lib_ in $2; do for boost_tag_ in -$boost_cv_lib_tag ''; do for boost_ver_ in -$boost_cv_lib_version ''; do for boost_mt_ in $boost_mt -mt ''; do for boost_rtopt_ in $boost_rtopt '' -d; do for boost_lib in \ - boost_$1$boost_tag_$boost_mt_$boost_rtopt_$boost_ver_ \ - boost_$1$boost_tag_$boost_rtopt_$boost_ver_ \ - boost_$1$boost_tag_$boost_mt_$boost_ver_ \ - boost_$1$boost_tag_$boost_ver_ + boost_$boost_lib_$boost_tag_$boost_mt_$boost_rtopt_$boost_ver_ \ + boost_$boost_lib_$boost_tag_$boost_rtopt_$boost_ver_ \ + boost_$boost_lib_$boost_tag_$boost_mt_$boost_ver_ \ + boost_$boost_lib_$boost_tag_$boost_ver_ do # Avoid testing twice the same lib case $boost_failed_libs in #( - *@$boost_lib@*) continue;; + (*@$boost_lib@*) continue;; esac # If with_boost is empty, we'll search in /lib first, which is not quite # right so instead we'll try to a location based on where the headers are. @@ -389,14 +443,17 @@ for boost_rtopt_ in $boost_rtopt '' -d; do /opt/local/lib* /usr/local/lib* /opt/lib* /usr/lib* \ "$with_boost" C:/Boost/lib /lib* do - test -e "$boost_ldpath" || continue + # Don't waste time with directories that don't exist. + if test x"$boost_ldpath" != x && test ! -e "$boost_ldpath"; then + continue + fi boost_save_LDFLAGS=$LDFLAGS # Are we looking for a static library? case $boost_ldpath:$boost_rtopt_ in #( - *?*:*s*) # Yes (Non empty boost_ldpath + s in rt opt) + (*?*:*s*) # Yes (Non empty boost_ldpath + s in rt opt) Boost_lib_LIBS="$boost_ldpath/lib$boost_lib.$libext" test -e "$Boost_lib_LIBS" || continue;; #( - *) # No: use -lboost_foo to find the shared library. + (*) # No: use -lboost_foo to find the shared library. Boost_lib_LIBS="-l$boost_lib";; esac boost_save_LIBS=$LIBS @@ -410,27 +467,35 @@ dnl generated only once above (before we start the for loops). LDFLAGS=$boost_save_LDFLAGS LIBS=$boost_save_LIBS if test x"$Boost_lib" = xyes; then - # Check or used cached result of whether or not using -R or -rpath makes sense. - # Some implementations of ld, such as for Mac OSX, require -rpath but - # -R is the flag known to work on other systems. - # https://github.com/tsuna/boost.m4/issues/19 + # Check or used cached result of whether or not using -R or + # -rpath makes sense. Some implementations of ld, such as for + # Mac OSX, require -rpath but -R is the flag known to work on + # other systems. https://github.com/tsuna/boost.m4/issues/19 AC_CACHE_VAL([boost_cv_rpath_link_ldflag], - [for boost_cv_rpath_link_ldflag in -Wl,-R, -Wl,-rpath,; do - LDFLAGS="$boost_save_LDFLAGS -L$boost_ldpath $boost_cv_rpath_link_ldflag$boost_ldpath" - LIBS="$boost_save_LIBS $Boost_lib_LIBS" - _BOOST_AC_LINK_IFELSE([], - [boost_rpath_link_ldflag_found=yes - break], - [boost_rpath_link_ldflag_found=no]) - done + [case $boost_ldpath in + '') # Nothing to do. + boost_cv_rpath_link_ldflag= + boost_rpath_link_ldflag_found=yes;; + *) + for boost_cv_rpath_link_ldflag in -Wl,-R, -Wl,-rpath,; do + LDFLAGS="$boost_save_LDFLAGS -L$boost_ldpath $boost_cv_rpath_link_ldflag$boost_ldpath" + LIBS="$boost_save_LIBS $Boost_lib_LIBS" + _BOOST_AC_LINK_IFELSE([], + [boost_rpath_link_ldflag_found=yes + break], + [boost_rpath_link_ldflag_found=no]) + done + ;; + esac AS_IF([test "x$boost_rpath_link_ldflag_found" != "xyes"], [AC_MSG_ERROR([Unable to determine whether to use -R or -rpath])]) LDFLAGS=$boost_save_LDFLAGS LIBS=$boost_save_LIBS ]) - Boost_lib_LDFLAGS="-L$boost_ldpath $boost_cv_rpath_link_ldflag$boost_ldpath" + test x"$boost_ldpath" != x && + Boost_lib_LDFLAGS="-L$boost_ldpath $boost_cv_rpath_link_ldflag$boost_ldpath" Boost_lib_LDPATH="$boost_ldpath" - break 6 + break 7 else boost_failed_libs="$boost_failed_libs@$boost_lib@" fi @@ -440,25 +505,10 @@ done done done done +done # boost_lib_ rm -f conftest.$ac_objext ]) -case $Boost_lib in #( - no) _AC_MSG_LOG_CONFTEST - AC_MSG_ERROR([cannot find the flags to link with Boost $1]) - ;; -esac -AC_SUBST(AS_TR_CPP([BOOST_$1_LDFLAGS]), [$Boost_lib_LDFLAGS])dnl -AC_SUBST(AS_TR_CPP([BOOST_$1_LDPATH]), [$Boost_lib_LDPATH])dnl -AC_SUBST([BOOST_LDPATH], [$Boost_lib_LDPATH])dnl -AC_SUBST(AS_TR_CPP([BOOST_$1_LIBS]), [$Boost_lib_LIBS])dnl -CPPFLAGS=$boost_save_CPPFLAGS -AS_VAR_POPDEF([Boost_lib])dnl -AS_VAR_POPDEF([Boost_lib_LDFLAGS])dnl -AS_VAR_POPDEF([Boost_lib_LDPATH])dnl -AS_VAR_POPDEF([Boost_lib_LIBS])dnl -AC_LANG_POP([C++])dnl -fi -])# BOOST_FIND_LIB + # --------------------------------------- # @@ -500,20 +550,20 @@ BOOST_FIND_HEADER([boost/asio.hpp])]) # BOOST_BIND() # ------------ -# Look for Boost.Bind +# Look for Boost.Bind. BOOST_DEFUN([Bind], [BOOST_FIND_HEADER([boost/bind.hpp])]) # BOOST_CHRONO() -# ------------------ -# Look for Boost.Chrono +# -------------- +# Look for Boost.Chrono. BOOST_DEFUN([Chrono], [# Do we have to check for Boost.System? This link-time dependency was # added as of 1.35.0. If we have a version <1.35, we must not attempt to # find Boost.System as it didn't exist by then. if test $boost_major_version -ge 135; then -BOOST_SYSTEM([$1]) + BOOST_SYSTEM([$1]) fi # end of the Boost.System check. boost_filesystem_save_LIBS=$LIBS boost_filesystem_save_LDFLAGS=$LDFLAGS @@ -522,16 +572,42 @@ LIBS="$LIBS $BOOST_SYSTEM_LIBS" LDFLAGS="$LDFLAGS $BOOST_SYSTEM_LDFLAGS" BOOST_FIND_LIB([chrono], [$1], [boost/chrono.hpp], - [boost::chrono::high_resolution_clock d;]) + [boost::chrono::thread_clock d;]) if test $enable_static_boost = yes && test $boost_major_version -ge 135; then - AC_SUBST([BOOST_FILESYSTEM_LIBS], ["$BOOST_FILESYSTEM_LIBS $BOOST_SYSTEM_LIBS"]) + BOOST_CHRONO_LIBS="$BOOST_CHRONO_LIBS $BOOST_SYSTEM_LIBS" fi LIBS=$boost_filesystem_save_LIBS LDFLAGS=$boost_filesystem_save_LDFLAGS - ])# BOOST_CHRONO +# BOOST_CONTEXT([PREFERRED-RT-OPT]) +# ----------------------------------- +# Look for Boost.Context. For the documentation of PREFERRED-RT-OPT, see the +# documentation of BOOST_FIND_LIB above. This library was introduced in Boost +# 1.51.0 +BOOST_DEFUN([Context], +[BOOST_FIND_LIB([context], [$1], + [boost/context/all.hpp],[[ +// creates a stack +void * stack_pointer = new void*[4096]; +std::size_t const size = sizeof(void*[4096]); + +// context fc uses f() as context function +// fcontext_t is placed on top of context stack +// a pointer to fcontext_t is returned +fc = ctx::make_fcontext(stack_pointer, size, f); +return ctx::jump_fcontext(&fcm, fc, 3) == 6;]],[dnl +namespace ctx = boost::context; +// context +static ctx::fcontext_t fcm, *fc; +// context-function +static void f(intptr_t i) { + ctx::jump_fcontext(fc, &fcm, i * 2); +}]) +])# BOOST_CONTEXT + + # BOOST_CONVERSION() # ------------------ # Look for Boost.Conversion (cast / lexical_cast) @@ -541,6 +617,44 @@ BOOST_FIND_HEADER([boost/lexical_cast.hpp]) ])# BOOST_CONVERSION +# BOOST_COROUTINE([PREFERRED-RT-OPT]) +# ----------------------------------- +# Look for Boost.Coroutine. For the documentation of PREFERRED-RT-OPT, see the +# documentation of BOOST_FIND_LIB above. This library was introduced in Boost +# 1.53.0 +BOOST_DEFUN([Coroutine], +[ +boost_coroutine_save_LIBS=$LIBS +boost_coroutine_save_LDFLAGS=$LDFLAGS +# Link-time dependency from coroutine to context +BOOST_CONTEXT([$1]) +# Starting from Boost 1.55 a dependency on Boost.System is added +if test $boost_major_version -ge 155; then + BOOST_SYSTEM([$1]) +fi +m4_pattern_allow([^BOOST_(CONTEXT|SYSTEM)_(LIBS|LDFLAGS)]) +LIBS="$LIBS $BOOST_CONTEXT_LIBS $BOOST_SYSTEM_LIBS" +LDFLAGS="$LDFLAGS $BOOST_CONTEXT_LDFLAGS" + +BOOST_FIND_LIB([coroutine], [$1], + [boost/coroutine/coroutine.hpp], + [boost::coroutines::coroutine< int(int) > coro; coro.empty();]) + +# Link-time dependency from coroutine to context, existed only in 1.53, in 1.54 +# coroutine doesn't use context from its headers but from its library. +if test $boost_major_version -eq 153 || test $enable_static_boost = yes && test $boost_major_version -ge 154; then + BOOST_COROUTINE_LIBS="$BOOST_COROUTINE_LIBS $BOOST_CONTEXT_LIBS" + BOOST_COROUTINE_LDFLAGS="$BOOST_COROUTINE_LDFLAGS $BOOST_CONTEXT_LDFLAGS" +fi +if test $enable_static_boost = yes && test $boost_major_version -ge 155; then + BOOST_COROUTINE_LIBS="$BOOST_COROUTINE_LIBS $BOOST_SYSTEM_LIBS" + BOOST_COROUTINE_LDFLAGS="$BOOST_COROUTINE_LDFLAGS $BOOST_SYSTEM_LDFLAGS" +fi +LIBS=$boost_coroutine_save_LIBS +LDFLAGS=$boost_coroutine_save_LDFLAGS +])# BOOST_COROUTINE + + # BOOST_CRC() # ----------- # Look for Boost.CRC @@ -571,7 +685,7 @@ BOOST_DEFUN([Filesystem], # added as of 1.35.0. If we have a version <1.35, we must not attempt to # find Boost.System as it didn't exist by then. if test $boost_major_version -ge 135; then -BOOST_SYSTEM([$1]) + BOOST_SYSTEM([$1]) fi # end of the Boost.System check. boost_filesystem_save_LIBS=$LIBS boost_filesystem_save_LDFLAGS=$LDFLAGS @@ -581,23 +695,34 @@ LDFLAGS="$LDFLAGS $BOOST_SYSTEM_LDFLAGS" BOOST_FIND_LIB([filesystem], [$1], [boost/filesystem/path.hpp], [boost::filesystem::path p;]) if test $enable_static_boost = yes && test $boost_major_version -ge 135; then - AC_SUBST([BOOST_FILESYSTEM_LIBS], ["$BOOST_FILESYSTEM_LIBS $BOOST_SYSTEM_LIBS"]) + BOOST_FILESYSTEM_LIBS="$BOOST_FILESYSTEM_LIBS $BOOST_SYSTEM_LIBS" fi LIBS=$boost_filesystem_save_LIBS LDFLAGS=$boost_filesystem_save_LDFLAGS ])# BOOST_FILESYSTEM +# BOOST_FLYWEIGHT() +# ----------------- +# Look for Boost.Flyweight. +BOOST_DEFUN([Flyweight], +[dnl There's a hidden dependency on pthreads. +AC_REQUIRE([_BOOST_PTHREAD_FLAG])dnl +BOOST_FIND_HEADER([boost/flyweight.hpp]) +AC_SUBST([BOOST_FLYWEIGHT_LIBS], [$boost_cv_pthread_flag]) +]) + + # BOOST_FOREACH() # --------------- -# Look for Boost.Foreach +# Look for Boost.Foreach. BOOST_DEFUN([Foreach], [BOOST_FIND_HEADER([boost/foreach.hpp])]) # BOOST_FORMAT() # -------------- -# Look for Boost.Format +# Look for Boost.Format. # Note: we can't check for boost/format/format_fwd.hpp because the header isn't # standalone. It can't be compiled because it triggers the following error: # boost/format/detail/config_macros.hpp:88: error: 'locale' in namespace 'std' @@ -656,9 +781,18 @@ BOOST_DEFUN([Lambda], [BOOST_FIND_HEADER([boost/lambda/lambda.hpp])]) +# BOOST_LOCALE() +# -------------- +# Look for Boost.Locale +BOOST_DEFUN([Locale], +[BOOST_FIND_LIB([locale], [$1], + [boost/locale.hpp], + [[boost::locale::generator gen; std::locale::global(gen(""));]]) +])# BOOST_LOCALE + # BOOST_LOG([PREFERRED-RT-OPT]) # ----------------------------- -# Look for Boost.Log For the documentation of PREFERRED-RT-OPT, see the +# Look for Boost.Log. For the documentation of PREFERRED-RT-OPT, see the # documentation of BOOST_FIND_LIB above. BOOST_DEFUN([Log], [BOOST_FIND_LIB([log], [$1], @@ -669,12 +803,12 @@ BOOST_DEFUN([Log], # BOOST_LOG_SETUP([PREFERRED-RT-OPT]) # ----------------------------------- -# Look for Boost.Log For the documentation of PREFERRED-RT-OPT, see the +# Look for Boost.Log. For the documentation of PREFERRED-RT-OPT, see the # documentation of BOOST_FIND_LIB above. BOOST_DEFUN([Log_Setup], [AC_REQUIRE([BOOST_LOG])dnl BOOST_FIND_LIB([log_setup], [$1], - [boost/log/utility/init/from_settings.hpp], + [boost/log/utility/setup/from_settings.hpp], [boost::log::basic_settings bs; bs.empty();]) ])# BOOST_LOG_SETUP @@ -691,6 +825,29 @@ BOOST_DEFUN([Math], [BOOST_FIND_HEADER([boost/math/special_functions.hpp])]) +# BOOST_MPI([PREFERRED-RT-OPT]) +# ------------------------------- +# Look for Boost MPI. For the documentation of PREFERRED-RT-OPT, see the +# documentation of BOOST_FIND_LIB above. Uses MPICXX variable if it is +# set, otherwise tries CXX +# +BOOST_DEFUN([MPI], +[boost_save_CXX=${CXX} +boost_save_CXXCPP=${CXXCPP} +if test x"${MPICXX}" != x; then + CXX=${MPICXX} + CXXCPP="${MPICXX} -E" +fi +BOOST_FIND_LIB([mpi], [$1], + [boost/mpi.hpp], + [int argc = 0; + char **argv = 0; + boost::mpi::environment env(argc,argv);]) +CXX=${boost_save_CXX} +CXXCPP=${boost_save_CXXCPP} +])# BOOST_MPI + + # BOOST_MULTIARRAY() # ------------------ # Look for Boost.MultiArray @@ -728,6 +885,12 @@ BOOST_DEFUN([Preprocessor], [BOOST_FIND_HEADER([boost/preprocessor/repeat.hpp])]) +# BOOST_RANGE() +# -------------------- +# Look for Boost.Range +BOOST_DEFUN([Range], +[BOOST_FIND_HEADER([boost/range/adaptors.hpp])]) + # BOOST_UNORDERED() # ----------------- # Look for Boost.Unordered @@ -774,9 +937,9 @@ BOOST_DEFUN([Python], _BOOST_PYTHON_CONFIG([LDFLAGS], [ldflags]) _BOOST_PYTHON_CONFIG([LIBS], [libs]) m4_pattern_allow([^BOOST_PYTHON_MODULE$])dnl -BOOST_FIND_LIB([python], [$1], - [boost/python.hpp], - [], [BOOST_PYTHON_MODULE(empty) {}]) +BOOST_FIND_LIBS([python], [python python3], [$1], + [boost/python.hpp], + [], [BOOST_PYTHON_MODULE(empty) {}]) CPPFLAGS=$boost_python_save_CPPFLAGS LDFLAGS=$boost_python_save_LDFLAGS LIBS=$boost_python_save_LIBS @@ -882,19 +1045,18 @@ BOOST_FIND_LIB([unit_test_framework], [$1], ])# BOOST_TEST -# BOOST_THREADS([PREFERRED-RT-OPT]) +# BOOST_THREAD([PREFERRED-RT-OPT]) # --------------------------------- # Look for Boost.Thread. For the documentation of PREFERRED-RT-OPT, see the # documentation of BOOST_FIND_LIB above. -# FIXME: Provide an alias "BOOST_THREAD". -BOOST_DEFUN([Threads], +BOOST_DEFUN([Thread], [dnl Having the pthread flag is required at least on GCC3 where dnl boost/thread.hpp would complain if we try to compile without dnl -pthread on GNU/Linux. AC_REQUIRE([_BOOST_PTHREAD_FLAG])dnl -boost_threads_save_LIBS=$LIBS -boost_threads_save_LDFLAGS=$LDFLAGS -boost_threads_save_CPPFLAGS=$CPPFLAGS +boost_thread_save_LIBS=$LIBS +boost_thread_save_LDFLAGS=$LDFLAGS +boost_thread_save_CPPFLAGS=$CPPFLAGS # Link-time dependency from thread to system was added as of 1.49.0. if test $boost_major_version -ge 149; then BOOST_SYSTEM([$1]) @@ -902,36 +1064,35 @@ fi # end of the Boost.System check. m4_pattern_allow([^BOOST_SYSTEM_(LIBS|LDFLAGS)$])dnl LIBS="$LIBS $BOOST_SYSTEM_LIBS $boost_cv_pthread_flag" LDFLAGS="$LDFLAGS $BOOST_SYSTEM_LDFLAGS" -# Yes, we *need* to put the -pthread thing in CPPFLAGS because with GCC3, -# boost/thread.hpp will trigger a #error if -pthread isn't used: -# boost/config/requires_threads.hpp:47:5: #error "Compiler threading support -# is not turned on. Please set the correct command line options for -# threading: -pthread (Linux), -pthreads (Solaris) or -mthreads (Mingw32)" CPPFLAGS="$CPPFLAGS $boost_cv_pthread_flag" # When compiling for the Windows platform, the threads library is named -# differently. +# differently. This suffix doesn't exist in new versions of Boost, or +# possibly new versions of GCC on mingw I am assuming it's Boost's change for +# now and I am setting version to 1.48, for lack of knowledge as to when this +# change occurred. +if test $boost_major_version -lt 148; then + case $host_os in + (*mingw*) boost_thread_lib_ext=_win32;; + esac +fi +BOOST_FIND_LIBS([thread], [thread$boost_thread_lib_ext], + [$1], + [boost/thread.hpp], [boost::thread t; boost::mutex m;]) + case $host_os in - (*mingw*) - BOOST_FIND_LIB([thread_win32], [$1], - [boost/thread.hpp], [boost::thread t; boost::mutex m;]) - BOOST_THREAD_LDFLAGS=$BOOST_THREAD_WIN32_LDFLAGS - BOOST_THREAD_LDPATH=$BOOST_THREAD_WIN32_LDPATH - BOOST_THREAD_LIBS=$BOOST_THREAD_WIN32_LIBS - ;; - (*) - BOOST_FIND_LIB([thread], [$1], - [boost/thread.hpp], [boost::thread t; boost::mutex m;]) - ;; + (*mingw*) boost_thread_w32_socket_link=-lws2_32;; esac -BOOST_THREAD_LIBS="$BOOST_THREAD_LIBS $BOOST_SYSTEM_LIBS $boost_cv_pthread_flag" +BOOST_THREAD_LIBS="$BOOST_THREAD_LIBS $BOOST_SYSTEM_LIBS $boost_cv_pthread_flag $boost_thread_w32_socket_link" BOOST_THREAD_LDFLAGS="$BOOST_SYSTEM_LDFLAGS" BOOST_CPPFLAGS="$BOOST_CPPFLAGS $boost_cv_pthread_flag" -LIBS=$boost_threads_save_LIBS -LDFLAGS=$boost_threads_save_LDFLAGS -CPPFLAGS=$boost_threads_save_CPPFLAGS -])# BOOST_THREADS +LIBS=$boost_thread_save_LIBS +LDFLAGS=$boost_thread_save_LDFLAGS +CPPFLAGS=$boost_thread_save_CPPFLAGS +])# BOOST_THREAD + +AU_ALIAS([BOOST_THREADS], [BOOST_THREAD]) # BOOST_TOKENIZER() @@ -979,7 +1140,8 @@ BOOST_DEFUN([Variant], [BOOST_FIND_HEADER([boost/variant/variant_fwd.hpp]) BOOST_FIND_HEADER([boost/variant.hpp])]) -# BOOST_POINTERCONTAINER() + +# BOOST_POINTER_CONTAINER() # ------------------------ # Look for Boost.PointerContainer BOOST_DEFUN([Pointer_Container], @@ -989,12 +1151,13 @@ BOOST_FIND_HEADER([boost/ptr_container/ptr_vector.hpp]) BOOST_FIND_HEADER([boost/ptr_container/ptr_array.hpp]) BOOST_FIND_HEADER([boost/ptr_container/ptr_set.hpp]) BOOST_FIND_HEADER([boost/ptr_container/ptr_map.hpp]) -])# BOOST_POINTERCONTAINER +])# BOOST_POINTER_CONTAINER + # BOOST_WAVE([PREFERRED-RT-OPT]) # ------------------------------ # NOTE: If you intend to use Wave/Spirit with thread support, make sure you -# call BOOST_THREADS first. +# call BOOST_THREAD first. # Look for Boost.Wave. For the documentation of PREFERRED-RT-OPT, see the # documentation of BOOST_FIND_LIB above. BOOST_DEFUN([Wave], @@ -1029,8 +1192,16 @@ BOOST_DEFUN([Xpressive], # _BOOST_PTHREAD_FLAG() # --------------------- -# Internal helper for BOOST_THREADS. Based on ACX_PTHREAD: -# http://autoconf-archive.cryp.to/acx_pthread.html +# Internal helper for BOOST_THREAD. Computes boost_cv_pthread_flag +# which must be used in CPPFLAGS and LIBS. +# +# Yes, we *need* to put the -pthread thing in CPPFLAGS because with GCC3, +# boost/thread.hpp will trigger a #error if -pthread isn't used: +# boost/config/requires_threads.hpp:47:5: #error "Compiler threading support +# is not turned on. Please set the correct command line options for +# threading: -pthread (Linux), -pthreads (Solaris) or -mthreads (Mingw32)" +# +# Based on ACX_PTHREAD: http://autoconf-archive.cryp.to/acx_pthread.html AC_DEFUN([_BOOST_PTHREAD_FLAG], [AC_REQUIRE([AC_PROG_CXX])dnl AC_REQUIRE([AC_CANONICAL_HOST])dnl @@ -1098,6 +1269,14 @@ AC_LANG_POP([C++])dnl m4_define([_BOOST_gcc_test], ["defined __GNUC__ && __GNUC__ == $1 && __GNUC_MINOR__ == $2 && !defined __ICC @ gcc$1$2"])dnl +# _BOOST_mingw_test(MAJOR, MINOR) +# ----------------------------- +# Internal helper for _BOOST_FIND_COMPILER_TAG. +m4_define([_BOOST_mingw_test], +["defined __GNUC__ && __GNUC__ == $1 && __GNUC_MINOR__ == $2 && !defined __ICC && \ + (defined WIN32 || defined WINNT || defined _WIN32 || defined __WIN32 \ + || defined __WIN32__ || defined __WINNT || defined __WINNT__) @ mgw$1$2"])dnl + # _BOOST_FIND_COMPILER_TAG() # -------------------------- @@ -1107,7 +1286,8 @@ m4_define([_BOOST_gcc_test], AC_DEFUN([_BOOST_FIND_COMPILER_TAG], [AC_REQUIRE([AC_PROG_CXX])dnl AC_REQUIRE([AC_CANONICAL_HOST])dnl -AC_CACHE_CHECK([for the toolset name used by Boost for $CXX], [boost_cv_lib_tag], +AC_CACHE_CHECK([for the toolset name used by Boost for $CXX], + [boost_cv_lib_tag], [boost_cv_lib_tag=unknown if test x$boost_cv_inc_path != xno; then AC_LANG_PUSH([C++])dnl @@ -1125,14 +1305,33 @@ if test x$boost_cv_inc_path != xno; then # I'm not sure about my test for `il' (be careful: Intel's ICC pre-defines # the same defines as GCC's). for i in \ + _BOOST_mingw_test(5, 2) \ + _BOOST_gcc_test(5, 2) \ + _BOOST_mingw_test(5, 1) \ + _BOOST_gcc_test(5, 1) \ + _BOOST_mingw_test(5, 0) \ + _BOOST_gcc_test(5, 0) \ + _BOOST_mingw_test(4, 10) \ + _BOOST_gcc_test(4, 10) \ + _BOOST_mingw_test(4, 9) \ + _BOOST_gcc_test(4, 9) \ + _BOOST_mingw_test(4, 8) \ _BOOST_gcc_test(4, 8) \ + _BOOST_mingw_test(4, 7) \ _BOOST_gcc_test(4, 7) \ + _BOOST_mingw_test(4, 6) \ _BOOST_gcc_test(4, 6) \ + _BOOST_mingw_test(4, 5) \ _BOOST_gcc_test(4, 5) \ + _BOOST_mingw_test(4, 4) \ _BOOST_gcc_test(4, 4) \ + _BOOST_mingw_test(4, 3) \ _BOOST_gcc_test(4, 3) \ + _BOOST_mingw_test(4, 2) \ _BOOST_gcc_test(4, 2) \ + _BOOST_mingw_test(4, 1) \ _BOOST_gcc_test(4, 1) \ + _BOOST_mingw_test(4, 0) \ _BOOST_gcc_test(4, 0) \ "defined __GNUC__ && __GNUC__ == 3 && !defined __ICC \ && (defined WIN32 || defined WINNT || defined _WIN32 || defined __WIN32 \ @@ -1198,6 +1397,7 @@ fi])dnl end of AC_CACHE_CHECK # Thread) flavors of Boost. Sets boost_guess_use_mt accordingly. AC_DEFUN([_BOOST_GUESS_WHETHER_TO_USE_MT], [# Check whether we do better use `mt' even though we weren't ask to. +AC_LANG_PUSH([C++])dnl AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ #if defined _REENTRANT || defined _MT || defined __MT__ /* use -mt */ @@ -1205,6 +1405,7 @@ AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ # error MT not needed #endif ]])], [boost_guess_use_mt=:], [boost_guess_use_mt=false]) +AC_LANG_POP([C++])dnl ]) # _BOOST_AC_LINK_IFELSE(PROGRAM, [ACTION-IF-TRUE], [ACTION-IF-FALSE]) @@ -1228,11 +1429,11 @@ boost_use_source=: test -f conftest.$ac_objext && ac_ext=$ac_objext && boost_use_source=false && _AS_ECHO_LOG([re-using the existing conftest.$ac_objext]) AS_IF([_AC_DO_STDERR($ac_link) && { - test -z "$ac_[]_AC_LANG_ABBREV[]_werror_flag" || - test ! -s conftest.err + test -z "$ac_[]_AC_LANG_ABBREV[]_werror_flag" || + test ! -s conftest.err } && test -s conftest$ac_exeext && { - test "$cross_compiling" = yes || - $as_executable_p conftest$ac_exeext + test "$cross_compiling" = yes || + $as_executable_p conftest$ac_exeext dnl FIXME: use AS_TEST_X instead when 2.61 is widespread enough. }], [$2], diff --git a/m4/config_summary.m4 b/m4/config_summary.m4 index d3a6d4007..c5009e0fc 100644 --- a/m4/config_summary.m4 +++ b/m4/config_summary.m4 @@ -26,6 +26,7 @@ echo echo Debug mode.................... : $enable_debug echo C++ compiler.................. : $CXX echo C++ compiler flags............ : $CXXFLAGS +echo C Preprocessor flags.......... : $QUESO_CPPFLAGS echo ' ' echo GSL_LIBS...................... : $GSL_LIBS ###echo GRVY DIR...................... : $GRVY_PREFIX diff --git a/src/Makefile.am b/src/Makefile.am index cba4bf411..e82b3f82e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,4 +1,4 @@ -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += -I$(top_builddir)/src/core/inc # For queso.h AM_CPPFLAGS += $(BOOST_CPPFLAGS) diff --git a/src/apps/queso-config.in b/src/apps/queso-config.in new file mode 100644 index 000000000..a1768bf21 --- /dev/null +++ b/src/apps/queso-config.in @@ -0,0 +1,94 @@ +#!/bin/sh + +# +# values substituted from configure +# +host=@host@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ +builddir=@abs_top_builddir@ +has_been_installed=no + +# +# Define the usage() function +# +usage () +{ + echo "usage: $0 --cppflags --cxxflags --include --libs" + echo " $0 --cxx" + echo " $0 --cc" + echo " $0 --version" + echo " $0 --host" + echo " $0 --ldflags" + exit +} + +# +# Need at least one command-line argument +# +if [ "$#" = "0" ] ; then + usage $0 +fi + +# +# Process the command-line arguments, build up +# return_val +# +return_val="" + +while [ "x$1" != "x" ]; do + case "$1" in + "--cxx") + return_val="@CXX@ $return_val" + ;; + + "--cc") + return_val="@CC@ $return_val" + ;; + + "--cppflags") + return_val="@QUESO_CPPFLAGS@ $return_val" + ;; + + "--cxxflags") + return_val="@CXXFLAGS@ $return_val" + ;; + + "--cflags") + return_val="${CFLAGS} $return_val" + ;; + + "--include") + return_val="-I${includedir} $return_val" + ;; + + "--libs") + return_val="-lqueso $return_val" + ;; + + "--ldflags") + return_val="-Wl,-rpath,${libdir} -L${libdir} $return_val" + ;; + + "--version") + return_val="@VERSION@" + ;; + + "--host") + return_val="$host" + ;; + + *) + echo "Unknown argument: $1" + usage $0 + esac + shift +done + +echo $return_val + +# Local Variables: +# mode: shell-script +# End: diff --git a/src/basic/src/ScalarSequence.C b/src/basic/src/ScalarSequence.C index 532ee6ef2..0a5b0c2cb 100644 --- a/src/basic/src/ScalarSequence.C +++ b/src/basic/src/ScalarSequence.C @@ -2593,14 +2593,61 @@ ScalarSequence::subWriteContents( false, // A 'true' causes problems when the user chooses (via options // in the input file) to use just one file for all outputs. filePtrSet)) { - this->subWriteContents(initialPos, - numPos, - *filePtrSet.ofsVar, - fileType); + + if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT || + fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) { + this->subWriteContents(initialPos, + numPos, + *filePtrSet.ofsVar, + fileType); + } +#ifdef QUESO_HAS_HDF5 + else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { + + // Create dataspace + hsize_t dims[1] = { this->subSequenceSize() }; + hid_t dataspace_id = H5Screate_simple(1, dims, dims); + queso_require_greater_equal_msg( + dataspace_id, + 0, + "error create dataspace with id: " << dataspace_id); + + // Create dataset + hid_t dataset_id = H5Dcreate(filePtrSet.h5Var, + "data", + H5T_IEEE_F64LE, + dataspace_id, + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); + + queso_require_greater_equal_msg( + dataset_id, + 0, + "error creating dataset with id: " << dataset_id); + + // Write the dataset + herr_t status = H5Dwrite( + dataset_id, + H5T_NATIVE_DOUBLE, // The type in memory + H5S_ALL, // The dataspace in memory + dataspace_id, // The file dataspace + H5P_DEFAULT, // Xfer property list + &m_seq[0]); + + queso_require_greater_equal_msg( + status, + 0, + "error writing dataset to file with id: " << filePtrSet.h5Var); + + // Clean up + H5Dclose(dataset_id); + H5Sclose(dataspace_id); + } +#endif + m_env.closeFile(filePtrSet,fileType); } - - return; } // -------------------------------------------------- template @@ -2723,13 +2770,12 @@ ScalarSequence::unifiedWriteContents( if (r == 0) { hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE); //std::cout << "In ScalarSequence::unifiedWriteContents(): h5 case, data type created" << std::endl; - hsize_t dimsf[2]; - dimsf[0] = numParams; - dimsf[1] = chainSize; - hid_t dataspace = H5Screate_simple(2, dimsf, NULL); // HDF5_rank = 2 + hsize_t dimsf[1]; + dimsf[0] = chainSize; + hid_t dataspace = H5Screate_simple(1, dimsf, NULL); // HDF5_rank = 2 //std::cout << "In ScalarSequence::unifiedWriteContents(): h5 case, data space created" << std::endl; hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var, - "seq_of_vectors", + "data", datatype, dataspace, H5P_DEFAULT, // Link creation property list @@ -2742,21 +2788,6 @@ ScalarSequence::unifiedWriteContents( iRC = gettimeofday(&timevalBegin,NULL); if (iRC) {}; // just to remove compiler warning - //double* dataOut[numParams]; // avoid compiler warning - std::vector dataOut((size_t) numParams,NULL); - dataOut[0] = (double*) malloc(numParams*chainSize*sizeof(double)); - for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1' - dataOut[i] = dataOut[i-1] + chainSize; // Yes, just 'chainSize', not 'chainSize*sizeof(double)' - } - //std::cout << "In ScalarSequence::unifiedWriteContents(): h5 case, memory allocated" << std::endl; - for (unsigned int j = 0; j < chainSize; ++j) { - T tmpScalar = m_seq[j]; - for (unsigned int i = 0; i < numParams; ++i) { - dataOut[i][j] = tmpScalar; - } - } - //std::cout << "In ScalarSequence::unifiedWriteContents(): h5 case, memory filled" << std::endl; - herr_t status; //std::cout << "\n In ScalarSequence::unifiedWriteContents(), pos 002 \n" << std::endl; status = H5Dwrite(dataset, @@ -2764,7 +2795,7 @@ ScalarSequence::unifiedWriteContents( H5S_ALL, H5S_ALL, H5P_DEFAULT, - (void*) dataOut[0]); + &m_seq[0]); if (status) {}; // just to remove compiler warning //std::cout << "\n In ScalarSequence::unifiedWriteContents(), pos 003 \n" << std::endl; @@ -2791,10 +2822,6 @@ ScalarSequence::unifiedWriteContents( //std::cout << "In ScalarSequence::unifiedWriteContents(): h5 case, data space closed" << std::endl; H5Tclose(datatype); //std::cout << "In ScalarSequence::unifiedWriteContents(): h5 case, data type closed" << std::endl; - //free(dataOut[0]); // related to the changes above for compiler warning - for (unsigned int tmpIndex = 0; tmpIndex < dataOut.size(); tmpIndex++) { - free (dataOut[tmpIndex]); - } } else { queso_error_msg("hdf file type not supported for multiple sub-environments yet"); @@ -3043,7 +3070,7 @@ ScalarSequence::unifiedReadContents( else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { if (r == 0) { hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var, - "seq_of_vectors", + "data", H5P_DEFAULT); // Dataset access property list hid_t datatype = H5Dget_type(dataset); H5T_class_t t_class = H5Tget_class(datatype); diff --git a/src/basic/src/SequenceOfVectors.C b/src/basic/src/SequenceOfVectors.C index a97e2c243..415a89c0b 100644 --- a/src/basic/src/SequenceOfVectors.C +++ b/src/basic/src/SequenceOfVectors.C @@ -1313,20 +1313,96 @@ SequenceOfVectors::subWriteContents( unsigned int initialPos, unsigned int numPos, FilePtrSetStruct& filePtrSet, - const std::string& fileType) const // "m or hdf" + const std::string& fileType) const { - queso_require_msg(filePtrSet.ofsVar, "filePtrSet.ofsVar should not be NULL"); - if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT || fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) { + queso_require_msg(filePtrSet.ofsVar, "filePtrSet.ofsVar should not be NULL"); this->subWriteContents(initialPos, numPos, *filePtrSet.ofsVar, fileType); } +#ifdef QUESO_HAS_HDF5 else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { - queso_error_msg("hdf file type not supported yet"); + + // Check the file is still legit + queso_require_greater_equal_msg( + filePtrSet.h5Var, + 0, + "filePtrSet.h5Var should not be non-negative"); + + // Sanity check extent + queso_require_less_equal_msg( + (initialPos+numPos), + this->subSequenceSize(), + "invalid routine input parameters"); + + unsigned int numParams = m_vectorSpace.dimLocal(); + unsigned int chainSize = this->subSequenceSize(); + hsize_t dims[2] = { chainSize, numParams }; + + // Create dataspace + hid_t dataspace_id = H5Screate_simple(2, dims, dims); + + // Sanity check dataspace + queso_require_greater_equal_msg( + dataspace_id, + 0, + "error creating dataspace with id: " << dataspace_id); + + // Create dataset + hid_t dataset_id = H5Dcreate(filePtrSet.h5Var, + "data", + H5T_IEEE_F64LE, + dataspace_id, + H5P_DEFAULT, + H5P_DEFAULT, + H5P_DEFAULT); + + // Sanity check dataset + queso_require_greater_equal_msg( + dataset_id, + 0, + "error creating dataset with id: " << dataset_id); + + // This is so egregiously awfully badly terribly horrific I want to die. + // + // And, of course, if any of the subsequent H5* sanity check fail we + // throw an exception and leak a metric fuckton of memory. + double * data = (double *)malloc(numParams * chainSize * sizeof(double)); + + for (unsigned int i = 0; i < chainSize; i++) { + V tmpVec(*(m_seq[i])); + for (unsigned int j = 0; j < numParams; j++) { + data[numParams*i+j] = tmpVec[j]; + } + } + + // Write the dataset + herr_t status = H5Dwrite( + dataset_id, + H5T_NATIVE_DOUBLE, // The type in memory + H5S_ALL, // The dataspace in memory + dataspace_id, // The file dataspace + H5P_DEFAULT, // Xfer property list + data); + + // Sanity check the write + queso_require_greater_equal_msg( + status, + 0, + "error writing dataset to file with id: " << filePtrSet.h5Var); + + // Clean up + free(data); + H5Dclose(dataset_id); + H5Sclose(dataspace_id); + + // Should we close the file too? It's unclear. + } +#endif else { queso_error_msg("invalid file type"); } @@ -1545,12 +1621,12 @@ SequenceOfVectors::unifiedWriteContents( hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE); //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, data type created" << std::endl; hsize_t dimsf[2]; - dimsf[0] = numParams; - dimsf[1] = chainSize; + dimsf[0] = chainSize; + dimsf[1] = numParams; hid_t dataspace = H5Screate_simple(2, dimsf, NULL); // HDF5_rank = 2 //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, data space created" << std::endl; hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var, - "seq_of_vectors", + "data", datatype, dataspace, H5P_DEFAULT, // Link creation property list @@ -1563,20 +1639,15 @@ SequenceOfVectors::unifiedWriteContents( iRC = gettimeofday(&timevalBegin,NULL); if (iRC) {}; // just to remover compiler warning - //double* dataOut[numParams]; // avoid compiler warning - std::vector dataOut((size_t) numParams,NULL); - dataOut[0] = (double*) malloc(numParams*chainSize*sizeof(double)); - for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1' - dataOut[i] = dataOut[i-1] + chainSize; // Yes, just 'chainSize', not 'chainSize*sizeof(double)' - } - //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, memory allocated" << std::endl; - for (unsigned int j = 0; j < chainSize; ++j) { - V tmpVec(*(m_seq[j])); - for (unsigned int i = 0; i < numParams; ++i) { - dataOut[i][j] = tmpVec[i]; + double * data; + data = (double *)malloc(numParams * chainSize * sizeof(double)); + + for (unsigned int i = 0; i < chainSize; ++i) { + V tmpVec(*(m_seq[i])); + for (unsigned int j = 0; j < numParams; ++j) { + data[numParams*i+j] = tmpVec[j]; } } - //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, memory filled" << std::endl; herr_t status; status = H5Dwrite(dataset, @@ -1584,7 +1655,7 @@ SequenceOfVectors::unifiedWriteContents( H5S_ALL, H5S_ALL, H5P_DEFAULT, - (void*) dataOut[0]); + data); if (status) {}; // just to remover compiler warning //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, data written" << std::endl; @@ -1609,10 +1680,7 @@ SequenceOfVectors::unifiedWriteContents( //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, data space closed" << std::endl; H5Tclose(datatype); //std::cout << "In SequenceOfVectors::unifiedWriteContents(): h5 case, data type closed" << std::endl; - //free(dataOut[0]); // related to the changes above for compiler warning - for (unsigned int tmpIndex = 0; tmpIndex < dataOut.size(); tmpIndex++) { - free (dataOut[tmpIndex]); - } + free(data); } else { queso_error_msg("hdf file type not supported for multiple sub-environments yet"); @@ -1856,7 +1924,7 @@ SequenceOfVectors::unifiedReadContents( else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { if (r == 0) { hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var, - "seq_of_vectors", + "data", H5P_DEFAULT); // Dataset access property list hid_t datatype = H5Dget_type(dataset); H5T_class_t t_class = H5Tget_class(datatype); diff --git a/src/core/src/Environment.C b/src/core/src/Environment.C index f325a81d8..9c42ee125 100644 --- a/src/core/src/Environment.C +++ b/src/core/src/Environment.C @@ -112,6 +112,10 @@ FilePtrSetStruct::FilePtrSetStruct() : ofsVar(NULL), ifsVar(NULL) +#ifdef QUESO_HAS_HDF5 + , // lol + h5Var(-1) +#endif { } @@ -600,9 +604,22 @@ BaseEnvironment::openOutputFile( filePtrSet.ofsVar = new std::ofstream((baseFileName+"_sub"+this->subIdString()+"."+fileType).c_str(), std::ofstream::out /*| std::ofstream::in*/ | std::ofstream::app); } +#ifdef QUESO_HAS_HDF5 else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { - queso_error_msg("hdf file type not supported yet"); + std::string fullFileName = + baseFileName+"_sub"+this->subIdString()+"."+fileType; + + // Use H5F_ACC_EXCL because not overwriting, so fail on existing file + filePtrSet.h5Var = H5Fcreate(fullFileName.c_str(), + H5F_ACC_EXCL, + H5P_DEFAULT, + H5P_DEFAULT); + + queso_require_greater_equal_msg( + filePtrSet.h5Var, 0, + "error opening file `" << fullFileName << "`"); } +#endif else { queso_error_msg("invalid file type"); } @@ -637,23 +654,14 @@ BaseEnvironment::openOutputFile( } } // only for matlab formats } - if (filePtrSet.ofsVar != NULL) { - if ((m_subDisplayFile) && (this->displayVerbosity() > 10)) { // output debug - *this->subDisplayFile() << "In BaseEnvironment::openOutputFile()" - << ", subId = " << this->subId() - << ": succeeded on opening output file with base name '" << baseFileName << "." << fileType - << "'" - << ", writeOver = " << writeOver - << std::endl; - } - } - else { - std::cerr << "In BaseEnvironment::openOutputFile()" - << ": failed to open output file with base name '" << baseFileName << "." << fileType - << "'" - << std::endl; + + // Check the file actually opened + if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) || + (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) { + queso_require_msg( + (filePtrSet.ofsVar && filePtrSet.ofsVar->is_open()), + "failed to open output file"); } - queso_require_msg((filePtrSet.ofsVar && filePtrSet.ofsVar->is_open()), "failed to open output file"); } else { returnValue = false; diff --git a/src/core/src/GslVector.C b/src/core/src/GslVector.C index bbadf36ee..245f9ac8b 100644 --- a/src/core/src/GslVector.C +++ b/src/core/src/GslVector.C @@ -51,7 +51,6 @@ GslVector::GslVector(const BaseEnvironment& env, const Map& map) // << "\n map.NumMyElements() = " << map.NumMyElements() // << std::endl; - //std::cout << "Leaving GslVector::constructor(1)" << std::endl; } GslVector::GslVector(const BaseEnvironment& env, const Map& map, double value) @@ -1091,7 +1090,6 @@ GslVector::abs() const } return abs_of_this_vec; - } std::ostream& diff --git a/src/stats/src/BetaJointPdf.C b/src/stats/src/BetaJointPdf.C index dcf692819..80c8abc7a 100644 --- a/src/stats/src/BetaJointPdf.C +++ b/src/stats/src/BetaJointPdf.C @@ -68,11 +68,21 @@ BetaJointPdf::actualValue( V* hessianEffect) const { queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input"); + queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for hessianMatrix and hessianEffect calculations"); - queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); + // No need to multiply by exp(m_logOfNormalizationFactor) because 'lnValue()' is called + double value = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect)); - // No need to multiply by exp(m_logOfNormalizationFactor) because 'lnValue()' is called [PDF-05] - return exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect)); + if (gradVector) { + // DM: This transformation may have issues near the boundary. I haven't + // fully explored this yet. + // + // This evaluation is also normalised (if m_normalizationStyle is zero) and + // so the gradient in physical space contains the normalisation constant. + (*gradVector) *= value; + } + + return value; } //-------------------------------------------------- template @@ -84,18 +94,28 @@ BetaJointPdf::lnValue( M* hessianMatrix, V* hessianEffect) const { - queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); + queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); double aux = 0.; double result = 0.; for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) { if (m_normalizationStyle == 0) { - //aux = log(gsl_ran_beta_pdf(domainVector[i],m_alpha[i],m_beta[i])); aux = log(m_env.basicPdfs()->betaPdfActualValue(domainVector[i],m_alpha[i],m_beta[i])); } else { aux = (m_alpha[i]-1.)*log(domainVector[i]) + (m_beta[i]-1.)*log(1.-domainVector[i]); } + + // The log of the beta pdf is + // f(x) = (alpha - 1) * log(x) + (beta - 1) * log(1 - x) + // Therefore + // df/dx (x) = ((alpha - 1) / x) + ((1 - beta) / (1 - x)) + if (gradVector) { + // We're computing grad of log of p which is p' / p + (*gradVector)[i] = (m_alpha[i] - 1.0) / domainVector[i] + + (1.0 - m_beta[i]) / (1.0 - domainVector[i]); + } + if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) { *m_env.subDisplayFile() << "In BetaJointPdf::lnValue()" << ", m_normalizationStyle = " << m_normalizationStyle @@ -107,7 +127,7 @@ BetaJointPdf::lnValue( } result += aux; } - result += m_logOfNormalizationFactor; // [PDF-05] + result += m_logOfNormalizationFactor; return result; } diff --git a/src/stats/src/GaussianJointPdf.C b/src/stats/src/GaussianJointPdf.C index fc7c8561b..2a4f287e3 100644 --- a/src/stats/src/GaussianJointPdf.C +++ b/src/stats/src/GaussianJointPdf.C @@ -139,17 +139,23 @@ GaussianJointPdf::actualValue( queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input"); - queso_require_msg(!(gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); + queso_require_msg(!(hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); double returnValue = 0.; - if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04 + if (this->m_domainSet.contains(domainVector) == false) { + // What should the gradient be here? returnValue = 0.; } else { + // Already normalised (so the gradient will have the normalisation constant + // in it) returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect)); + + if (gradVector) { + (*gradVector) *= returnValue; + } } - //returnValue *= exp(m_logOfNormalizationFactor); // No need, because 'lnValue()' is called right above [PDF-03] if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) { *m_env.subDisplayFile() << "Leaving GaussianJointPdf::actualValue()" @@ -180,41 +186,72 @@ GaussianJointPdf::lnValue( << std::endl; } - queso_require_msg(!(gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); - - if (domainDirection) {}; // just to remove compiler warning + queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); double returnValue = 0.; double lnDeterminant = 0.; - if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04 + if (this->m_domainSet.contains(domainVector) == false) { + // What should the gradient be here? returnValue = -INFINITY; } else { V diffVec(domainVector - this->lawExpVector()); if (m_diagonalCovMatrix) { returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents(); + + // Compute the gradient of log of the pdf. + // The log of a Gaussian pdf is: + // f(x) = - 1/2 (x - \mu)^T \Sigma^{-1} (x - \mu) + // Therefore + // \frac{df}{dx}(x) = - (x - \mu)^T \Sigma^{-1} + // = - (\Sigma^{-1}^T (x - \mu))^T + // = - (\Sigma^{-1} (x - \mu))^T + // = - \Sigma^{-1} (x - \mu) (row/column vector doesn't matter) + // + // So if \Sigma is diagonal we have a component-wise product of two + // vectors (x - \mu) and the diagonal elements of \Sigma^{-1} + if (gradVector) { + (*gradVector) = diffVec; // Copy + (*gradVector) /= this->lawVarVector(); + (*gradVector) *= -1.0; + } + if (m_normalizationStyle == 0) { unsigned int iMax = this->lawVarVector().sizeLocal(); for (unsigned int i = 0; i < iMax; ++i) { - lnDeterminant += log(this->lawVarVector()[i]); + lnDeterminant += std::log(this->lawVarVector()[i]); } } } else { V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec); returnValue = (diffVec*tmpVec).sumOfComponents(); + + // Compute the gradient of log of the pdf. + // The log of a Gaussian pdf is: + // f(x) = - 1/2 (x - \mu)^T \Sigma^{-1} (x - \mu) + // Therefore + // \frac{df}{dx}(x) = - (x - \mu)^T \Sigma^{-1} + // = - (\Sigma^{-1}^T (x - \mu))^T + // = - (\Sigma^{-1} (x - \mu))^T + // = - \Sigma^{-1} (x - \mu) (row/column vector doesn't matter) + if (gradVector) { + (*gradVector) = tmpVec; + (*gradVector) *= -1.0; + } + if (m_normalizationStyle == 0) { lnDeterminant = this->m_lawCovMatrix->lnDeterminant(); } } if (m_normalizationStyle == 0) { - returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2*M_PI); // normalization of pdf + returnValue += ((double) this->lawVarVector().sizeLocal()) * std::log(2*M_PI); // normalization of pdf returnValue += lnDeterminant; // normalization of pdf } returnValue *= -0.5; } - returnValue += m_logOfNormalizationFactor; // [PDF-03] + returnValue += m_logOfNormalizationFactor; if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) { *m_env.subDisplayFile() << "Leaving GaussianJointPdf::lnValue()" diff --git a/src/stats/src/JeffreysJointPdf.C b/src/stats/src/JeffreysJointPdf.C index 6cabae1d7..d0d169866 100644 --- a/src/stats/src/JeffreysJointPdf.C +++ b/src/stats/src/JeffreysJointPdf.C @@ -126,7 +126,7 @@ JeffreysJointPdf::lnValue( } else { pdf = pdf * (1.0 / domainVector[i]); - result = -log(pdf); + result = std::log(pdf); } } if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) { diff --git a/src/stats/src/LogNormalJointPdf.C b/src/stats/src/LogNormalJointPdf.C index 4e8459a9a..148356da5 100644 --- a/src/stats/src/LogNormalJointPdf.C +++ b/src/stats/src/LogNormalJointPdf.C @@ -103,22 +103,28 @@ LogNormalJointPdf::actualValue( queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input"); - queso_require_msg(!(gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); + queso_require_msg(!(hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); double returnValue = 0.; V zeroVector(domainVector); zeroVector.cwSet(0.); if (domainVector.atLeastOneComponentSmallerOrEqualThan(zeroVector)) { + // What should the gradient be here? returnValue = 0.; } - else if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04 + else if (this->m_domainSet.contains(domainVector) == false) { + // What should the gradient be here? returnValue = 0.; } else { + // Already normalised returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect)); + + if (gradVector) { + (*gradVector) *= returnValue; + } } - //returnValue *= exp(m_logOfNormalizationFactor); // No need, because 'lnValue()' is called right above // [PDF-10] if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) { *m_env.subDisplayFile() << "Leaving LogNormalJointPdf::actualValue()" @@ -147,18 +153,18 @@ LogNormalJointPdf::lnValue( << std::endl; } - queso_require_msg(!(gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); - - if (domainDirection) {}; // just to remove compiler warning + queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations"); double returnValue = 0.; V zeroVector(domainVector); zeroVector.cwSet(0.); if (domainVector.atLeastOneComponentSmallerOrEqualThan(zeroVector)) { + // What should the gradient be here? returnValue = -INFINITY; } - else if (this->m_domainSet.contains(domainVector) == false) { // prudenci 2011-Oct-04 + else if (this->m_domainSet.contains(domainVector) == false) { + // What should the gradient be here? returnValue = -INFINITY; } else { @@ -166,9 +172,20 @@ LogNormalJointPdf::lnValue( V diffVec(zeroVector); for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) { diffVec[i] = std::log(domainVector[i]) - this->lawExpVector()[i]; + + // Compute the gradient of log of the PDF + // The log of a log normal pdf is: + // f(x) = -log(x \sigma sqrt(2 \pi)) - ((log(x) - \mu)^2 / (2 \sigma^2)) + // Therefore + // \frac{df}{dx}(x) = -1/x - (log(x) - \mu) / (x \sigma^2) + if (gradVector) { + (*gradVector)[i] = -(1.0 / domainVector[i]) - + diffVec[i] / (domainVector[i] * this->lawVarVector()[i]); + } } returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents(); returnValue *= -0.5; + if (m_normalizationStyle == 0) { for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) { returnValue -= std::log(domainVector[i] * std::sqrt(2. * M_PI * this->lawVarVector()[i])); // Contribution of 1/(x\sqrt{2\pi\sigma^2}) @@ -178,7 +195,7 @@ LogNormalJointPdf::lnValue( else { queso_error_msg("situation with a non-diagonal covariance matrix makes no sense"); } - returnValue += m_logOfNormalizationFactor; // [PDF-10] + returnValue += m_logOfNormalizationFactor; } if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) { diff --git a/src/stats/src/MetropolisHastingsSG.C b/src/stats/src/MetropolisHastingsSG.C index 001de3b95..614c34023 100644 --- a/src/stats/src/MetropolisHastingsSG.C +++ b/src/stats/src/MetropolisHastingsSG.C @@ -1215,8 +1215,8 @@ MetropolisHastingsSG::generateSequence( m_optionsObj->m_rawChainDataOutputFileType); } - // Compute raw unified MLE - if (workingLogLikelihoodValues) { + // Compute raw unified MLE only in inter0Comm + if (workingLogLikelihoodValues && (m_env.subRank() == 0)) { SequenceOfVectors rawUnifiedMLEpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawUnifiedMLEseq"); double rawUnifiedMLEvalue = workingChain.unifiedPositionsOfMaximum(*workingLogLikelihoodValues, @@ -1240,8 +1240,8 @@ MetropolisHastingsSG::generateSequence( } } - // Compute raw unified MAP - if (workingLogTargetValues) { + // Compute raw unified MAP only in inter0Comm + if (workingLogTargetValues && (m_env.subRank() == 0)) { SequenceOfVectors rawUnifiedMAPpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawUnifiedMAPseq"); double rawUnifiedMAPvalue = workingChain.unifiedPositionsOfMaximum(*workingLogTargetValues, rawUnifiedMAPpositions); diff --git a/src/stats/src/StatisticalForwardProblem.C b/src/stats/src/StatisticalForwardProblem.C index 6b835b751..5b18d7d04 100644 --- a/src/stats/src/StatisticalForwardProblem.C +++ b/src/stats/src/StatisticalForwardProblem.C @@ -279,17 +279,21 @@ StatisticalForwardProblem::solveWithMonteCarlo( << ": instantiating cov and corr matrices" << std::endl; } - pqCovarianceMatrix = new P_M(m_env, - m_paramRv.imageSet().vectorSpace().map(), // number of rows - m_qoiRv.imageSet().vectorSpace().dimGlobal()); // number of cols - pqCorrelationMatrix = new P_M(m_env, - m_paramRv.imageSet().vectorSpace().map(), // number of rows - m_qoiRv.imageSet().vectorSpace().dimGlobal()); // number of cols - ComputeCovCorrMatricesBetweenVectorSequences(*m_paramChain, - *m_qoiChain, - std::min(m_paramRv.realizer().subPeriod(),m_qoiRv.realizer().subPeriod()), // FIX ME: might be INFINITY - *pqCovarianceMatrix, - *pqCorrelationMatrix); + + // Only compute correlations on the inter0Comm communicator + if (m_env.subRank() == 0) { + pqCovarianceMatrix = new P_M(m_env, + m_paramRv.imageSet().vectorSpace().map(), // number of rows + m_qoiRv.imageSet().vectorSpace().dimGlobal()); // number of cols + pqCorrelationMatrix = new P_M(m_env, + m_paramRv.imageSet().vectorSpace().map(), // number of rows + m_qoiRv.imageSet().vectorSpace().dimGlobal()); // number of cols + ComputeCovCorrMatricesBetweenVectorSequences(*m_paramChain, + *m_qoiChain, + std::min(m_paramRv.realizer().subPeriod(),m_qoiRv.realizer().subPeriod()), // FIX ME: might be INFINITY + *pqCovarianceMatrix, + *pqCorrelationMatrix); + } } // Write data out diff --git a/test/Makefile.am b/test/Makefile.am index 7d99d95f8..de23f7101 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -54,10 +54,16 @@ check_PROGRAMS += test_NoInputFile check_PROGRAMS += test_optimizer_options check_PROGRAMS += test_SharedPtr check_PROGRAMS += test_serialEnv +check_PROGRAMS += test_jeffreys_pdf +check_PROGRAMS += test_gaussian_pdf_gradient +check_PROGRAMS += test_log_normal_pdf_gradient +check_PROGRAMS += test_beta_pdf_gradient +check_PROGRAMS += test_intercomm0_gravity +check_PROGRAMS += test_seq_of_vec_hdf5_write LDADD = $(top_builddir)/src/libqueso.la -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += $(BOOST_CPPFLAGS) $(GSL_CFLAGS) $(ANN_CFLAGS) @@ -133,6 +139,20 @@ test_NoInputFile_SOURCES = test_StatisticalInverseProblem/test_NoInputFile.C test_optimizer_options_SOURCES = test_optimizer/test_optimizer_options.C test_SharedPtr_SOURCES = pointers/test_SharedPtr.C test_serialEnv_SOURCES = test_Environment/test_serialEnv.C +test_jeffreys_pdf_SOURCES = test_pdfs/test_jeffreys_pdf.C +test_gaussian_pdf_gradient_SOURCES = test_pdfs/test_gaussian_pdf_gradient.C +test_log_normal_pdf_gradient_SOURCES = test_pdfs/test_log_normal_pdf_gradient.C +test_beta_pdf_gradient_SOURCES = test_pdfs/test_beta_pdf_gradient.C + +test_intercomm0_gravity_SOURCES = +test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_likelihood.C +test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_likelihood.h +test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_main.C +test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_qoi.C +test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_qoi.h +test_intercomm0_gravity_CPPFLAGS = -I$(srcdir)/test_intercomm0 $(AM_CPPFLAGS) + +test_seq_of_vec_hdf5_write_SOURCES = test_SequenceOfVectors/test_HDF5Write.C # Files to freedom stamp srcstamp = @@ -184,6 +204,12 @@ srcstamp += $(test_NoInputFile_SOURCES) srcstamp += $(test_optimizer_options_SOURCES) srcstamp += $(test_SharedPtr_SOURCES) srcstamp += $(test_serialEnv_SOURCES) +srcstamp += $(test_jeffreys_pdf_SOURCES) +srcstamp += $(test_gaussian_pdf_gradient_SOURCES) +srcstamp += $(test_log_normal_pdf_gradient_SOURCES) +srcstamp += $(test_beta_pdf_gradient_SOURCES) +srcstamp += $(test_intercomm0_gravity_SOURCES) +srcstamp += $(test_seq_of_vec_hdf5_write_SOURCES) TESTS = TESTS += $(top_builddir)/test/test_uqEnvironmentNonFatal @@ -235,6 +261,12 @@ TESTS += $(top_builddir)/test/test_NoInputFile TESTS += $(top_builddir)/test/test_optimizer_options TESTS += $(top_builddir)/test/test_SharedPtr TESTS += $(top_builddir)/test/test_serialEnv +TESTS += $(top_builddir)/test/test_jeffreys_pdf +TESTS += $(top_builddir)/test/test_gaussian_pdf_gradient +TESTS += $(top_builddir)/test/test_log_normal_pdf_gradient +TESTS += $(top_builddir)/test/test_beta_pdf_gradient +TESTS += $(top_builddir)/test/test_intercomm0/test_intercomm0_gravity_run.sh +TESTS += $(top_builddir)/test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED @@ -261,6 +293,7 @@ EXTRA_DIST += test_Regression/jeffreys_input.txt.in EXTRA_DIST += test_Regression/test_jeffreys_samples_diff.sh.in EXTRA_DIST += test_Regression/test_jeffreys_samples.m.in EXTRA_DIST += test_Regression/adaptedcov_input.txt.in +EXTRA_DIST += test_gaussian_likelihoods/gaussian_consistency_input.txt.in EXTRA_DIST += test_gaussian_likelihoods/queso_input.txt.in EXTRA_DIST += test_InterpolationSurrogate/queso_input.txt.in EXTRA_DIST += test_SequenceOfVectors/test_unifiedPositionsOfMaximum.sh.in @@ -269,6 +302,10 @@ EXTRA_DIST += test_InputOptionsParser/test_options_bad.txt.in EXTRA_DIST += test_InputOptionsParser/test_options_default.txt.in EXTRA_DIST += test_optimizer/input_test_optimizer_options.in EXTRA_DIST += test_Environment/input_test_serialEnv.in +EXTRA_DIST += test_intercomm0/gravity_1proc.txt.in +EXTRA_DIST += test_intercomm0/gravity_2proc.txt.in +EXTRA_DIST += test_intercomm0/test_intercomm0_gravity_run.sh.in +EXTRA_DIST += test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh.in CLEANFILES = CLEANFILES += test_Environment/debug_output_sub0.txt @@ -286,6 +323,8 @@ clean-local: rm -rf $(top_builddir)/test/test_output_interp_surrogates rm -rf $(top_builddir)/test/output_test_optimizer_options rm -rf $(top_builddir)/test/output_test_serialEnv + rm -rf $(top_builddir)/test/output_test_intercomm0_gravity_1 + rm -rf $(top_builddir)/test/output_test_intercomm0_gravity_2 if CODE_COVERAGE_ENABLED CLEANFILES += *.gcda *.gcno diff --git a/test/gsl_tests/Makefile.am b/test/gsl_tests/Makefile.am index e7ba2eb14..53071693c 100644 --- a/test/gsl_tests/Makefile.am +++ b/test/gsl_tests/Makefile.am @@ -9,7 +9,7 @@ check_PROGRAMS += multiple_rhs_matrix_solve check_PROGRAMS += power_method check_PROGRAMS += inverse_power_method -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += $(BOOST_CPPFLAGS) $(GSL_CFLAGS) $(ANN_CFLAGS) diff --git a/test/gsl_tests/input.in b/test/gsl_tests/input.in index 8458509ca..41d207208 100644 --- a/test/gsl_tests/input.in +++ b/test/gsl_tests/input.in @@ -2,10 +2,10 @@ # UQ Environment ############################################### #env_help = anything -env_numSubEnvironments = 1 -env_subDisplayFileName = outputData/display +env_numSubEnvironments = 1 +env_subDisplayFileName = . env_subDisplayAllowAll = 0 env_subDisplayAllowedSet = 0 -env_displayVerbosity = 0 +env_displayVerbosity = 0 env_syncVerbosity = 0 env_seed = -1 diff --git a/test/t01_valid_cycle/Makefile.am b/test/t01_valid_cycle/Makefile.am index c74875cb8..0797bdc49 100644 --- a/test/t01_valid_cycle/Makefile.am +++ b/test/t01_valid_cycle/Makefile.am @@ -2,7 +2,7 @@ if GIT_CLONE BUILT_SOURCES = .license.stamp endif -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += $(BOOST_CPPFLAGS) $(GSL_CFLAGS) $(ANN_CFLAGS) diff --git a/test/t02_sip_sfp/Makefile.am b/test/t02_sip_sfp/Makefile.am index a1f3fd20a..6c97464c0 100644 --- a/test/t02_sip_sfp/Makefile.am +++ b/test/t02_sip_sfp/Makefile.am @@ -2,7 +2,7 @@ if GIT_CLONE BUILT_SOURCES = .license.stamp endif -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += $(BOOST_CPPFLAGS) $(GSL_CFLAGS) $(ANN_CFLAGS) diff --git a/test/t03_sequence/Makefile.am b/test/t03_sequence/Makefile.am index ec5c5010b..44ec5b603 100644 --- a/test/t03_sequence/Makefile.am +++ b/test/t03_sequence/Makefile.am @@ -2,7 +2,7 @@ if GIT_CLONE BUILT_SOURCES = .license.stamp endif -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += $(BOOST_CPPFLAGS) $(GSL_CFLAGS) $(ANN_CFLAGS) diff --git a/test/t03_sequence/sequence/Makefile.am b/test/t03_sequence/sequence/Makefile.am index b36fa8f2b..36e6cdf37 100644 --- a/test/t03_sequence/sequence/Makefile.am +++ b/test/t03_sequence/sequence/Makefile.am @@ -2,7 +2,7 @@ if SVN_CHECKOUT BUILT_SOURCES = .license.stamp endif -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += -I$(TRILINOS_INCLUDE) diff --git a/test/t04_bimodal/Makefile.am b/test/t04_bimodal/Makefile.am index de4f1a78e..e7c9af064 100644 --- a/test/t04_bimodal/Makefile.am +++ b/test/t04_bimodal/Makefile.am @@ -2,7 +2,7 @@ if GIT_CLONE BUILT_SOURCES = .license.stamp endif -AM_CPPFLAGS = +AM_CPPFLAGS = $(QUESO_CPPFLAGS) AM_CPPFLAGS += -I. AM_CPPFLAGS += -I$(top_builddir)/inc AM_CPPFLAGS += $(BOOST_CPPFLAGS) $(GSL_CFLAGS) $(ANN_CFLAGS) diff --git a/test/test_SequenceOfVectors/test_HDF5Write.C b/test/test_SequenceOfVectors/test_HDF5Write.C new file mode 100644 index 000000000..ab337c882 --- /dev/null +++ b/test/test_SequenceOfVectors/test_HDF5Write.C @@ -0,0 +1,66 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char **argv) { +#ifndef QUESO_HAS_HDF5 + // If we don't have HDF5, skip this test + return 77; +#else +#ifdef QUESO_HAS_MPI + MPI_Init(&argc, &argv); +#endif + + QUESO::EnvOptionsValues options; + options.m_numSubEnvironments = 1; + options.m_subDisplayFileName = "."; + options.m_subDisplayAllowAll = 0; + options.m_subDisplayAllowedSet.insert(0); + options.m_seed = 1.0; + options.m_checkingLevel = 1; + options.m_displayVerbosity = 0; + +#ifdef QUESO_HAS_MPI + QUESO::FullEnvironment env(MPI_COMM_WORLD, "", "", &options); +#else + QUESO::FullEnvironment env("", "", &options); +#endif + + // Create a vector space + QUESO::VectorSpace<> vec_space(env, "", 2, NULL); + + // Create some things to put in the sequence + QUESO::GslVector v1(vec_space.zeroVector()); + v1[0] = 1.0; + v1[1] = 2.0; + + QUESO::GslVector v2(vec_space.zeroVector()); + v2[0] = 3.0; + v2[1] = 4.0; + + // Create a sequence of vectors + QUESO::SequenceOfVectors<> vec_seq(vec_space, 2, ""); + + vec_seq.setPositionValues(0, v1); + vec_seq.setPositionValues(1, v2); + + // Now write + if (env.fullRank() == 0) { + std::set allowedIds; + allowedIds.insert(0); + vec_seq.subWriteContents(0, 2, "output_test_hdf5", "h5", allowedIds); + } + +#ifdef QUESO_HAS_MPI + MPI_Finalize(); +#endif + return 0; +#endif // QUESO_HAS_HDF +} diff --git a/test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh.in b/test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh.in new file mode 100644 index 000000000..11561de16 --- /dev/null +++ b/test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh.in @@ -0,0 +1,6 @@ +#!/bin/bash +set -eu +set -o pipefail + +../libtool --mode=execute ./test_seq_of_vec_hdf5_write +rm output_test_hdf5_sub0.h5 diff --git a/test/test_infinite/test_operator.C b/test/test_infinite/test_operator.C index a1e4f5a92..bf59fef9a 100644 --- a/test/test_infinite/test_operator.C +++ b/test/test_infinite/test_operator.C @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -71,7 +72,7 @@ int main(int argc, char **argv) for (i = 0; i < builder.num_req_eigenpairs; i++) { eig_sys.get_eigenpair(i); norms[i] = eig_sys.calculate_norm(*eig_sys.solution, - libMeshEnums::L2); + libMesh::SystemNorm(libMeshEnums::L2)); if (abs(norms[i] - 1.0) > TEST_TOL) { return 1; } diff --git a/test/test_intercomm0/gravity_1proc.txt.in b/test/test_intercomm0/gravity_1proc.txt.in new file mode 100644 index 000000000..d9fc7612f --- /dev/null +++ b/test/test_intercomm0/gravity_1proc.txt.in @@ -0,0 +1,73 @@ +############################################### +# UQ Environment +############################################### +env_numSubEnvironments = 1 +env_subDisplayFileName = output_test_intercomm0_gravity_1/display_env +env_subDisplayAllowAll = 0 +env_subDisplayAllowedSet = 0 1 2 3 4 5 6 7 +env_displayVerbosity = 2 +env_seed = -1 + +############################################### +# Statistical inverse problem (ip) +############################################### +ip_computeSolution = 1 +ip_dataOutputFileName = output_test_intercomm0_gravity_1/sip_gravity +ip_dataOutputAllowedSet = 0 1 + +############################################### +# Information for Metropolis-Hastings algorithm +############################################### +ip_mh_dataOutputFileName = output_test_intercomm0_gravity_1/sip_gravity +ip_mh_dataOutputAllowedSet = 0 1 + +ip_mh_rawChain_dataInputFileName = . +ip_mh_rawChain_size = 20000 +ip_mh_rawChain_generateExtra = 0 +ip_mh_rawChain_displayPeriod = 2000 +ip_mh_rawChain_measureRunTimes = 1 +ip_mh_rawChain_dataOutputFileName = output_test_intercomm0_gravity_1/sip_gravity_raw_chain +ip_mh_rawChain_dataOutputAllowedSet = 0 1 2 3 4 5 6 7 + +ip_mh_displayCandidates = 0 +ip_mh_putOutOfBoundsInChain = 0 +ip_mh_dr_maxNumExtraStages = 3 +ip_mh_dr_listOfScalesForExtraStages = 5. 10. 20. +ip_mh_am_initialNonAdaptInterval = 0 +ip_mh_am_adaptInterval = 100 +ip_mh_am_eta = 1.98 #(2.4^2)/d, d is the dimension of the problem +ip_mh_am_epsilon = 1.e-5 + +ip_mh_filteredChain_generate = 1 +ip_mh_filteredChain_discardedPortion = 0. +ip_mh_filteredChain_lag = 20 +ip_mh_filteredChain_dataOutputFileName = output_test_intercomm0_gravity_1/sip_gravity_filtered_chain +ip_mh_filteredChain_dataOutputAllowedSet = 0 1 + +############################################### +# Statistical forward problem (fp) +############################################### +fp_help = anything +fp_computeSolution = 1 +fp_computeCovariances = 1 +fp_computeCorrelations = 1 +fp_dataOutputFileName = output_test_intercomm0_gravity_1/sfp_gravity +fp_dataOutputAllowedSet = 0 1 + +############################################### +# 'fp_': information for Monte Carlo algorithm +############################################### +fp_mc_help = anything +fp_mc_dataOutputFileName = output_test_intercomm0_gravity_1/sfp_gravity +fp_mc_dataOutputAllowedSet = 0 1 + +fp_mc_pseq_dataOutputFileName = output_test_intercomm0_gravity_1/sfp_gravity_p_seq +fp_mc_pseq_dataOutputAllowedSet = 0 1 + +fp_mc_qseq_dataInputFileName = . +fp_mc_qseq_size = 16384 +fp_mc_qseq_displayPeriod = 20000 +fp_mc_qseq_measureRunTimes = 1 +fp_mc_qseq_dataOutputFileName = output_test_intercomm0_gravity_1/sfp_gravity_qoi_seq +fp_mc_qseq_dataOutputAllowedSet = 0 1 + diff --git a/test/test_intercomm0/gravity_2proc.txt.in b/test/test_intercomm0/gravity_2proc.txt.in new file mode 100644 index 000000000..93c019c1d --- /dev/null +++ b/test/test_intercomm0/gravity_2proc.txt.in @@ -0,0 +1,73 @@ +############################################### +# UQ Environment +############################################### +env_numSubEnvironments = 1 +env_subDisplayFileName = output_test_intercomm0_gravity_2/display_env +env_subDisplayAllowAll = 0 +env_subDisplayAllowedSet = 0 1 2 3 4 5 6 7 +env_displayVerbosity = 2 +env_seed = -1 + +############################################### +# Statistical inverse problem (ip) +############################################### +ip_computeSolution = 1 +ip_dataOutputFileName = output_test_intercomm0_gravity_2/sip_gravity +ip_dataOutputAllowedSet = 0 1 + +############################################### +# Information for Metropolis-Hastings algorithm +############################################### +ip_mh_dataOutputFileName = output_test_intercomm0_gravity_2/sip_gravity +ip_mh_dataOutputAllowedSet = 0 1 + +ip_mh_rawChain_dataInputFileName = . +ip_mh_rawChain_size = 20000 +ip_mh_rawChain_generateExtra = 0 +ip_mh_rawChain_displayPeriod = 2000 +ip_mh_rawChain_measureRunTimes = 1 +ip_mh_rawChain_dataOutputFileName = output_test_intercomm0_gravity_2/sip_gravity_raw_chain +ip_mh_rawChain_dataOutputAllowedSet = 0 1 2 3 4 5 6 7 + +ip_mh_displayCandidates = 0 +ip_mh_putOutOfBoundsInChain = 0 +ip_mh_dr_maxNumExtraStages = 3 +ip_mh_dr_listOfScalesForExtraStages = 5. 10. 20. +ip_mh_am_initialNonAdaptInterval = 0 +ip_mh_am_adaptInterval = 100 +ip_mh_am_eta = 1.98 #(2.4^2)/d, d is the dimension of the problem +ip_mh_am_epsilon = 1.e-5 + +ip_mh_filteredChain_generate = 1 +ip_mh_filteredChain_discardedPortion = 0. +ip_mh_filteredChain_lag = 20 +ip_mh_filteredChain_dataOutputFileName = output_test_intercomm0_gravity_2/sip_gravity_filtered_chain +ip_mh_filteredChain_dataOutputAllowedSet = 0 1 + +############################################### +# Statistical forward problem (fp) +############################################### +fp_help = anything +fp_computeSolution = 1 +fp_computeCovariances = 1 +fp_computeCorrelations = 1 +fp_dataOutputFileName = output_test_intercomm0_gravity_2/sfp_gravity +fp_dataOutputAllowedSet = 0 1 + +############################################### +# 'fp_': information for Monte Carlo algorithm +############################################### +fp_mc_help = anything +fp_mc_dataOutputFileName = output_test_intercomm0_gravity_2/sfp_gravity +fp_mc_dataOutputAllowedSet = 0 1 + +fp_mc_pseq_dataOutputFileName = output_test_intercomm0_gravity_2/sfp_gravity_p_seq +fp_mc_pseq_dataOutputAllowedSet = 0 1 + +fp_mc_qseq_dataInputFileName = . +fp_mc_qseq_size = 16384 +fp_mc_qseq_displayPeriod = 20000 +fp_mc_qseq_measureRunTimes = 1 +fp_mc_qseq_dataOutputFileName = output_test_intercomm0_gravity_2/sfp_gravity_qoi_seq +fp_mc_qseq_dataOutputAllowedSet = 0 1 + diff --git a/test/test_intercomm0/gravity_likelihood.C b/test/test_intercomm0/gravity_likelihood.C new file mode 100644 index 000000000..a691fe6d6 --- /dev/null +++ b/test/test_intercomm0/gravity_likelihood.C @@ -0,0 +1,93 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +template +Likelihood::Likelihood(const char * prefix, + const QUESO::VectorSet & domain) + : QUESO::BaseScalarFunction(prefix, domain), + m_heights(0), + m_times(0), + m_stdDevs(0) +{ + double const heights[] = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, + 120, 130, 140}; + + double const times [] = {1.41, 2.14, 2.49, 2.87, 3.22, 3.49, 3.81, 4.07, + 4.32, 4.47, 4.75, 4.99, 5.16, 5.26}; + + double const stdDevs[] = {0.020, 0.120, 0.020, 0.010, 0.030, 0.010, 0.030, + 0.030, 0.030, 0.050, 0.010, 0.040, 0.010, 0.09}; + + std::size_t const n = sizeof(heights) / sizeof(*heights); + m_heights.assign(heights, heights + n); + m_times.assign(times, times + n); + m_stdDevs.assign(stdDevs, stdDevs + n); +} + +template +Likelihood::~Likelihood() +{ + // Deconstruct here +} + +template +double +Likelihood::lnValue(const V & domainVector, const V * domainDirection, + V * gradVector, M * hessianMatrix, V * hessianEffect) const +{ + double g = domainVector[0]; + + double misfitValue = 0.0; + for (unsigned int i = 0; i < m_heights.size(); ++i) { + double modelTime = std::sqrt(2.0 * m_heights[i] / g); + double ratio = (modelTime - m_times[i]) / m_stdDevs[i]; + misfitValue += ratio * ratio; + } + + return -0.5 * misfitValue; +} + +template +double +Likelihood::actualValue(const V & domainVector, + const V * domainDirection, V * gradVector, M * hessianMatrix, + V * hessianEffect) const +{ + return std::exp(this->lnValue(domainVector, domainDirection, gradVector, + hessianMatrix, hessianEffect)); +} + +template class Likelihood; diff --git a/test/test_intercomm0/gravity_likelihood.h b/test/test_intercomm0/gravity_likelihood.h new file mode 100644 index 000000000..467dd4746 --- /dev/null +++ b/test/test_intercomm0/gravity_likelihood.h @@ -0,0 +1,48 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#ifndef QUESO_EXAMPLE_GRAVITY_LIKELIHOOD_H +#define QUESO_EXAMPLE_GRAVITY_LIKELIHOOD_H + +#include +#include + +template +class Likelihood : public QUESO::BaseScalarFunction +{ +public: + Likelihood(const char * prefix, const QUESO::VectorSet & domain); + virtual ~Likelihood(); + virtual double lnValue(const V & domainVector, const V * domainDirection, + V * gradVector, M * hessianMatrix, V * hessianEffect) const; + virtual double actualValue(const V & domainVector, const V * domainDirection, + V * gradVector, M * hessianMatrix, V * hessianEffect) const; + +private: + std::vector m_heights; // heights + std::vector m_times; // times + std::vector m_stdDevs; // uncertainties in time measurements +}; + +#endif diff --git a/test/test_intercomm0/gravity_main.C b/test/test_intercomm0/gravity_main.C new file mode 100644 index 000000000..742a03537 --- /dev/null +++ b/test/test_intercomm0/gravity_main.C @@ -0,0 +1,174 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +/* + * Brief description of this file: + * + * The SIP consists on calibrating the magnitude 'g' of acceleration gravity + * using measurements of the time that it takes for an object in free fall to + * reach the ground from a given height and zero initial velocity. The + * solution of the SIP is the posterior probability density function (PDF) of + * 'g'. + * + * The SFP consists of calculating the maximum distance traveled by an object + * in projectile motion. The posterior PDF of 'g' from the SIP might be used + * as input to the SFP. + * + * The code consists of 7 files: + * - 'gravity_main.C' (this file) + * - 'gravity_likelihood.C' (necessary for the SIP) + * - 'gravity_likelihood.h' + * - 'gravity_qoi.C' (necessary for the SFP) + * - 'gravity_qoi.h' + */ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +int main(int argc, char* argv[]) +{ +#ifndef QUESO_HAS_MPI + // Skip this test if we're not in parallel + return 77; +#else + + MPI_Init(&argc, &argv); + + // Initialize QUESO environment + QUESO::FullEnvironment env(MPI_COMM_WORLD, argv[1], "", NULL); + + //================================================================ + // Statistical inverse problem (SIP): find posterior PDF for 'g' + //================================================================ + + //------------------------------------------------------ + // SIP Step 1 of 6: Instantiate the parameter space + //------------------------------------------------------ + QUESO::VectorSpace<> paramSpace(env, "param_", 1, NULL); + + //------------------------------------------------------ + // SIP Step 2 of 6: Instantiate the parameter domain + //------------------------------------------------------ + QUESO::GslVector paramMinValues(paramSpace.zeroVector()); + QUESO::GslVector paramMaxValues(paramSpace.zeroVector()); + + paramMinValues[0] = 8.; + paramMaxValues[0] = 11.; + + QUESO::BoxSubset<> paramDomain("param_", paramSpace, paramMinValues, + paramMaxValues); + + //------------------------------------------------------ + // SIP Step 3 of 6: Instantiate the likelihood function + // object to be used by QUESO. + //------------------------------------------------------ + Likelihood<> lhood("like_", paramDomain); + + //------------------------------------------------------ + // SIP Step 4 of 6: Define the prior RV + //------------------------------------------------------ + QUESO::UniformVectorRV<> priorRv("prior_", paramDomain); + + //------------------------------------------------------ + // SIP Step 5 of 6: Instantiate the inverse problem + //------------------------------------------------------ + // Extra prefix before the default "rv_" prefix + QUESO::GenericVectorRV<> postRv("post_", paramSpace); + + // No extra prefix before the default "ip_" prefix + QUESO::StatisticalInverseProblem<> ip("", NULL, priorRv, lhood, postRv); + + //------------------------------------------------------ + // SIP Step 6 of 6: Solve the inverse problem, that is, + // set the 'pdf' and the 'realizer' of the posterior RV + //------------------------------------------------------ + QUESO::GslVector paramInitials(paramSpace.zeroVector()); + priorRv.realizer().realization(paramInitials); + + QUESO::GslMatrix proposalCovMatrix(paramSpace.zeroVector()); + proposalCovMatrix(0,0) = std::pow(std::abs(paramInitials[0]) / 20.0, 2.0); + + ip.solveWithBayesMetropolisHastings(NULL, paramInitials, &proposalCovMatrix); + + //================================================================ + // Statistical forward problem (SFP): find the max distance + // traveled by an object in projectile motion; input pdf for 'g' + // is the solution of the SIP above. + //================================================================ + + //------------------------------------------------------ + // SFP Step 1 of 6: Instantiate the parameter *and* qoi spaces. + // SFP input RV = FIP posterior RV, so SFP parameter space + // has been already defined. + //------------------------------------------------------ + QUESO::VectorSpace<> qoiSpace(env, "qoi_", 1, NULL); + + //------------------------------------------------------ + // SFP Step 2 of 6: Instantiate the parameter domain + //------------------------------------------------------ + + // Not necessary because input RV of the SFP = output RV of SIP. + // Thus, the parameter domain has been already defined. + + //------------------------------------------------------ + // SFP Step 3 of 6: Instantiate the qoi object + // to be used by QUESO. + //------------------------------------------------------ + Qoi<> qoi("qoi_", paramDomain, qoiSpace); + + //------------------------------------------------------ + // SFP Step 4 of 6: Define the input RV + //------------------------------------------------------ + + // Not necessary because input RV of SFP = output RV of SIP + // (postRv). + + //------------------------------------------------------ + // SFP Step 5 of 6: Instantiate the forward problem + //------------------------------------------------------ + QUESO::GenericVectorRV<> qoiRv("qoi_", qoiSpace); + + QUESO::StatisticalForwardProblem<> fp("", NULL, postRv, qoi, qoiRv); + + //------------------------------------------------------ + // SFP Step 6 of 6: Solve the forward problem + //------------------------------------------------------ + fp.solveWithMonteCarlo(NULL); + + MPI_Finalize(); + + return 0; +#endif // QUESO_HAS_MPI +} diff --git a/test/test_intercomm0/gravity_qoi.C b/test/test_intercomm0/gravity_qoi.C new file mode 100644 index 000000000..f814a8cb0 --- /dev/null +++ b/test/test_intercomm0/gravity_qoi.C @@ -0,0 +1,79 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +/* + * Brief description of this file: + * + * This file contains the code for the user defined qoi routine. + */ + +#include + +#include +#include +#include + +template +Qoi::Qoi(const char * prefix, + const QUESO::VectorSet & domainSet, + const QUESO::VectorSet & imageSet) + : QUESO::BaseVectorFunction(prefix, domainSet, imageSet), + m_angle(M_PI / 4.0), + m_initialVelocity(5.0), + m_initialHeight(0.0) +{ +} + +template +Qoi::~Qoi() +{ + // Deconstruct here +} + +template +void +Qoi::compute(const P_V & domainVector, + const P_V * domainDirection, + Q_V & imageVector, QUESO::DistArray * gradVectors, + QUESO::DistArray * hessianMatrices, + QUESO::DistArray * hessianEffects) const +{ + if (domainVector.sizeLocal() != 1) { + queso_error_msg("domainVector does not have size 1"); + } + if (imageVector.sizeLocal() != 1) { + queso_error_msg("imageVector does not have size 1"); + } + + double g = domainVector[0]; // Sample of the RV 'gravity acceleration' + double distanceTraveled = 0.0; + double aux = m_initialVelocity * std::sin(m_angle); + distanceTraveled = (m_initialVelocity * std::cos(m_angle) / g) * + (aux + std::sqrt(std::pow(aux, 2) + 2.0 * g * m_initialHeight)); + + imageVector[0] = distanceTraveled; +} + +template class Qoi; diff --git a/test/test_intercomm0/gravity_qoi.h b/test/test_intercomm0/gravity_qoi.h new file mode 100644 index 000000000..a5ad6604a --- /dev/null +++ b/test/test_intercomm0/gravity_qoi.h @@ -0,0 +1,60 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +/* + * Brief description of this file: + * + * This is the header file from gravity_qoi.C. + */ + +#ifndef QUESO_EXAMPLE_GRAVITY_QOI_H +#define QUESO_EXAMPLE_GRAVITY_QOI_H + +#include +#include + +template +class Qoi : public QUESO::BaseVectorFunction +{ +public: + Qoi(const char * prefix, const QUESO::VectorSet & domainSet, + const QUESO::VectorSet & imageSet); + virtual ~Qoi(); + virtual void compute(const P_V & domainVector, const P_V * domainDirection, + Q_V & imageVector, QUESO::DistArray * gradVectors, + QUESO::DistArray * hessianMatrices, + QUESO::DistArray * hessianEffects) const; + + void setAngle(double angle); + void setInitialVelocity(double velocity); + void setInitialHeight(double height); + +private: + double m_angle; + double m_initialVelocity; + double m_initialHeight; +}; + +#endif diff --git a/test/test_intercomm0/test_intercomm0_gravity_run.sh.in b/test/test_intercomm0/test_intercomm0_gravity_run.sh.in new file mode 100644 index 000000000..f2b05c18a --- /dev/null +++ b/test/test_intercomm0/test_intercomm0_gravity_run.sh.in @@ -0,0 +1,18 @@ +#!/bin/bash +set -eu +set -o pipefail + +if (( @HAVE_MPI@ )); then + mpirun -np 1 ../libtool --mode=execute ./test_intercomm0_gravity \ + test_intercomm0/gravity_1proc.txt + + mpirun -np 2 ../libtool --mode=execute ./test_intercomm0_gravity \ + test_intercomm0/gravity_2proc.txt + + for i in output_test_intercomm0_gravity_1/*.m; do + j=`basename $i`; + diff output_test_intercomm0_gravity_1/$j output_test_intercomm0_gravity_2/$j; + done +else + exit 77 +fi diff --git a/test/test_pdfs/test_beta_pdf_gradient.C b/test/test_pdfs/test_beta_pdf_gradient.C new file mode 100644 index 000000000..058e22d71 --- /dev/null +++ b/test/test_pdfs/test_beta_pdf_gradient.C @@ -0,0 +1,80 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include +#include +#include + +#define TOL 1e-14 + +int main(int argc, char ** argv) +{ +#ifdef QUESO_HAS_MPI + MPI_Init(&argc, &argv); + + QUESO::FullEnvironment env(MPI_COMM_WORLD, "", "", NULL); +#else + QUESO::FullEnvironment env("", "", NULL); +#endif + + QUESO::VectorSpace<> paramSpace(env, "param_", 1, NULL); + + QUESO::GslVector paramMins(paramSpace.zeroVector()); + paramMins.cwSet(0.0); + + QUESO::GslVector paramMaxs(paramSpace.zeroVector()); + paramMaxs.cwSet(1.0); + + QUESO::BoxSubset<> paramDomain("param_", paramSpace, paramMins, paramMaxs); + + // We should test the other cases of alpha and beta + QUESO::GslVector alpha(paramSpace.zeroVector()); + alpha[0] = 2.0; + + QUESO::GslVector beta(paramSpace.zeroVector()); + beta[0] = 3.0; + + QUESO::BetaJointPdf<> pdf("", paramDomain, alpha, beta); + + // This point is the mode of the beta distribution defined above + QUESO::GslVector point(paramSpace.zeroVector()); + point[0] = (alpha[0] - 1.0) / (alpha[0] + beta[0] - 2.0); + + QUESO::GslVector lnGradVector(paramSpace.zeroVector()); + QUESO::GslVector gradVector(paramSpace.zeroVector()); + + // Here we test the gradient of the log pdf and the gradient of the pdf are + // zero, as we expect. + double lnPdfValue1 = pdf.lnValue(point, NULL, &lnGradVector, NULL, NULL); + double pdfValue1 = pdf.actualValue(point, NULL, &gradVector, NULL, NULL); + + queso_require_less_equal_msg(std::abs(lnGradVector[0]), TOL, "grad log beta pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(gradVector[0]), TOL, "grad beta pdf values are incorrect"); + +#ifdef QUESO_HAS_MPI + MPI_Finalize(); +#endif + + return 0; +} diff --git a/test/test_pdfs/test_gaussian_pdf_gradient.C b/test/test_pdfs/test_gaussian_pdf_gradient.C new file mode 100644 index 000000000..54aa9b6e9 --- /dev/null +++ b/test/test_pdfs/test_gaussian_pdf_gradient.C @@ -0,0 +1,120 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include +#include +#include +#include + +#define TOL 1e-14 + +int main(int argc, char ** argv) +{ +#ifdef QUESO_HAS_MPI + MPI_Init(&argc, &argv); + + QUESO::FullEnvironment env(MPI_COMM_WORLD, "", "", NULL); +#else + QUESO::FullEnvironment env("", "", NULL); +#endif + + unsigned int dim = 3; + QUESO::VectorSpace<> paramSpace(env, "param_", dim, NULL); + + QUESO::GslVector paramMins(paramSpace.zeroVector()); + paramMins.cwSet(-INFINITY); + + QUESO::GslVector paramMaxs(paramSpace.zeroVector()); + paramMaxs.cwSet(INFINITY); + + QUESO::BoxSubset<> paramDomain("param_", paramSpace, paramMins, paramMaxs); + + QUESO::GslVector mean(paramSpace.zeroVector()); + QUESO::GslMatrix var(paramSpace.zeroVector()); + mean[0] = 2.0; + mean[1] = 3.0; + mean[2] = 4.0; + var(0,0) = 5.0; + var(1,1) = 6.0; + var(2,2) = 7.0; + + // Construct a Gaussian PDF + QUESO::GaussianJointPdf<> pdf("", paramDomain, mean, var); + + // Vectors to store gradient calculations + QUESO::GslVector lnGradVector(paramSpace.zeroVector()); + QUESO::GslVector gradVector(paramSpace.zeroVector()); + + // Where to evaluate the gradient. Evaluating at the mean (the mode for a + // Gaussian) should give a gradient consisting of a vector of zeros. + QUESO::GslVector point(mean); + + // We are testing that the gradient of log of the pdf is all zeros + double lnPdfValue1 = pdf.lnValue(point, NULL, &lnGradVector, NULL, NULL); + queso_require_less_equal_msg(std::abs(lnGradVector[0]), TOL, "grad log gaussian pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(lnGradVector[1]), TOL, "grad log gaussian pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(lnGradVector[2]), TOL, "grad log gaussian pdf values are incorrect"); + + // We are testing that the of the pdf is all zeros + double pdfValue1 = pdf.actualValue(point, NULL, &gradVector, NULL, NULL); + queso_require_less_equal_msg(std::abs(gradVector[0]), TOL, "grad guassian pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(gradVector[1]), TOL, "grad guassian pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(gradVector[2]), TOL, "grad guassian pdf values are incorrect"); + + + // Now construct another Gaussian. This time we're constructing a Gaussian + // that we know will have a gradient consisting entirely of ones (in log + // space). + mean[0] = 0.0; + mean[1] = 0.0; + mean[2] = 0.0; + + var(0,0) = 1.0; + var(0,1) = 0.8; + var(0,2) = 0.7; + var(1,0) = 0.8; + var(1,1) = 2.0; + var(1,2) = 0.6; + var(2,0) = 0.7; + var(2,1) = 0.6; + var(2,2) = 3.0; + + point[0] = -2.5; + point[1] = -3.4; + point[2] = -4.3; + + QUESO::GaussianJointPdf<> pdf2("", paramDomain, mean, var); + + lnPdfValue1 = pdf2.lnValue(point, NULL, &lnGradVector, NULL, NULL); + + queso_require_less_equal_msg(std::abs(lnGradVector[0] - 1.0), TOL, "grad log gaussian pdf2 values are incorrect"); + queso_require_less_equal_msg(std::abs(lnGradVector[1] - 1.0), TOL, "grad log gaussian pdf2 values are incorrect"); + queso_require_less_equal_msg(std::abs(lnGradVector[2] - 1.0), TOL, "grad log gaussian pdf2 values are incorrect"); + +#ifdef QUESO_HAS_MPI + MPI_Finalize(); +#endif + + return 0; +} diff --git a/test/test_pdfs/test_jeffreys_pdf.C b/test/test_pdfs/test_jeffreys_pdf.C new file mode 100644 index 000000000..45a68d8e6 --- /dev/null +++ b/test/test_pdfs/test_jeffreys_pdf.C @@ -0,0 +1,70 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include +#include +#include + +#define TOL 1e-15 + +int main(int argc, char ** argv) +{ +#ifdef QUESO_HAS_MPI + MPI_Init(&argc, &argv); + + QUESO::FullEnvironment env(MPI_COMM_WORLD, "", "", NULL); +#else + QUESO::FullEnvironment env("", "", NULL); +#endif + + QUESO::VectorSpace<> paramSpace(env, "param_", 1, NULL); + + QUESO::GslVector paramMins(paramSpace.zeroVector()); + paramMins.cwSet(0.0); + + QUESO::GslVector paramMaxs(paramSpace.zeroVector()); + paramMaxs.cwSet(INFINITY); + + QUESO::BoxSubset<> paramDomain("param_", paramSpace, paramMins, paramMaxs); + + QUESO::JeffreysJointPdf<> pdf("", paramDomain); + + QUESO::GslVector point(paramSpace.zeroVector()); + point[0] = 2.8; + + double lnPdfValue1 = pdf.lnValue(point, NULL, NULL, NULL, NULL); + double lnPdfValue2 = std::log(1.0 / point[0]); + + double pdfValue1 = pdf.actualValue(point, NULL, NULL, NULL, NULL); + double pdfValue2 = std::exp(lnPdfValue2); + + queso_require_less_equal_msg(std::abs(lnPdfValue1 - lnPdfValue2), TOL, "log jeffreys pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(pdfValue1 - pdfValue2), TOL, "jeffreys pdf values are incorrect"); + +#ifdef QUESO_HAS_MPI + MPI_Finalize(); +#endif + + return 0; +} diff --git a/test/test_pdfs/test_log_normal_pdf_gradient.C b/test/test_pdfs/test_log_normal_pdf_gradient.C new file mode 100644 index 000000000..85c7bcd89 --- /dev/null +++ b/test/test_pdfs/test_log_normal_pdf_gradient.C @@ -0,0 +1,89 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2015 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include +#include +#include + +#define TOL 1e-14 + +int main(int argc, char ** argv) +{ +#ifdef QUESO_HAS_MPI + MPI_Init(&argc, &argv); + + QUESO::FullEnvironment env(MPI_COMM_WORLD, "", "", NULL); +#else + QUESO::FullEnvironment env("", "", NULL); +#endif + + unsigned int dim = 3; + QUESO::VectorSpace<> paramSpace(env, "param_", dim, NULL); + + QUESO::GslVector paramMins(paramSpace.zeroVector()); + paramMins.cwSet(-INFINITY); + + QUESO::GslVector paramMaxs(paramSpace.zeroVector()); + paramMaxs.cwSet(INFINITY); + + QUESO::BoxSubset<> paramDomain("param_", paramSpace, paramMins, paramMaxs); + + QUESO::GslVector mean(paramSpace.zeroVector()); + QUESO::GslVector var(paramSpace.zeroVector()); + mean[0] = 2.0; + mean[1] = 3.0; + mean[2] = 4.0; + var[0] = 5.0; + var[1] = 6.0; + var[2] = 7.0; + + QUESO::LogNormalJointPdf<> pdf("", paramDomain, mean, var); + + // This point is the mode of the log normal distribution defined above + QUESO::GslVector point(paramSpace.zeroVector()); + point[0] = std::exp(mean[0] - var[0]); + point[1] = std::exp(mean[1] - var[1]); + point[2] = std::exp(mean[2] - var[2]); + + QUESO::GslVector lnGradVector(paramSpace.zeroVector()); + QUESO::GslVector gradVector(paramSpace.zeroVector()); + + // Here we test the gradient of the log of the pdf and the gradient of the + // pdf are zero at the mode, as we expect. + double lnPdfValue1 = pdf.lnValue(point, NULL, &lnGradVector, NULL, NULL); + queso_require_less_equal_msg(std::abs(lnGradVector[0]), TOL, "grad log log_normal pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(lnGradVector[1]), TOL, "grad log log_normal pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(lnGradVector[2]), TOL, "grad log log_normal pdf values are incorrect"); + + double pdfValue1 = pdf.actualValue(point, NULL, &gradVector, NULL, NULL); + queso_require_less_equal_msg(std::abs(gradVector[0]), TOL, "grad log_normal pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(gradVector[1]), TOL, "grad log_normal pdf values are incorrect"); + queso_require_less_equal_msg(std::abs(gradVector[2]), TOL, "grad log_normal pdf values are incorrect"); + +#ifdef QUESO_HAS_MPI + MPI_Finalize(); +#endif + + return 0; +}