From cda033f8b3198f2a3c56e17954f3bfa91fc1f39d Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Wed, 18 Oct 2017 09:13:33 -0500 Subject: [PATCH 1/6] Fix emulator correlation strength indexing bug We were assuming only one basis vector when indexing the emulator correlation strength. --- src/gp/src/GPMSA.C | 72 ++++++++++++++++++++++++---------------------- 1 file changed, 37 insertions(+), 35 deletions(-) diff --git a/src/gp/src/GPMSA.C b/src/gp/src/GPMSA.C index 8925f2dba..106d858bc 100644 --- a/src/gp/src/GPMSA.C +++ b/src/gp/src/GPMSA.C @@ -245,48 +245,50 @@ GPMSAEmulator::lnValue(const V & domainVector, (this->m_simulationParameters)[j-this->m_numExperiments]; } - // Emulator component // = first term in (1) - prodScenario = 1.0; unsigned int emulatorCorrStrStart = dimParameter + (this->num_svd_terms < num_nonzero_eigenvalues) + num_svd_terms; - for (unsigned int k = 0; k < dimScenario; k++) { - const double & emulator_corr_strength = - domainVector[emulatorCorrStrStart+k]; - double scenario_param1 = - m_opts.normalized_scenario_parameter(k, (*scenario1)[k]); - double scenario_param2 = - m_opts.normalized_scenario_parameter(k, (*scenario2)[k]); - prodScenario *= std::pow(emulator_corr_strength, - 4.0 * (scenario_param1 - scenario_param2) * - (scenario_param1 - scenario_param2)); - } - - queso_assert (!queso_isnan(prodScenario)); - - // = second term in (1) - prodParameter = 1.0; - for (unsigned int k = 0; k < dimParameter; k++) { - queso_assert (!queso_isnan(domainVector[emulatorCorrStrStart+dimScenario+k])); - queso_assert (!queso_isnan((*parameter1)[k])); - queso_assert (!queso_isnan((*parameter2)[k])); - const double & emulator_corr_strength = - domainVector[emulatorCorrStrStart+dimScenario+k]; - double uncertain_param1 = - m_opts.normalized_uncertain_parameter(k, (*parameter1)[k]); - double uncertain_param2 = - m_opts.normalized_uncertain_parameter(k, (*parameter2)[k]); - prodParameter *= std::pow( - emulator_corr_strength, - 4.0 * (uncertain_param1 - uncertain_param2) * - (uncertain_param1 - uncertain_param2)); - } - - queso_assert (!queso_isnan(prodParameter)); // Sigma_eta in scalar case, // [Sigma_u, Sigma_uw; Sigma_uw^T, Sigma_w] in vector case for (unsigned int basis = 0; basis != num_svd_terms; ++basis) { + + // Emulator component // = first term in (1) + prodScenario = 1.0; + for (unsigned int k = 0; k < dimScenario; k++) { + const double & emulator_corr_strength = + domainVector[emulatorCorrStrStart+(dimScenario+dimParameter)*basis+k]; + double scenario_param1 = + m_opts.normalized_scenario_parameter(k, (*scenario1)[k]); + double scenario_param2 = + m_opts.normalized_scenario_parameter(k, (*scenario2)[k]); + prodScenario *= std::pow(emulator_corr_strength, + 4.0 * (scenario_param1 - scenario_param2) * + (scenario_param1 - scenario_param2)); + } + + queso_assert (!queso_isnan(prodScenario)); + + // = second term in (1) + prodParameter = 1.0; + for (unsigned int k = 0; k < dimParameter; k++) { + queso_assert (!queso_isnan(domainVector[emulatorCorrStrStart+(dimScenario+dimParameter)*basis+dimScenario+k])); + queso_assert (!queso_isnan((*parameter1)[k])); + queso_assert (!queso_isnan((*parameter2)[k])); + const double & emulator_corr_strength = + domainVector[emulatorCorrStrStart+(dimScenario+dimParameter)*basis+dimScenario+k]; + double uncertain_param1 = + m_opts.normalized_uncertain_parameter(k, (*parameter1)[k]); + double uncertain_param2 = + m_opts.normalized_uncertain_parameter(k, (*parameter2)[k]); + prodParameter *= std::pow( + emulator_corr_strength, + 4.0 * (uncertain_param1 - uncertain_param2) * + (uncertain_param1 - uncertain_param2)); + } + + queso_assert (!queso_isnan(prodParameter)); + // coefficient in (1) // The relevant precision for Sigma_eta is lambda_eta; for // Sigma_uw etc. it's lambda_wi and we skip lambda_eta From 68aac288402a1f51be902679a8ab6698eb0cd843 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Wed, 18 Oct 2017 09:38:05 -0500 Subject: [PATCH 2/6] Fix discrepancy correlation param indexing This relied on the correct calculation of how many emulator correlation strength parameters there were; we had assumed there was one svd term. --- src/gp/src/GPMSA.C | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/gp/src/GPMSA.C b/src/gp/src/GPMSA.C index 106d858bc..5b94f9419 100644 --- a/src/gp/src/GPMSA.C +++ b/src/gp/src/GPMSA.C @@ -317,7 +317,7 @@ GPMSAEmulator::lnValue(const V & domainVector, typename SharedPtr::Type cross_scenario1 = (this->m_experimentScenarios)[i]; typename SharedPtr::Type cross_scenario2 = (this->m_experimentScenarios)[j]; unsigned int discrepancyCorrStrStart = - dimParameter + num_svd_terms + dimParameter + dimScenario + num_discrepancy_groups + + dimParameter + num_svd_terms + (dimParameter + dimScenario) * num_svd_terms + num_discrepancy_groups + (this->num_svd_terms::lnValue(const V & domainVector, unsigned int discrepancyPrecisionStart = dimParameter + (num_svd_terms Date: Wed, 18 Oct 2017 09:40:06 -0500 Subject: [PATCH 3/6] Fix indexing error in dimSum variable Again, this relied on the calculation of the number of emulator correlation strength parameters. We had assumed there was one svd basis vector which is an incorrect assumption. --- src/gp/src/GPMSA.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/gp/src/GPMSA.C b/src/gp/src/GPMSA.C index 5b94f9419..01f8c667e 100644 --- a/src/gp/src/GPMSA.C +++ b/src/gp/src/GPMSA.C @@ -176,8 +176,7 @@ GPMSAEmulator::lnValue(const V & domainVector, m_opts.m_calibrateObservationalPrecision + num_svd_terms + dimParameter + - dimParameter + - dimScenario + + num_svd_terms * (dimParameter + dimScenario) + num_discrepancy_groups + (num_discrepancy_groups * dimScenario); // yum From 5eff358220365af46da2a021e68a4bc68162f1bb Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 31 Oct 2017 17:47:27 -0500 Subject: [PATCH 4/6] Fix emulator data precision indexing There are num_svd_terms of these things, not 1. --- src/gp/inc/GPMSA.h | 1 + src/gp/src/GPMSA.C | 30 ++++++++++++++++++++---------- 2 files changed, 21 insertions(+), 10 deletions(-) diff --git a/src/gp/inc/GPMSA.h b/src/gp/inc/GPMSA.h index 0ceaad179..be1541572 100644 --- a/src/gp/inc/GPMSA.h +++ b/src/gp/inc/GPMSA.h @@ -489,6 +489,7 @@ class GPMSAFactory typename ScopedPtr >::Type discrepancyCorrelationDomain; // Emulator data precision + typename ScopedPtr >::Type emulatorDataPrecisionSpace; typename ScopedPtr::Type emulatorDataPrecisionMin; typename ScopedPtr::Type emulatorDataPrecisionMax; typename ScopedPtr >::Type emulatorDataPrecisionDomain; diff --git a/src/gp/src/GPMSA.C b/src/gp/src/GPMSA.C index 01f8c667e..f035efb60 100644 --- a/src/gp/src/GPMSA.C +++ b/src/gp/src/GPMSA.C @@ -171,7 +171,7 @@ GPMSAEmulator::lnValue(const V & domainVector, unsigned int dimParameter = (this->m_parameterSpace).dimLocal(); // Length of prior+hyperprior inputs - unsigned int dimSum = 1 + + unsigned int dimSum = num_svd_terms + (this->num_svd_terms < num_nonzero_eigenvalues) + m_opts.m_calibrateObservationalPrecision + num_svd_terms + @@ -1686,14 +1686,21 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) *(this->m_discrepancyCorrelationStrengthBetaVec))); // Emulator data precision + this->emulatorDataPrecisionSpace.reset + (new VectorSpace + (this->m_env, + "", + num_svd_terms, + NULL)); + this->emulatorDataPrecisionMin.reset - (new V(this->oneDSpace->zeroVector())); + (new V(this->emulatorDataPrecisionSpace->zeroVector())); this->emulatorDataPrecisionMax.reset - (new V(this->oneDSpace->zeroVector())); + (new V(this->emulatorDataPrecisionSpace->zeroVector())); this->m_emulatorDataPrecisionShapeVec.reset - (new V(this->oneDSpace->zeroVector())); + (new V(this->emulatorDataPrecisionSpace->zeroVector())); this->m_emulatorDataPrecisionScaleVec.reset - (new V(this->oneDSpace->zeroVector())); + (new V(this->emulatorDataPrecisionSpace->zeroVector())); this->emulatorDataPrecisionMin->cwSet(60.0); this->emulatorDataPrecisionMax->cwSet(1e5); this->m_emulatorDataPrecisionShapeVec->cwSet(emulatorDataPrecisionShape); @@ -1702,7 +1709,7 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) this->emulatorDataPrecisionDomain.reset (new BoxSubset ("", - *(this->oneDSpace), + *(this->emulatorDataPrecisionSpace), *(this->emulatorDataPrecisionMin), *(this->emulatorDataPrecisionMax))); @@ -1715,7 +1722,7 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) // Now form full prior const unsigned int dimHyper = - 1 + + num_svd_terms + (this->num_svd_terms < num_nonzero_eigenvalues) + this->m_opts->m_calibrateObservationalPrecision + num_svd_terms + @@ -1766,9 +1773,12 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) } const int emulator_data_precision_index = - dimHyper - 1 - this->m_opts->m_calibrateObservationalPrecision; - (*(this->hyperparamMins))[emulator_data_precision_index] = 60.0; // Min emulator data precision - (*(this->hyperparamMaxs))[emulator_data_precision_index] = 1e5; // Max emulator data precision + dimHyper - num_svd_terms - this->m_opts->m_calibrateObservationalPrecision; + for (unsigned int emulator_basis = 0; emulator_basis < num_svd_terms; + emulator_basis++) { + (*(this->hyperparamMins))[emulator_data_precision_index+emulator_basis] = 60.0; // Min emulator data precision + (*(this->hyperparamMaxs))[emulator_data_precision_index+emulator_basis] = 1e5; // Max emulator data precision + } if (this->m_opts->m_calibrateObservationalPrecision) { (*(this->hyperparamMins))[dimHyper-1] = 0.3; // Min observation error precision From 83f9393a038f88f6fae79467af8a9b074bc7ede1 Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 31 Oct 2017 17:48:02 -0500 Subject: [PATCH 5/6] Throw in an assert for good measure --- src/gp/src/GPMSA.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/gp/src/GPMSA.C b/src/gp/src/GPMSA.C index f035efb60..e4c7f9260 100644 --- a/src/gp/src/GPMSA.C +++ b/src/gp/src/GPMSA.C @@ -180,6 +180,8 @@ GPMSAEmulator::lnValue(const V & domainVector, num_discrepancy_groups + (num_discrepancy_groups * dimScenario); // yum + queso_assert_equal_to(dimSum, domainVector.sizeLocal()); + // Offset for Sigma_eta equivalent in vector case const unsigned int offset1 = (numSimulationOutputs == 1) ? 0 : m_numExperiments * num_discrepancy_bases; From 40d390ffe11d535c2d0172a759fc31ccdc4ce4af Mon Sep 17 00:00:00 2001 From: Damon McDougall Date: Tue, 31 Oct 2017 17:48:30 -0500 Subject: [PATCH 6/6] Fix emulator correlation strength indexing We were underestimating how many of these variables there were. There aren't always dimScenario+dimParameter. There's dimScenario+dimParameter *for each* emulator precision variable. --- src/gp/src/GPMSA.C | 8 +- .../test_gpmsa_autoscaled_samples.m | 2002 ++++++++--------- test/test_gpmsa/test_gpmsa_vector.C | 60 +- test/test_gpmsa/test_gpmsa_vector_samples.m | 2002 ++++++++--------- 4 files changed, 2041 insertions(+), 2031 deletions(-) diff --git a/src/gp/src/GPMSA.C b/src/gp/src/GPMSA.C index e4c7f9260..b076aa00f 100644 --- a/src/gp/src/GPMSA.C +++ b/src/gp/src/GPMSA.C @@ -1557,7 +1557,7 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) (new VectorSpace (this->m_env, "", - dimScenario + dimParameter, + (dimScenario + dimParameter) * num_svd_terms, NULL)); this->emulatorCorrelationMin.reset @@ -1728,8 +1728,7 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) (this->num_svd_terms < num_nonzero_eigenvalues) + this->m_opts->m_calibrateObservationalPrecision + num_svd_terms + - dimParameter + - dimScenario + + (dimParameter + dimScenario) * num_svd_terms + num_discrepancy_groups + (num_discrepancy_groups*dimScenario); // yum @@ -1764,8 +1763,7 @@ GPMSAFactory::setUpHyperpriors(const M & Wy) // Starting index of the discrepancy precisions in the hyperparameter vector unsigned int discrepancyPrecisionIdx = (num_svd_terms