From 65adc7240ffd908daf89ea007bf19167c345cb1f Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Tue, 1 Dec 2015 16:45:52 -0600 Subject: [PATCH 01/47] Fix compatibility with recent libMesh This is the same bug as was fixed in 25fd4c51 - I just missed this instance because I hadn't been using a SLEPc-enabled libMesh build. --- test/test_infinite/test_operator.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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; } From d2943ffe6364ca6c3f2052263745dac5b5410a62 Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Wed, 16 Dec 2015 13:43:28 -0600 Subject: [PATCH 02/47] Update to newest https://github.com/tsuna/boost.m4 This fixes configuration for me on Ubuntu 15.10 --- m4/common/boost.m4 | 475 ++++++++++++++++++++++++++++++++------------- 1 file changed, 338 insertions(+), 137 deletions(-) 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], From bfbf52f7fa05dd9548f774553badd25018331889 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Thu, 4 Feb 2016 13:14:53 -0600 Subject: [PATCH 03/47] Fix sign error in lnValue of JeffreysJointPdf The Jeffreys prior is 1.0 / x. Log of it is log(1.0 / x) = -log(x), but what was returned was the negative of this. --- src/stats/src/JeffreysJointPdf.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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)) { From 7ea0d45f55a83dfe042f36399ad8d3f62253d9cf Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Thu, 11 Feb 2016 15:28:23 -0600 Subject: [PATCH 04/47] Add jeffreys pdf evaluation test --- .gitignore | 1 + test/Makefile.am | 4 ++ test/test_pdfs/test_jeffreys_pdf.C | 70 ++++++++++++++++++++++++++++++ 3 files changed, 75 insertions(+) create mode 100644 test/test_pdfs/test_jeffreys_pdf.C diff --git a/.gitignore b/.gitignore index ba4be7dd5..38c50aa1e 100644 --- a/.gitignore +++ b/.gitignore @@ -214,3 +214,4 @@ test/test_SharedPtr test/output_test_serialEnv test/test_Environment/input_test_serialEnv test/test_serialEnv +test/test_jeffreys_pdf diff --git a/test/Makefile.am b/test/Makefile.am index 7d99d95f8..0b8a1b1ef 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -54,6 +54,7 @@ check_PROGRAMS += test_NoInputFile check_PROGRAMS += test_optimizer_options check_PROGRAMS += test_SharedPtr check_PROGRAMS += test_serialEnv +check_PROGRAMS += test_jeffreys_pdf LDADD = $(top_builddir)/src/libqueso.la @@ -133,6 +134,7 @@ 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 # Files to freedom stamp srcstamp = @@ -184,6 +186,7 @@ srcstamp += $(test_NoInputFile_SOURCES) srcstamp += $(test_optimizer_options_SOURCES) srcstamp += $(test_SharedPtr_SOURCES) srcstamp += $(test_serialEnv_SOURCES) +srcstamp += $(test_jeffreys_pdf_SOURCES) TESTS = TESTS += $(top_builddir)/test/test_uqEnvironmentNonFatal @@ -235,6 +238,7 @@ 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 XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED 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; +} From b49d938adbb74158cad42d6cad801f808940137b Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Fri, 12 Feb 2016 13:25:18 -0600 Subject: [PATCH 05/47] Add gradients of Gaussian distribution (with tests) --- .gitignore | 1 + src/stats/src/GaussianJointPdf.C | 38 +++++++--- test/Makefile.am | 4 + test/test_pdfs/test_gaussian_pdf_gradient.C | 81 +++++++++++++++++++++ 4 files changed, 114 insertions(+), 10 deletions(-) create mode 100644 test/test_pdfs/test_gaussian_pdf_gradient.C diff --git a/.gitignore b/.gitignore index 38c50aa1e..f7c669464 100644 --- a/.gitignore +++ b/.gitignore @@ -215,3 +215,4 @@ test/output_test_serialEnv test/test_Environment/input_test_serialEnv test/test_serialEnv test/test_jeffreys_pdf +test/test_gaussian_pdf_gradient diff --git a/src/stats/src/GaussianJointPdf.C b/src/stats/src/GaussianJointPdf.C index fc7c8561b..0cfdeefd4 100644 --- a/src/stats/src/GaussianJointPdf.C +++ b/src/stats/src/GaussianJointPdf.C @@ -139,17 +139,22 @@ 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 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 +185,54 @@ 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 + 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(); + + 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/test/Makefile.am b/test/Makefile.am index 0b8a1b1ef..d62cd0203 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -55,6 +55,7 @@ check_PROGRAMS += test_optimizer_options check_PROGRAMS += test_SharedPtr check_PROGRAMS += test_serialEnv check_PROGRAMS += test_jeffreys_pdf +check_PROGRAMS += test_gaussian_pdf_gradient LDADD = $(top_builddir)/src/libqueso.la @@ -135,6 +136,7 @@ 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 # Files to freedom stamp srcstamp = @@ -187,6 +189,7 @@ srcstamp += $(test_optimizer_options_SOURCES) srcstamp += $(test_SharedPtr_SOURCES) srcstamp += $(test_serialEnv_SOURCES) srcstamp += $(test_jeffreys_pdf_SOURCES) +srcstamp += $(test_gaussian_pdf_gradient_SOURCES) TESTS = TESTS += $(top_builddir)/test/test_uqEnvironmentNonFatal @@ -239,6 +242,7 @@ 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 XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED 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..f760c5d2b --- /dev/null +++ b/test/test_pdfs/test_gaussian_pdf_gradient.C @@ -0,0 +1,81 @@ +//-----------------------------------------------------------------------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 + + 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::GaussianJointPdf<> pdf("", paramDomain, mean, var); + + QUESO::GslVector lnGradVector(paramSpace.zeroVector()); + QUESO::GslVector gradVector(paramSpace.zeroVector()); + + double lnPdfValue1 = pdf.lnValue(mean, 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"); + + double pdfValue1 = pdf.actualValue(mean, 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"); + +#ifdef QUESO_HAS_MPI + MPI_Finalize(); +#endif + + return 0; +} From e7da9f6bfc35cf2ac08561c1135ec883d9a146da Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Mon, 15 Feb 2016 14:18:02 -0600 Subject: [PATCH 06/47] Make a comment clearer --- src/stats/src/GaussianJointPdf.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stats/src/GaussianJointPdf.C b/src/stats/src/GaussianJointPdf.C index 0cfdeefd4..5e80b6df5 100644 --- a/src/stats/src/GaussianJointPdf.C +++ b/src/stats/src/GaussianJointPdf.C @@ -148,7 +148,8 @@ GaussianJointPdf::actualValue( returnValue = 0.; } else { - // Already normalised + // Already normalised (so the gradient will have the normalisation constant + // in it) returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect)); if (gradVector) { From 411e7cb5514f3f944fd312fc067223ad1ea2a093 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 16 Feb 2016 14:05:02 -0600 Subject: [PATCH 07/47] Adding another (analytical) gradient test for the Gaussian --- test/test_pdfs/test_gaussian_pdf_gradient.C | 43 ++++++++++++++++++--- 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/test/test_pdfs/test_gaussian_pdf_gradient.C b/test/test_pdfs/test_gaussian_pdf_gradient.C index f760c5d2b..e809ecdda 100644 --- a/test/test_pdfs/test_gaussian_pdf_gradient.C +++ b/test/test_pdfs/test_gaussian_pdf_gradient.C @@ -24,6 +24,7 @@ #include #include +#include #include #define TOL 1e-15 @@ -50,29 +51,59 @@ int main(int argc, char ** argv) QUESO::BoxSubset<> paramDomain("param_", paramSpace, paramMins, paramMaxs); QUESO::GslVector mean(paramSpace.zeroVector()); - QUESO::GslVector var(paramSpace.zeroVector()); + QUESO::GslMatrix 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; + var(0,0) = 5.0; + var(1,1) = 6.0; + var(2,2) = 7.0; + // Check the zero-gradient case QUESO::GaussianJointPdf<> pdf("", paramDomain, mean, var); QUESO::GslVector lnGradVector(paramSpace.zeroVector()); QUESO::GslVector gradVector(paramSpace.zeroVector()); + QUESO::GslVector point(mean); - double lnPdfValue1 = pdf.lnValue(mean, NULL, &lnGradVector, NULL, NULL); + 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"); - double pdfValue1 = pdf.actualValue(mean, NULL, &gradVector, NULL, NULL); + 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"); + + // Check the 1-gradient case (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 From 762da4afc80f05e2a202b5288e43fa9f48f8c279 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Fri, 4 Mar 2016 16:58:03 -0600 Subject: [PATCH 08/47] Adding some mathemtics comments to explain derivative calculation --- src/stats/src/GaussianJointPdf.C | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/stats/src/GaussianJointPdf.C b/src/stats/src/GaussianJointPdf.C index 5e80b6df5..2a4f287e3 100644 --- a/src/stats/src/GaussianJointPdf.C +++ b/src/stats/src/GaussianJointPdf.C @@ -200,7 +200,17 @@ GaussianJointPdf::lnValue( if (m_diagonalCovMatrix) { returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents(); - // Compute the gradient of log of the pdf + // 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(); @@ -218,6 +228,14 @@ GaussianJointPdf::lnValue( 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; From 4cd2da8bbbeec4f3466dcd53b3c987bd513bfd86 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Fri, 4 Mar 2016 16:59:01 -0600 Subject: [PATCH 09/47] Loosen up tolerance for testing on different archs --- test/test_pdfs/test_gaussian_pdf_gradient.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_pdfs/test_gaussian_pdf_gradient.C b/test/test_pdfs/test_gaussian_pdf_gradient.C index e809ecdda..565cc1ab9 100644 --- a/test/test_pdfs/test_gaussian_pdf_gradient.C +++ b/test/test_pdfs/test_gaussian_pdf_gradient.C @@ -27,7 +27,7 @@ #include #include -#define TOL 1e-15 +#define TOL 1e-14 int main(int argc, char ** argv) { From a667005af2f26e669e42d05b78cc1a1723bd4505 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Fri, 4 Mar 2016 17:04:40 -0600 Subject: [PATCH 10/47] Comment the Gaussian gradient test better --- test/test_pdfs/test_gaussian_pdf_gradient.C | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/test/test_pdfs/test_gaussian_pdf_gradient.C b/test/test_pdfs/test_gaussian_pdf_gradient.C index 565cc1ab9..54aa9b6e9 100644 --- a/test/test_pdfs/test_gaussian_pdf_gradient.C +++ b/test/test_pdfs/test_gaussian_pdf_gradient.C @@ -59,25 +59,33 @@ int main(int argc, char ** argv) var(1,1) = 6.0; var(2,2) = 7.0; - // Check the zero-gradient case + // 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"); - // Check the 1-gradient case (in log space) + // 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; From 76942120290a30a7d049bd48c8f7778de801e9cc Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Fri, 12 Feb 2016 17:16:02 -0600 Subject: [PATCH 11/47] Code for gradient of lognormals --- .gitignore | 1 + src/stats/src/LogNormalJointPdf.C | 28 ++++-- test/Makefile.am | 4 + test/test_pdfs/test_log_normal_pdf_gradient.C | 86 +++++++++++++++++++ 4 files changed, 111 insertions(+), 8 deletions(-) create mode 100644 test/test_pdfs/test_log_normal_pdf_gradient.C diff --git a/.gitignore b/.gitignore index f7c669464..647b688e7 100644 --- a/.gitignore +++ b/.gitignore @@ -216,3 +216,4 @@ test/test_Environment/input_test_serialEnv test/test_serialEnv test/test_jeffreys_pdf test/test_gaussian_pdf_gradient +test/test_log_normal_pdf_gradient diff --git a/src/stats/src/LogNormalJointPdf.C b/src/stats/src/LogNormalJointPdf.C index 4e8459a9a..0d6bb4ec0 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,15 @@ LogNormalJointPdf::lnValue( V diffVec(zeroVector); for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) { diffVec[i] = std::log(domainVector[i]) - this->lawExpVector()[i]; + + 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 +190,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/test/Makefile.am b/test/Makefile.am index d62cd0203..a55b3c173 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -56,6 +56,7 @@ 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 LDADD = $(top_builddir)/src/libqueso.la @@ -137,6 +138,7 @@ 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 # Files to freedom stamp srcstamp = @@ -190,6 +192,7 @@ 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) TESTS = TESTS += $(top_builddir)/test/test_uqEnvironmentNonFatal @@ -243,6 +246,7 @@ 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 XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED 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..c199bf378 --- /dev/null +++ b/test/test_pdfs/test_log_normal_pdf_gradient.C @@ -0,0 +1,86 @@ +//-----------------------------------------------------------------------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 + + 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); + + 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()); + + 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; +} From ba420cdc0567cc7b9d760b7678726439f0cb2f41 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sun, 27 Mar 2016 16:41:00 -0400 Subject: [PATCH 12/47] Add some comments explaining grad of lognormal --- src/stats/src/LogNormalJointPdf.C | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/stats/src/LogNormalJointPdf.C b/src/stats/src/LogNormalJointPdf.C index 0d6bb4ec0..148356da5 100644 --- a/src/stats/src/LogNormalJointPdf.C +++ b/src/stats/src/LogNormalJointPdf.C @@ -173,6 +173,11 @@ LogNormalJointPdf::lnValue( 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]); From b8c985097b84d0c7d007661f55236dbbc430c501 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sun, 27 Mar 2016 16:42:24 -0400 Subject: [PATCH 13/47] Loosen up test tolerance --- test/test_pdfs/test_log_normal_pdf_gradient.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_pdfs/test_log_normal_pdf_gradient.C b/test/test_pdfs/test_log_normal_pdf_gradient.C index c199bf378..188f0e5f4 100644 --- a/test/test_pdfs/test_log_normal_pdf_gradient.C +++ b/test/test_pdfs/test_log_normal_pdf_gradient.C @@ -26,7 +26,7 @@ #include #include -#define TOL 1e-15 +#define TOL 1e-14 int main(int argc, char ** argv) { From fc5daa24cbe720795fbfaf1141a15c14a0da4ddc Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sun, 27 Mar 2016 16:44:48 -0400 Subject: [PATCH 14/47] Add comments describing what the test does --- test/test_pdfs/test_log_normal_pdf_gradient.C | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_pdfs/test_log_normal_pdf_gradient.C b/test/test_pdfs/test_log_normal_pdf_gradient.C index 188f0e5f4..85c7bcd89 100644 --- a/test/test_pdfs/test_log_normal_pdf_gradient.C +++ b/test/test_pdfs/test_log_normal_pdf_gradient.C @@ -60,6 +60,7 @@ int main(int argc, char ** argv) 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]); @@ -68,6 +69,8 @@ int main(int argc, char ** argv) 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"); From c962a9f6d206fc8f248c30b9e357f60a5407c489 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Thu, 11 Feb 2016 18:05:48 -0600 Subject: [PATCH 15/47] Add gradient code for beta distribution --- src/stats/src/BetaJointPdf.C | 27 +++++++-- test/Makefile.am | 4 ++ test/test_pdfs/test_beta_pdf_gradient.C | 77 +++++++++++++++++++++++++ 3 files changed, 103 insertions(+), 5 deletions(-) create mode 100644 test/test_pdfs/test_beta_pdf_gradient.C diff --git a/src/stats/src/BetaJointPdf.C b/src/stats/src/BetaJointPdf.C index dcf692819..5c4f1be3b 100644 --- a/src/stats/src/BetaJointPdf.C +++ b/src/stats/src/BetaJointPdf.C @@ -68,11 +68,16 @@ 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) { + (*gradVector) *= value; + } + + return value; } //-------------------------------------------------- template @@ -84,18 +89,30 @@ 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]); } + + if (gradVector) { + (*gradVector)[i] = (m_alpha[i] - 1.0) * + std::pow(domainVector[i], m_alpha[i] - 2.0) * + std::pow(1.0 - domainVector[i], m_beta[i] - 1.0) + + (1.0 - m_beta[i]) * + std::pow(domainVector[i], m_alpha[i] - 1.0) * + std::pow(1.0 - domainVector[i], m_beta[i] - 2.0); + + // We're computing grad of log of p which is p' / p + (*gradVector)[i] = (*gradVector)[i] / std::exp(aux); + } + if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) { *m_env.subDisplayFile() << "In BetaJointPdf::lnValue()" << ", m_normalizationStyle = " << m_normalizationStyle diff --git a/test/Makefile.am b/test/Makefile.am index a55b3c173..298347908 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -57,6 +57,7 @@ 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 LDADD = $(top_builddir)/src/libqueso.la @@ -139,6 +140,7 @@ 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 # Files to freedom stamp srcstamp = @@ -193,6 +195,7 @@ 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) TESTS = TESTS += $(top_builddir)/test/test_uqEnvironmentNonFatal @@ -247,6 +250,7 @@ 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 XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED 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..b4b9d0740 --- /dev/null +++ b/test/test_pdfs/test_beta_pdf_gradient.C @@ -0,0 +1,77 @@ +//-----------------------------------------------------------------------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(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); + + 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()); + + 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; +} From 3d8bdac806d419d24922a4bdbd6cbdbceaf8e791 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Mon, 15 Feb 2016 13:45:27 -0600 Subject: [PATCH 16/47] Use already-cancelled gradient formula in log space --- src/stats/src/BetaJointPdf.C | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/stats/src/BetaJointPdf.C b/src/stats/src/BetaJointPdf.C index 5c4f1be3b..5545faf16 100644 --- a/src/stats/src/BetaJointPdf.C +++ b/src/stats/src/BetaJointPdf.C @@ -74,6 +74,11 @@ BetaJointPdf::actualValue( double value = std::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; } @@ -102,15 +107,9 @@ BetaJointPdf::lnValue( } if (gradVector) { - (*gradVector)[i] = (m_alpha[i] - 1.0) * - std::pow(domainVector[i], m_alpha[i] - 2.0) * - std::pow(1.0 - domainVector[i], m_beta[i] - 1.0) + - (1.0 - m_beta[i]) * - std::pow(domainVector[i], m_alpha[i] - 1.0) * - std::pow(1.0 - domainVector[i], m_beta[i] - 2.0); - // We're computing grad of log of p which is p' / p - (*gradVector)[i] = (*gradVector)[i] / std::exp(aux); + (*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)) { From c8230ab87018de0ceb6a924e3e1a0d33770e9f9c Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Wed, 17 Feb 2016 11:29:30 -0600 Subject: [PATCH 17/47] Update ignores --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 647b688e7..14f25904a 100644 --- a/.gitignore +++ b/.gitignore @@ -217,3 +217,4 @@ test/test_serialEnv test/test_jeffreys_pdf test/test_gaussian_pdf_gradient test/test_log_normal_pdf_gradient +test/test_beta_pdf_gradient From ccc3da95434b97c44918deb571fb303e46ae22c7 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sun, 27 Mar 2016 16:47:53 -0400 Subject: [PATCH 18/47] Add comment detailing derivative of log beta pdf --- src/stats/src/BetaJointPdf.C | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/stats/src/BetaJointPdf.C b/src/stats/src/BetaJointPdf.C index 5545faf16..80c8abc7a 100644 --- a/src/stats/src/BetaJointPdf.C +++ b/src/stats/src/BetaJointPdf.C @@ -106,6 +106,10 @@ BetaJointPdf::lnValue( 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] + @@ -123,7 +127,7 @@ BetaJointPdf::lnValue( } result += aux; } - result += m_logOfNormalizationFactor; // [PDF-05] + result += m_logOfNormalizationFactor; return result; } From ec8dbb7ba9c2cfbeff43069ce49a43507ee6448f Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sun, 27 Mar 2016 16:48:15 -0400 Subject: [PATCH 19/47] Loosen up test tolerance --- test/test_pdfs/test_beta_pdf_gradient.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_pdfs/test_beta_pdf_gradient.C b/test/test_pdfs/test_beta_pdf_gradient.C index b4b9d0740..4bc40f9c3 100644 --- a/test/test_pdfs/test_beta_pdf_gradient.C +++ b/test/test_pdfs/test_beta_pdf_gradient.C @@ -26,7 +26,7 @@ #include #include -#define TOL 1e-15 +#define TOL 1e-14 int main(int argc, char ** argv) { From eca8c5c96b43746ec6aae37d4e4fe477216d13af Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sun, 27 Mar 2016 16:49:35 -0400 Subject: [PATCH 20/47] Add some comments explaining what the test does --- test/test_pdfs/test_beta_pdf_gradient.C | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_pdfs/test_beta_pdf_gradient.C b/test/test_pdfs/test_beta_pdf_gradient.C index 4bc40f9c3..058e22d71 100644 --- a/test/test_pdfs/test_beta_pdf_gradient.C +++ b/test/test_pdfs/test_beta_pdf_gradient.C @@ -57,12 +57,15 @@ int main(int argc, char ** argv) 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); From 03cfacfb59ddc181d17faa78f2c00b41d36d79af Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Fri, 4 Mar 2016 13:58:02 -0600 Subject: [PATCH 21/47] Fix issue involving unified MLE on non-inter0Comm This patch adds some safeguards to avoid code that does inter0Comm computing from being executed by processes that don't belong to the inter0Comm. Specifically, computing the unified MLE (and positions) and computing correlation matrices are to only be done on the inter0Comm, since this is where all the samples live anyway. --- src/stats/src/MetropolisHastingsSG.C | 8 +++---- src/stats/src/StatisticalForwardProblem.C | 26 +++++++++++++---------- 2 files changed, 19 insertions(+), 15 deletions(-) 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 From aedb12b2bc023ede5b8b0e5ea9ba8b05a1177e2e Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Mon, 28 Mar 2016 16:14:20 -0500 Subject: [PATCH 22/47] Adding parallel test for MLE on intercomm0 --- .gitignore | 6 + configure.ac | 9 + test/Makefile.am | 15 ++ test/test_intercomm0/gravity_1proc.txt.in | 73 ++++++++ test/test_intercomm0/gravity_2proc.txt.in | 73 ++++++++ test/test_intercomm0/gravity_likelihood.C | 93 ++++++++++ test/test_intercomm0/gravity_likelihood.h | 48 +++++ test/test_intercomm0/gravity_main.C | 174 ++++++++++++++++++ test/test_intercomm0/gravity_qoi.C | 79 ++++++++ test/test_intercomm0/gravity_qoi.h | 60 ++++++ .../test_intercomm0_gravity_run.sh.in | 18 ++ 11 files changed, 648 insertions(+) create mode 100644 test/test_intercomm0/gravity_1proc.txt.in create mode 100644 test/test_intercomm0/gravity_2proc.txt.in create mode 100644 test/test_intercomm0/gravity_likelihood.C create mode 100644 test/test_intercomm0/gravity_likelihood.h create mode 100644 test/test_intercomm0/gravity_main.C create mode 100644 test/test_intercomm0/gravity_qoi.C create mode 100644 test/test_intercomm0/gravity_qoi.h create mode 100644 test/test_intercomm0/test_intercomm0_gravity_run.sh.in diff --git a/.gitignore b/.gitignore index 14f25904a..5512bb5df 100644 --- a/.gitignore +++ b/.gitignore @@ -218,3 +218,9 @@ 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 diff --git a/configure.ac b/configure.ac index 8e2d95948..1321f5872 100644 --- a/configure.ac +++ b/configure.ac @@ -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,13 @@ 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 +]) + dnl ---------------------------------------------- dnl Collect files for licence header stamping here dnl ---------------------------------------------- diff --git a/test/Makefile.am b/test/Makefile.am index 298347908..3e1000a66 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -58,6 +58,7 @@ 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 LDADD = $(top_builddir)/src/libqueso.la @@ -142,6 +143,13 @@ 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_main.C +test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_qoi.C +test_intercomm0_gravity_CPPFLAGS = -Itest_intercomm0 $(AM_CPPFLAGS) + + # Files to freedom stamp srcstamp = srcstamp += $(test_uqEnvironmentNonFatal_SOURCES) @@ -196,6 +204,7 @@ 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) TESTS = TESTS += $(top_builddir)/test/test_uqEnvironmentNonFatal @@ -251,6 +260,7 @@ 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 XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED @@ -285,6 +295,9 @@ 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 CLEANFILES = CLEANFILES += test_Environment/debug_output_sub0.txt @@ -302,6 +315,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/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 From 6905451a227c4db67969394a0d29077ba5cd0028 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Sat, 26 Mar 2016 17:10:47 -0400 Subject: [PATCH 23/47] Add queso-config --- src/apps/queso-config.in | 108 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 src/apps/queso-config.in diff --git a/src/apps/queso-config.in b/src/apps/queso-config.in new file mode 100644 index 000000000..e1ba45e51 --- /dev/null +++ b/src/apps/queso-config.in @@ -0,0 +1,108 @@ +#!/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 --fc" + echo " $0 --fflags" + 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" + ;; + + "--f77") + return_val="@F77@ $return_val" + ;; + + "--fc") + return_val="@FC@ $return_val" + ;; + + "--cppflags") + return_val="${CPPFLAGS} $return_val" + ;; + + "--cxxflags") + return_val="@CXXFLAGS@ $return_val" + ;; + + "--cflags") + return_val="${CFLAGS} $return_val" + ;; + + "--fflags") + return_val="@FFLAGS@ $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: From 38d33e8314f4484d6e0f3ee574e49a2639c0aa6c Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Sat, 26 Mar 2016 17:12:09 -0400 Subject: [PATCH 24/47] Add queso-config to build system --- Makefile.am | 4 ++++ configure.ac | 2 ++ 2 files changed, 6 insertions(+) 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 1321f5872..6a6e94587 100644 --- a/configure.ac +++ b/configure.ac @@ -279,6 +279,8 @@ AC_CONFIG_FILES([ chmod +x test/test_intercomm0/test_intercomm0_gravity_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 ---------------------------------------------- From 2112e42342166195d0c5b88eca1e1713aa7ed534 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Mon, 28 Mar 2016 21:28:38 -0400 Subject: [PATCH 25/47] Display CPPFLAGS with queso-config, configure Internally we use QUESO_CPPFLAGS because somehow CPPFLAGS doesn't get propagated through to Automake. This required a s/QUESO_CPPFLAGS/AM_CPPFLAGS in the examples/Makefile.am. Now, everyone should be respecting QUESO_CPPFLAGS in the main, build, examples, and tests. We also now see it in the configure summary. --- configure.ac | 4 +- examples/Makefile.am | 71 +++++++++++++------------- m4/config_summary.m4 | 1 + src/Makefile.am | 2 +- src/apps/queso-config.in | 2 +- test/Makefile.am | 2 +- test/gsl_tests/Makefile.am | 2 +- test/t01_valid_cycle/Makefile.am | 2 +- test/t02_sip_sfp/Makefile.am | 2 +- test/t03_sequence/Makefile.am | 2 +- test/t03_sequence/sequence/Makefile.am | 2 +- test/t04_bimodal/Makefile.am | 2 +- 12 files changed, 48 insertions(+), 46 deletions(-) diff --git a/configure.ac b/configure.ac index 6a6e94587..f69364d39 100644 --- a/configure.ac +++ b/configure.ac @@ -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) 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/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 index e1ba45e51..5efe4d1d9 100644 --- a/src/apps/queso-config.in +++ b/src/apps/queso-config.in @@ -59,7 +59,7 @@ while [ "x$1" != "x" ]; do ;; "--cppflags") - return_val="${CPPFLAGS} $return_val" + return_val="@QUESO_CPPFLAGS@ $return_val" ;; "--cxxflags") diff --git a/test/Makefile.am b/test/Makefile.am index 3e1000a66..2aebde59a 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -62,7 +62,7 @@ check_PROGRAMS += test_intercomm0_gravity 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) 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/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) From ff4c890780da08121ac9ca92b28e2f142141f231 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Mon, 28 Mar 2016 21:58:58 -0400 Subject: [PATCH 26/47] Remove Fortran bits from queso-config --- src/apps/queso-config.in | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/apps/queso-config.in b/src/apps/queso-config.in index 5efe4d1d9..a1768bf21 100644 --- a/src/apps/queso-config.in +++ b/src/apps/queso-config.in @@ -19,8 +19,6 @@ usage () echo "usage: $0 --cppflags --cxxflags --include --libs" echo " $0 --cxx" echo " $0 --cc" - echo " $0 --fc" - echo " $0 --fflags" echo " $0 --version" echo " $0 --host" echo " $0 --ldflags" @@ -50,14 +48,6 @@ while [ "x$1" != "x" ]; do return_val="@CC@ $return_val" ;; - "--f77") - return_val="@F77@ $return_val" - ;; - - "--fc") - return_val="@FC@ $return_val" - ;; - "--cppflags") return_val="@QUESO_CPPFLAGS@ $return_val" ;; @@ -70,10 +60,6 @@ while [ "x$1" != "x" ]; do return_val="${CFLAGS} $return_val" ;; - "--fflags") - return_val="@FFLAGS@ $return_val" - ;; - "--include") return_val="-I${includedir} $return_val" ;; From 0a2e71df23f18096e29237bb9773b6ad1779c996 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Wed, 16 Mar 2016 19:48:43 -0500 Subject: [PATCH 27/47] Add HDF5 file creation within openOutputFile This also had to rejig some checks for non-NULLness on the ofstream object. They were checked even when the file type was set to h5. --- src/core/src/Environment.C | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/core/src/Environment.C b/src/core/src/Environment.C index f325a81d8..9e7ebad2b 100644 --- a/src/core/src/Environment.C +++ b/src/core/src/Environment.C @@ -600,9 +600,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"); } @@ -653,7 +666,15 @@ BaseEnvironment::openOutputFile( << "'" << std::endl; } - queso_require_msg((filePtrSet.ofsVar && filePtrSet.ofsVar->is_open()), "failed to open output file"); + + // 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"); + } + } else { returnValue = false; From 775bef1ef23816f5d8fcc9fa82a42b86c9b4e7ec Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Wed, 16 Mar 2016 19:54:29 -0500 Subject: [PATCH 28/47] Initialise hid_t variable when HDF5 is present --- src/core/src/Environment.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/core/src/Environment.C b/src/core/src/Environment.C index 9e7ebad2b..77b1f0027 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 { } From 80a83af4fe12395b9829f1ed54ffcbc8ced0b01d Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 19 Mar 2016 00:04:11 -0400 Subject: [PATCH 29/47] Inconsequential boredom --- src/core/src/GslVector.C | 2 -- 1 file changed, 2 deletions(-) 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& From e4c9a2ffbc822cf707911478a4b2f41228b25f2c Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 19 Mar 2016 00:04:51 -0400 Subject: [PATCH 30/47] Adding HDF5 output support for subchains This relies on an awful chainSize * numParams *COPY*. This is completely and utterly the wrong thing to do, but I wanted to get something working. Most of the logic here was actually copied from the unified HDF5 write, which doesn't seem to be working at the time of this commit. --- src/basic/src/SequenceOfVectors.C | 84 +++++++++++++++++++++++++++++-- 1 file changed, 80 insertions(+), 4 deletions(-) diff --git a/src/basic/src/SequenceOfVectors.C b/src/basic/src/SequenceOfVectors.C index a97e2c243..78e8755f2 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, + "samples", + 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"); } From 1d53d244a133571f1aa7c329a70d95e692b46d29 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 19 Mar 2016 00:18:29 -0400 Subject: [PATCH 31/47] Remove unneeded logging We throw an exception if we can't open the file anyway. --- src/core/src/Environment.C | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/src/core/src/Environment.C b/src/core/src/Environment.C index 77b1f0027..9c42ee125 100644 --- a/src/core/src/Environment.C +++ b/src/core/src/Environment.C @@ -654,22 +654,6 @@ 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) || @@ -678,7 +662,6 @@ BaseEnvironment::openOutputFile( (filePtrSet.ofsVar && filePtrSet.ofsVar->is_open()), "failed to open output file"); } - } else { returnValue = false; From 6ba03f51138ff411358b2064de0975da7aa68c28 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 26 Mar 2016 15:46:43 -0500 Subject: [PATCH 32/47] Only use the ofsVar ostream if we're writing ascii --- src/basic/src/ScalarSequence.C | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/basic/src/ScalarSequence.C b/src/basic/src/ScalarSequence.C index 532ee6ef2..80e565012 100644 --- a/src/basic/src/ScalarSequence.C +++ b/src/basic/src/ScalarSequence.C @@ -2593,10 +2593,18 @@ 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); + } + else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { + // Do nothing + } + m_env.closeFile(filePtrSet,fileType); } From 93a29129717b89263be6b041acaffd92e14eb6e1 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 26 Mar 2016 15:47:28 -0500 Subject: [PATCH 33/47] Remove unnecessary return statement --- src/basic/src/ScalarSequence.C | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/basic/src/ScalarSequence.C b/src/basic/src/ScalarSequence.C index 80e565012..1183291bc 100644 --- a/src/basic/src/ScalarSequence.C +++ b/src/basic/src/ScalarSequence.C @@ -2607,8 +2607,6 @@ ScalarSequence::subWriteContents( m_env.closeFile(filePtrSet,fileType); } - - return; } // -------------------------------------------------- template From 2c2ed6e9b6d775d06f718ec96b31b15a4649e47d Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 26 Mar 2016 16:47:57 -0500 Subject: [PATCH 34/47] Implement HDF5 output for scalar sequences --- src/basic/src/ScalarSequence.C | 43 +++++++++++++++++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/src/basic/src/ScalarSequence.C b/src/basic/src/ScalarSequence.C index 1183291bc..3bdef4324 100644 --- a/src/basic/src/ScalarSequence.C +++ b/src/basic/src/ScalarSequence.C @@ -2601,9 +2601,50 @@ ScalarSequence::subWriteContents( *filePtrSet.ofsVar, fileType); } +#ifdef QUESO_HAS_HDF5 else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) { - // Do nothing + + // 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, + "samples", + 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); } From 4084cb4669d2534ea45b507db05affee5d2e9510 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 26 Mar 2016 16:48:45 -0500 Subject: [PATCH 35/47] Make scalar output rank 1, not rank 2 --- src/basic/src/ScalarSequence.C | 28 ++++------------------------ 1 file changed, 4 insertions(+), 24 deletions(-) diff --git a/src/basic/src/ScalarSequence.C b/src/basic/src/ScalarSequence.C index 3bdef4324..e18eb3f94 100644 --- a/src/basic/src/ScalarSequence.C +++ b/src/basic/src/ScalarSequence.C @@ -2770,10 +2770,9 @@ 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", @@ -2789,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, @@ -2811,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; @@ -2838,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"); From 980d28e54ef173cba7b2ce16ebab5e721f71389a Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 26 Mar 2016 16:49:33 -0500 Subject: [PATCH 36/47] Write chainSize x numParams slab not the transpose --- src/basic/src/SequenceOfVectors.C | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/src/basic/src/SequenceOfVectors.C b/src/basic/src/SequenceOfVectors.C index 78e8755f2..567e62d85 100644 --- a/src/basic/src/SequenceOfVectors.C +++ b/src/basic/src/SequenceOfVectors.C @@ -1621,8 +1621,8 @@ 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, @@ -1639,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, @@ -1660,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; @@ -1685,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"); From 411d91b493d4eb5614db7962bd8ecf1a891c7fe4 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Sat, 26 Mar 2016 17:15:41 -0500 Subject: [PATCH 37/47] Call all written hdf5 datasets 'data' --- src/basic/src/ScalarSequence.C | 6 +++--- src/basic/src/SequenceOfVectors.C | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/basic/src/ScalarSequence.C b/src/basic/src/ScalarSequence.C index e18eb3f94..0a5b0c2cb 100644 --- a/src/basic/src/ScalarSequence.C +++ b/src/basic/src/ScalarSequence.C @@ -2614,7 +2614,7 @@ ScalarSequence::subWriteContents( // Create dataset hid_t dataset_id = H5Dcreate(filePtrSet.h5Var, - "samples", + "data", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, @@ -2775,7 +2775,7 @@ ScalarSequence::unifiedWriteContents( 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 @@ -3070,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 567e62d85..415a89c0b 100644 --- a/src/basic/src/SequenceOfVectors.C +++ b/src/basic/src/SequenceOfVectors.C @@ -1353,7 +1353,7 @@ SequenceOfVectors::subWriteContents( // Create dataset hid_t dataset_id = H5Dcreate(filePtrSet.h5Var, - "samples", + "data", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, @@ -1626,7 +1626,7 @@ SequenceOfVectors::unifiedWriteContents( 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 @@ -1924,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); From d64dd8616aee8d07470de5a7f335fbd7903d9664 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 29 Mar 2016 09:19:45 -0500 Subject: [PATCH 38/47] Add test for successful HDF5 output write --- .gitignore | 1 + test/Makefile.am | 5 ++ test/test_SequenceOfVectors/test_HDF5Write.C | 66 ++++++++++++++++++++ 3 files changed, 72 insertions(+) create mode 100644 test/test_SequenceOfVectors/test_HDF5Write.C diff --git a/.gitignore b/.gitignore index 5512bb5df..fd39c42a2 100644 --- a/.gitignore +++ b/.gitignore @@ -224,3 +224,4 @@ 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 diff --git a/test/Makefile.am b/test/Makefile.am index 3e1000a66..c146da048 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -59,6 +59,7 @@ 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 @@ -149,6 +150,7 @@ test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_main.C test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_qoi.C test_intercomm0_gravity_CPPFLAGS = -Itest_intercomm0 $(AM_CPPFLAGS) +test_seq_of_vec_hdf5_write_SOURCES = test_SequenceOfVectors/test_HDF5Write.C # Files to freedom stamp srcstamp = @@ -205,6 +207,7 @@ 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 @@ -261,6 +264,7 @@ 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_seq_of_vec_hdf5_write XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED @@ -304,6 +308,7 @@ CLEANFILES += test_Environment/debug_output_sub0.txt CLEANFILES += gslvector_out_sub0.m CLEANFILES += test_write_InterpolationSurrogateBuilder_1.dat CLEANFILES += test_write_InterpolationSurrogateBuilder_2.dat +CLEANFILES += output_test_hdf5_sub0.h5 clean-local: rm -rf $(top_builddir)/test/chain0 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 +} From 69d8e9da45ee84fcea580491dca437d0162275af Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 29 Mar 2016 13:15:57 -0500 Subject: [PATCH 39/47] Update changelog --- CHANGES | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGES b/CHANGES index be1920e46..ea6d7563b 100644 --- a/CHANGES +++ b/CHANGES @@ -3,6 +3,9 @@ QUESO: Quantification of Uncertainty for Estimation, Simulation, and Optimization. ----------------------------------------------------- +Version 0.55.0 + * Add support for HDF5 output + Version 0.54.0 * Fix memory leak in test_GslVector * Make MPI optional From bfb955d1116a9cdfb3146357b29e731fece737a7 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 29 Mar 2016 16:36:09 -0500 Subject: [PATCH 40/47] Bump version number --- configure.ac | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configure.ac b/configure.ac index f69364d39..749582e58 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" From 561a6acdd79128bd5c1c5a9d86ce463b8c0b4bac Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 29 Mar 2016 16:41:42 -0500 Subject: [PATCH 41/47] Update changelog --- CHANGES | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGES b/CHANGES index ea6d7563b..ac574160f 100644 --- a/CHANGES +++ b/CHANGES @@ -4,6 +4,9 @@ Simulation, and Optimization. ----------------------------------------------------- Version 0.55.0 + * 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 From 6e42a85b7a2ff29d1c4caa5f2affb7a72bd5476a Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Wed, 30 Mar 2016 10:13:41 -0500 Subject: [PATCH 42/47] Make HDF5 write test clean up output file --- .gitignore | 1 + configure.ac | 7 +++++++ test/Makefile.am | 4 ++-- .../test_seq_of_vec_hdf5_write_run.sh.in | 6 ++++++ 4 files changed, 16 insertions(+), 2 deletions(-) create mode 100644 test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh.in diff --git a/.gitignore b/.gitignore index fd39c42a2..5f67726c8 100644 --- a/.gitignore +++ b/.gitignore @@ -225,3 +225,4 @@ 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/configure.ac b/configure.ac index f69364d39..c32b53c0a 100644 --- a/configure.ac +++ b/configure.ac @@ -279,6 +279,13 @@ AC_CONFIG_FILES([ 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 ---------------------------------------------- diff --git a/test/Makefile.am b/test/Makefile.am index 8d80bbb45..1c97f7039 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -264,7 +264,7 @@ 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_seq_of_vec_hdf5_write +TESTS += $(top_builddir)/test/test_SequenceOfVectors/test_seq_of_vec_hdf5_write_run.sh XFAIL_TESTS = $(top_builddir)/test/test_SequenceOfVectorsErase if ! MPI_ENABLED @@ -302,13 +302,13 @@ 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 CLEANFILES += gslvector_out_sub0.m CLEANFILES += test_write_InterpolationSurrogateBuilder_1.dat CLEANFILES += test_write_InterpolationSurrogateBuilder_2.dat -CLEANFILES += output_test_hdf5_sub0.h5 clean-local: rm -rf $(top_builddir)/test/chain0 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 From 6524230b4df09859628026d0c2a10fd36f6440c1 Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Wed, 30 Mar 2016 09:26:29 -0500 Subject: [PATCH 43/47] Fix "make check" for out-of-source builds --- test/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Makefile.am b/test/Makefile.am index 8d80bbb45..5937c8814 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -148,7 +148,7 @@ test_intercomm0_gravity_SOURCES = test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_likelihood.C test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_main.C test_intercomm0_gravity_SOURCES += test_intercomm0/gravity_qoi.C -test_intercomm0_gravity_CPPFLAGS = -Itest_intercomm0 $(AM_CPPFLAGS) +test_intercomm0_gravity_CPPFLAGS = -I$(srcdir)/test_intercomm0 $(AM_CPPFLAGS) test_seq_of_vec_hdf5_write_SOURCES = test_SequenceOfVectors/test_HDF5Write.C From 19d42475a3a36ef04c16893c1185137b195d7313 Mon Sep 17 00:00:00 2001 From: Nicholas Malaya Date: Thu, 31 Mar 2016 17:46:07 -0500 Subject: [PATCH 44/47] queso: updating changelog to include date of release * Also cleaned up log so that each release has consistently form * Removed redundant release logging for 0.47.0 --- CHANGES | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/CHANGES b/CHANGES index ac574160f..2384e4108 100644 --- a/CHANGES +++ b/CHANGES @@ -3,13 +3,10 @@ QUESO: Quantification of Uncertainty for Estimation, Simulation, and Optimization. ----------------------------------------------------- -Version 0.55.0 - * 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 +Version 0.55.0 (in progress) * Add support for HDF5 output -Version 0.54.0 +Version 0.54.0 (Oct 8, 2015) * Fix memory leak in test_GslVector * Make MPI optional * Optimise GSLMatrix::operator() and GslVector::operator[] @@ -19,7 +16,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 @@ -27,42 +24,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 From 7081a36148135dfea659a5d7712cfe29bacd8500 Mon Sep 17 00:00:00 2001 From: Nicholas Malaya Date: Thu, 31 Mar 2016 17:48:43 -0500 Subject: [PATCH 45/47] missed that the changelog had been updated... --- CHANGES | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGES b/CHANGES index 2384e4108..bac940af3 100644 --- a/CHANGES +++ b/CHANGES @@ -4,6 +4,9 @@ Simulation, and Optimization. ----------------------------------------------------- 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) From c96b02a76caabc2c43192d8eb8c06a6c1b88b5f0 Mon Sep 17 00:00:00 2001 From: "Roy H. Stogner" Date: Sat, 2 Apr 2016 00:09:08 -0500 Subject: [PATCH 46/47] Add missing files to fix distcheck --- test/Makefile.am | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/Makefile.am b/test/Makefile.am index 6c1db6bfb..de23f7101 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -146,8 +146,10 @@ 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 @@ -291,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 From 8b8d02fa78b7e0c0448aeac2e36d62c17b6e955b Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Thu, 16 Jun 2016 17:24:34 -0500 Subject: [PATCH 47/47] Fix a race condition in gsl test output directory The output isn't even needed, so this patch just turns it off. --- test/gsl_tests/input.in | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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