-
Notifications
You must be signed in to change notification settings - Fork 40
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
GPMSA Example: Brian Williams' linear verification #580
Conversation
Add initial version of Bayes linear verification problem due to Brian Williams [1] Adams, B.M., Coleman, Kayla, Hooper, Russell W., Khuwaileh, Bassam A., Lewis, Allison, Smith, Ralph C., Swiler, Laura P., Turinsky, Paul J., Williams, Brian J., "User Guidelines and Best Practices for CASL VUQ Analysis Using Dakota, Sandia Technical Report 2016-11614 (also CASL-U-2016-1233-000), September 2016.
Codecov Report
@@ Coverage Diff @@
## dev #580 +/- ##
==========================================
+ Coverage 71.15% 71.33% +0.17%
==========================================
Files 167 167
Lines 14457 14453 -4
==========================================
+ Hits 10287 10310 +23
+ Misses 4170 4143 -27
Continue to review full report at Codecov.
|
Spoke with @roystgnr who suggested if possible we make some or all of this a test, with XFAIL status if needed. While the ultimate test @brianw525 has in mind likely won't be feasible for this (requires 100k or 250k samples) and is why I targeted examples/, we can make a scaled down version to use as at least a regression test. |
Do you want me to merge this example or wait for a test? |
Yes, but I defer to you all as to what's appropriate in the project. One of my goals prior to our next in-person meeting is to implement a few of Brian Williams' verification cases so we can actually iterate on working code in the meeting and run cases. Since the cases have to run many samples, perhaps for now they belong as examples and not tests and that's just okay? Also, since I'll be adding a few more cases similar to this, I am keenly interested in feedback on my implementation before I copy/paste and make other cases, but that's not required. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You can add your example to this file if you want it to be easily executed from our run_all.sh
script.
Caveat other remarks, this looks fine to me. Thanks!
############################################### | ||
#ip_mh_help = anything | ||
ip_mh_dataOutputFileName = outputData/sipOutput | ||
ip_mh_dataOutputAllowedSet = 0 1 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Vectors need single quotes around them ('0 1'
) becausedev
uses GetPot
to parse input files, and not boost
. You can still use boost
if you want, but the default is GetPot
.
ip_mh_rawChain_displayPeriod = 2000 | ||
ip_mh_rawChain_measureRunTimes = 1 | ||
ip_mh_rawChain_dataOutputFileName = outputData/ip_raw_chain | ||
ip_mh_rawChain_dataOutputAllowedSet = 0 1 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Needs quotes.
ip_mh_tk_useLocalHessian = 0 | ||
ip_mh_tk_useNewtonComponent = 0 | ||
ip_mh_dr_maxNumExtraStages = 0 | ||
ip_mh_dr_listOfScalesForExtraStages = 5. 10. 20. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Needs quotes.
@@ -0,0 +1,345 @@ | |||
// Bayes linear verification problem due to Brian Williams |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can you put the QUESO licence before this? You can copy it from one of the other examples.
Thanks for the helpful feedback. Brian Williams updated this example and gave me several others, so I'm thinking I'll also reorganize the directories and do some cleanup, maybe even build some utilities that can be shared among the cases, before we merge... |
No worries. Let me know when you're ready to pull the trigger. |
Address feedback from @dmcdougall in libqueso#580 * Add license * Quote vectors in text input file
Refactor linear verification data read to work for both scalar and multivariate cases (number of responses > 1). Rename input and data files in preparation for adding multivariate case.
Minor rework of main to drive scalar vs. multivariate cases.
Add multivariate five-response linear verification case (example toggles based on input filename for now), including observation error covariance matrix. Case is almost same as scalar, but just copied the code for now to make debugging more straightforward. Also add comments with @brianw525 recommended changes to priors and initial point (tried changing them, but then scalar case doesn't run)
@@ -11,7 +11,7 @@ | |||
// 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 | |||
// but WITHOUT ANY WARRANT Y; without even the implied warranty of |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
please fix thx
outputVecs, experimentScenarios, experimentVecs); | ||
|
||
// Experiment observation error | ||
for (unsigned int i = 0; i < numExperiments; i++) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
numExperiments
should be numExperiments * experimentSize
in the general case.
|
||
// Experiment observation error | ||
for (unsigned int i = 0; i < numExperiments; i++) | ||
(*experimentMat)(i, i) = std::pow(0.05, 2.0); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If we're normalising the experiments by the mean and variance of the simulations, we also need to normalise this matrix by the simulation variance.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since this is a matrix, we'd pre- and post-mulitply by 1/sqrt(diag(S))
where division is done component-wise.
// paramInitials[6] = std::exp(-0.025); // emul corr strength (scenario/beta?) | ||
// paramInitials[7] = std::exp(-0.025); // emul corr strength (scenario/beta?) | ||
// paramInitials[8] = std::exp(-0.025); // emul corr strength (scenario/beta?) | ||
// paramInitials[9] = std::exp(-0.025); // emul corr strength (scenario/beta?) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The commented out code is correct. Please replace elements 4 through 9 above.
// paramInitials[9] = std::exp(-0.025); // emul corr strength (scenario/beta?) | ||
|
||
paramInitials[10] = 1000.0; // discrepancy precision | ||
// paramInitials[10] = 20.0; // discrepancy precision |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ditto. 20
should be the default and for the verification.
paramInitials[13] = 0.1; // discrepancy corr strength x3 | ||
// paramInitials[11] = std::exp(-0.025); // discrepancy corr strength x1 | ||
// paramInitials[12] = std::exp(-0.025); // discrepancy corr strength x2 | ||
// paramInitials[13] = std::exp(-0.025); // discrepancy corr strength x3 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ditto.
// paramInitials[13] = std::exp(-0.025); // discrepancy corr strength x3 | ||
|
||
paramInitials[14] = 10000.0; // emul data precision | ||
// paramInitials[14] = 1000.0; // emul data precision |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ditto.
paramInitials[5] = 1.0; // weights0 precision | ||
paramInitials[6] = 1.0; // weights0 precision | ||
paramInitials[7] = 1.0; // weights0 precision | ||
paramInitials[8] = 1.0; // weights0 precision |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
By this point we already know how many basis vectors there are. QUESO just needs to expose enough to let the user query it.
* Use @brianw525 recommended prior and initial point settings. * Scale experiment observation error covariance by simulation variance.
examples/Makefile.am
Outdated
gpmsa_linverif_CPPFLAGS = -I$(top_srcdir)/examples/gp/linear_verif $(AM_CPPFLAGS) | ||
|
||
dist_gpmsa_linverif_DATA = | ||
dist_gpmsa_linverif_DATA += ${gpmsa_scalar_SOURCES} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This installs gpmsa_scalar.C, not gpmsa_linear_verif.C
Minor edits per reviewer feedback
@roystgnr @dmcdougall @brianw525 This starts to address #579 by implementing Brian's test problem in a QUESO example. It runs and doesn't crash (issues a bunch of warnings from one particular spot though.)
In terms of your development workflow, I wasn't sure if you want this to be managed on a branch until the test is verified vs. merged to devel. Feel free to guide me as to how to work.
Add initial version of Bayes linear verification problem due to Brian
Williams
[1] Adams, B.M., Coleman, Kayla, Hooper, Russell W., Khuwaileh, Bassam
A., Lewis, Allison, Smith, Ralph C., Swiler, Laura P., Turinsky, Paul
J., Williams, Brian J., "User Guidelines and Best Practices for CASL
VUQ Analysis Using Dakota, Sandia Technical Report 2016-11614 (also
CASL-U-2016-1233-000), September 2016.