1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24
25 #include <queso/Defines.h>
26
27 #ifndef QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
28 #include <boost/program_options.hpp>
29 #else
30 #define GETPOT_NAMESPACE QUESO // So we don't clash with other getpots
31 #include <queso/getpot.h>
32 #undef GETPOT_NAMESPACE
33 #endif // QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
34
35 #include <queso/GPMSAOptions.h>
36
37 #include <queso/GslVector.h>
38 #include <queso/SimulationOutputMesh.h>
39
40
41 // ODV = option default value
42 #define UQ_GPMSA_HELP ""
43 #define UQ_GPMSA_MAX_SIMULATOR_BASIS_VECTORS_ODV 0
44 #define UQ_GPMSA_SIMULATOR_BASIS_VARIANCE_TO_CAPTURE 1.0
45 #define UQ_GPMSA_TRUNCATION_ERROR_PRECISION_SHAPE_ODV 5.0
46 #define UQ_GPMSA_TRUNCATION_ERROR_PRECISION_SCALE_ODV 200.0
47 #define UQ_GPMSA_EMULATOR_PRECISION_SHAPE_ODV 5.0
48 #define UQ_GPMSA_EMULATOR_PRECISION_SCALE_ODV 0.2
49 #define UQ_GPMSA_OBSERVATIONAL_PRECISION_SHAPE_ODV 5.0
50 #define UQ_GPMSA_OBSERVATIONAL_PRECISION_SCALE_ODV 0.2
51 #define UQ_GPMSA_EMULATOR_CORRELATION_STRENGTH_ALPHA_ODV 1.0
52 #define UQ_GPMSA_EMULATOR_CORRELATION_STRENGTH_BETA_ODV 0.1
53 #define UQ_GPMSA_DISCREPANCY_PRECISION_SHAPE_ODV 1.0
54 #define UQ_GPMSA_DISCREPANCY_PRECISION_SCALE_ODV 1e4
55 #define UQ_GPMSA_DISCREPANCY_CORRELATION_STRENGTH_ALPHA_ODV 1.0
56 #define UQ_GPMSA_DISCREPANCY_CORRELATION_STRENGTH_BETA_ODV 0.1
57 #define UQ_GPMSA_EMULATOR_DATA_PRECISION_SHAPE_ODV 3.0
58 #define UQ_GPMSA_EMULATOR_DATA_PRECISION_SCALE_ODV 333.333
59 #define UQ_GPMSA_OBSERVATIONAL_PRECISION_RIDGE 1e-4
60 #define UQ_GPMSA_OBSERVATIONAL_COVARIANCE_RIDGE 0.0
61 #define UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE 1.0
62 static const bool UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC = false;
63 #define UQ_GPMSA_GAUSSIAN_DISCREPANCY_SUPPORT_THRESHOLD 0.05
64
65 namespace { // Anonymous namespace for helper functions
66
67 template <typename V>
68 void min_max_update(V & min, V & max, const V & new_data,
69 const char * warn_on_update = NULL)
70 {
71 unsigned int dim = min.sizeGlobal();
72 queso_assert_equal_to(dim, max.sizeGlobal());
73 queso_assert_equal_to(dim, new_data.sizeGlobal());
74
75 for (unsigned int p=0; p != dim; ++p)
76 {
77 if (warn_on_update)
78 {
79 if (min[p] > new_data[p])
80 queso_warning("Experimental " << warn_on_update << " " <<
81 new_data[p] << " at index " << p <<
82 " is below minimum simulation " <<
83 warn_on_update << " " << min[p]);
84 if (max[p] < new_data[p])
85 queso_warning("Experimental " << warn_on_update << " " <<
86 new_data[p] << " at index " << p <<
87 " is above maximum simulation " <<
88 warn_on_update << " " << max[p]);
89 }
90 else
91 {
92 min[p] = std::min(min[p], new_data[p]);
93 max[p] = std::max(max[p], new_data[p]);
94 }
95 }
96 }
97
98
99 template <typename V>
mean_var_update(unsigned int & n,V & mean,V & var,const V & new_data)100 void mean_var_update(unsigned int & n, V & mean, V & var, const V & new_data)
101 {
102 ++n;
103
104 const V delta (new_data - mean);
105
106 V delta_n (delta);
107 delta_n /= n;
108
109 mean += delta_n;
110
111 V delta2 (new_data - mean);
112 delta2 *= delta;
113
114 var += delta2;
115 }
116
117 } // end anonymous namespace
118
119
120
121 namespace QUESO {
122
GPMSAOptions(const BaseEnvironment & env,const char * prefix)123 GPMSAOptions::GPMSAOptions(
124 const BaseEnvironment & env,
125 const char * prefix) :
126 options_have_been_used(false)
127 {
128 this->set_defaults();
129 this->parse(env, prefix);
130 }
131
132
GPMSAOptions()133 GPMSAOptions::GPMSAOptions()
134 :
135 m_env(NULL),
136 #ifndef QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
137 m_parser(new BoostInputOptionsParser()),
138 #endif // QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
139 options_have_been_used(false)
140 {
141 this->set_defaults();
142 this->set_prefix("");
143 }
144
145
146 void
set_prefix(const char * prefix)147 GPMSAOptions::set_prefix(const char * prefix)
148 {
149 queso_require(!options_have_been_used);
150
151 m_prefix = std::string(prefix) + "gpmsa_";
152
153 m_option_help = m_prefix + "help";
154 m_option_maxEmulatorBasisVectors = m_prefix + "max_emulator_basis_vectors";
155 m_option_emulatorBasisVarianceToCapture = m_prefix + "emulator_basis_variance_to_capture";
156 m_option_truncationErrorPrecisionShape = m_prefix + "truncation_error_precision_shape";
157 m_option_truncationErrorPrecisionScale = m_prefix + "truncation_error_precision_scale";
158 m_option_emulatorPrecisionShape = m_prefix + "emulator_precision_shape";
159 m_option_emulatorPrecisionScale = m_prefix + "emulator_precision_scale";
160 m_option_calibrateObservationalPrecision = m_prefix + "calibrate_observational_precision";
161 m_option_observationalPrecisionShape = m_prefix + "observational_precision_shape";
162 m_option_observationalPrecisionScale = m_prefix + "observational_precision_scale";
163 m_option_emulatorCorrelationStrengthAlpha = m_prefix + "emulator_correlation_strength_alpha";
164 m_option_emulatorCorrelationStrengthBeta = m_prefix + "emulator_correlation_strength_beta";
165 m_option_discrepancyPrecisionShape = m_prefix + "discrepancy_precision_shape";
166 m_option_discrepancyPrecisionScale = m_prefix + "discrepancy_precision_scale";
167 m_option_discrepancyCorrelationStrengthAlpha = m_prefix + "discrepancy_correlation_strength_alpha";
168 m_option_discrepancyCorrelationStrengthBeta = m_prefix + "discrepancy_correlation_strength_beta";
169 m_option_emulatorDataPrecisionShape = m_prefix + "emulator_data_precision_shape";
170 m_option_emulatorDataPrecisionScale = m_prefix + "emulator_data_precision_scale";
171 m_option_observationalPrecisionRidge = m_prefix + "observational_precision_ridge";
172 m_option_observationalCovarianceRidge = m_prefix + "observational_covariance_ridge";
173 m_option_autoscaleMinMaxAll = m_prefix + "autoscale_min_max_all";
174 m_option_autoscaleMeanVarAll = m_prefix + "autoscale_mean_var_all";
175 m_option_gaussianDiscrepancyDistanceX = m_prefix + "gaussian_discrepancy_distance_x";
176 m_option_gaussianDiscrepancyDistanceY = m_prefix + "gaussian_discrepancy_distance_y";
177 m_option_gaussianDiscrepancyDistanceZ = m_prefix + "gaussian_discrepancy_distance_z";
178 m_option_gaussianDiscrepancyDistanceT = m_prefix + "gaussian_discrepancy_distance_t";
179 m_option_gaussianDiscrepancyPeriodicX = m_prefix + "gaussian_discrepancy_periodic_x";
180 m_option_gaussianDiscrepancyPeriodicY = m_prefix + "gaussian_discrepancy_periodic_y";
181 m_option_gaussianDiscrepancyPeriodicZ = m_prefix + "gaussian_discrepancy_periodic_z";
182 m_option_gaussianDiscrepancyPeriodicT = m_prefix + "gaussian_discrepancy_periodic_t";
183 m_option_gaussianDiscrepancySupportThreshold = m_prefix + "gaussian_discrepancy_support_threshold";
184 }
185
186
187
188 void
set_defaults()189 GPMSAOptions::set_defaults()
190 {
191 queso_require(!options_have_been_used);
192
193 m_help = UQ_GPMSA_HELP;
194 m_maxEmulatorBasisVectors = UQ_GPMSA_MAX_SIMULATOR_BASIS_VECTORS_ODV;
195 m_emulatorBasisVarianceToCapture = UQ_GPMSA_SIMULATOR_BASIS_VARIANCE_TO_CAPTURE;
196 m_truncationErrorPrecisionShape = UQ_GPMSA_TRUNCATION_ERROR_PRECISION_SHAPE_ODV;
197 m_truncationErrorPrecisionScale = UQ_GPMSA_TRUNCATION_ERROR_PRECISION_SCALE_ODV;
198 m_emulatorPrecisionShape = UQ_GPMSA_EMULATOR_PRECISION_SHAPE_ODV;
199 m_emulatorPrecisionScale = UQ_GPMSA_EMULATOR_PRECISION_SCALE_ODV;
200 m_calibrateObservationalPrecision = false;
201 m_observationalPrecisionShape = UQ_GPMSA_OBSERVATIONAL_PRECISION_SHAPE_ODV;
202 m_observationalPrecisionScale = UQ_GPMSA_OBSERVATIONAL_PRECISION_SCALE_ODV;
203 m_emulatorCorrelationStrengthAlpha = UQ_GPMSA_EMULATOR_CORRELATION_STRENGTH_ALPHA_ODV;
204 m_emulatorCorrelationStrengthBeta = UQ_GPMSA_EMULATOR_CORRELATION_STRENGTH_BETA_ODV;
205 m_discrepancyPrecisionShape = UQ_GPMSA_DISCREPANCY_PRECISION_SHAPE_ODV;
206 m_discrepancyPrecisionScale = UQ_GPMSA_DISCREPANCY_PRECISION_SCALE_ODV;
207 m_discrepancyCorrelationStrengthAlpha = UQ_GPMSA_DISCREPANCY_CORRELATION_STRENGTH_ALPHA_ODV;
208 m_discrepancyCorrelationStrengthBeta = UQ_GPMSA_DISCREPANCY_CORRELATION_STRENGTH_BETA_ODV;
209 m_emulatorDataPrecisionShape = UQ_GPMSA_EMULATOR_DATA_PRECISION_SHAPE_ODV;
210 m_emulatorDataPrecisionScale = UQ_GPMSA_EMULATOR_DATA_PRECISION_SCALE_ODV;
211 m_observationalPrecisionRidge = UQ_GPMSA_OBSERVATIONAL_PRECISION_RIDGE;
212 m_observationalCovarianceRidge = UQ_GPMSA_OBSERVATIONAL_COVARIANCE_RIDGE;
213
214 m_autoscaleMinMaxAll = false;
215 m_autoscaleMeanVarAll = false;
216
217 m_gaussianDiscrepancySupportThreshold = UQ_GPMSA_GAUSSIAN_DISCREPANCY_SUPPORT_THRESHOLD;
218
219 checkOptions();
220 }
221
222
223 void
parse(const BaseEnvironment & env,const char * prefix)224 GPMSAOptions::parse(const BaseEnvironment & env,
225 const char * prefix)
226 {
227 queso_require(!options_have_been_used);
228
229 m_env = &env;
230
231 if (m_env->optionsInputFileName() == "") {
232 queso_error_msg("Missing input file is required");
233 }
234
235 this->set_prefix(prefix);
236
237 #ifndef QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
238 m_parser.reset(new BoostInputOptionsParser(env.optionsInputFileName()));
239
240 m_parser->registerOption<std::string>
241 (m_option_help,
242 m_help,
243 "produce help message Gaussian process emulator");
244
245 m_parser->registerOption
246 (m_option_truncationErrorPrecisionShape,
247 m_truncationErrorPrecisionShape,
248 "shape hyperprior (Gamma) parameter for truncation error precision");
249 m_parser->registerOption
250 (m_option_truncationErrorPrecisionScale,
251 m_truncationErrorPrecisionScale,
252 "scale hyperprior (Gamma) parameter for truncation error precision");
253
254 m_parser->registerOption
255 (m_option_emulatorPrecisionShape,
256 m_emulatorPrecisionShape,
257 "shape hyperprior (Gamma) parameter for emulator precision");
258 m_parser->registerOption
259 (m_option_emulatorPrecisionScale,
260 m_emulatorPrecisionScale,
261 "scale hyperprior (Gamma) parameter for emulator precision");
262
263 m_parser->registerOption
264 (m_option_calibrateObservationalPrecision,
265 m_calibrateObservationalPrecision,
266 "whether to use a calibrated hyperparameter for observational precision");
267
268 m_parser->registerOption
269 (m_option_observationalPrecisionShape,
270 m_observationalPrecisionShape,
271 "shape hyperprior (Gamma) parameter for observational precision");
272 m_parser->registerOption
273 (m_option_observationalPrecisionScale,
274 m_observationalPrecisionScale,
275 "scale hyperprior (Gamma) parameter for observational precision");
276
277 m_parser->registerOption
278 (m_option_emulatorCorrelationStrengthAlpha,
279 m_emulatorCorrelationStrengthAlpha,
280 "alpha hyperprior (Beta) parameter for emulator correlation strength");
281 m_parser->registerOption
282 (m_option_emulatorCorrelationStrengthBeta,
283 m_emulatorCorrelationStrengthBeta,
284 "beta hyperprior (Beta) parameter for emulator correlation strength");
285
286 m_parser->registerOption
287 (m_option_discrepancyPrecisionShape,
288 m_discrepancyPrecisionShape,
289 "shape hyperprior (Gamma) parameter for discrepancy precision");
290 m_parser->registerOption
291 (m_option_discrepancyPrecisionScale,
292 m_discrepancyPrecisionScale,
293 "scale hyperprior (Gamma) parameter for discrepancy precision");
294
295 m_parser->registerOption
296 (m_option_discrepancyCorrelationStrengthAlpha,
297 m_discrepancyCorrelationStrengthAlpha,
298 "alpha hyperprior (Beta) parameter for discrepancy correlation strength");
299 m_parser->registerOption
300 (m_option_discrepancyCorrelationStrengthBeta,
301 m_discrepancyCorrelationStrengthBeta,
302 "beta hyperprior (Beta) parameter for discrepancy correlation strength");
303
304 m_parser->registerOption
305 (m_option_emulatorDataPrecisionShape,
306 m_emulatorDataPrecisionShape,
307 "shape hyperprior (Gamma) parameter for emulator data precision");
308 m_parser->registerOption
309 (m_option_emulatorDataPrecisionScale,
310 m_emulatorDataPrecisionScale,
311 "scale hyperprior (Gamma) parameter for emulator data precision");
312
313 m_parser->registerOption
314 (m_option_observationalPrecisionRidge,
315 m_observationalPrecisionRidge,
316 "ridge to add to observational precision matrix");
317
318 m_parser->registerOption
319 (m_option_observationalCovarianceRidge,
320 m_observationalCovarianceRidge,
321 "ridge to add to observational covariance matrix");
322
323 m_parser->registerOption
324 (m_option_autoscaleMinMaxAll,
325 m_autoscaleMinMaxAll,
326 "option to autoscale all parameters and outputs based on data range");
327 m_parser->registerOption
328 (m_option_autoscaleMeanVarAll,
329 m_autoscaleMeanVarAll,
330 "option to autoscale all parameters and outputs based on data statistics");
331
332 std::ostringstream defaultDiscrepancyDistanceString;
333 defaultDiscrepancyDistanceString << UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE;
334
335 m_parser->registerOption
336 (m_option_gaussianDiscrepancyDistanceX,
337 defaultDiscrepancyDistanceString.str(),
338 "x distance between neighboring gaussian discrepancy kernels");
339
340 m_parser->registerOption
341 (m_option_gaussianDiscrepancyDistanceY,
342 defaultDiscrepancyDistanceString.str(),
343 "y distance between neighboring gaussian discrepancy kernels");
344
345 m_parser->registerOption
346 (m_option_gaussianDiscrepancyDistanceZ,
347 defaultDiscrepancyDistanceString.str(),
348 "z distance between neighboring gaussian discrepancy kernels");
349
350 m_parser->registerOption
351 (m_option_gaussianDiscrepancyDistanceT,
352 defaultDiscrepancyDistanceString.str(),
353 "t distance between neighboring gaussian discrepancy kernels");
354
355 std::ostringstream defaultDiscrepancyPeriodicString;
356 defaultDiscrepancyPeriodicString << std::boolalpha << UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC;
357
358 m_parser->registerOption
359 (m_option_gaussianDiscrepancyPeriodicX,
360 defaultDiscrepancyPeriodicString.str(),
361 "whether gaussian discrepancy kernels are periodic in x");
362
363 m_parser->registerOption
364 (m_option_gaussianDiscrepancyPeriodicY,
365 defaultDiscrepancyPeriodicString.str(),
366 "whether gaussian discrepancy kernels are periodic in y");
367
368 m_parser->registerOption
369 (m_option_gaussianDiscrepancyPeriodicZ,
370 defaultDiscrepancyPeriodicString.str(),
371 "whether gaussian discrepancy kernels are periodic in z");
372
373 m_parser->registerOption
374 (m_option_gaussianDiscrepancyPeriodicT,
375 defaultDiscrepancyPeriodicString.str(),
376 "whether gaussian discrepancy kernels are periodic in t");
377
378 m_parser->registerOption
379 (m_option_gaussianDiscrepancySupportThreshold,
380 m_gaussianDiscrepancySupportThreshold,
381 "threshold below which to omit out-of-domain centered gaussian kernels");
382
383
384 m_parser->registerOption
385 (m_option_maxEmulatorBasisVectors,
386 m_maxEmulatorBasisVectors,
387 "max number of basis vectors to use in SVD of simulation output");
388
389 m_parser->scanInputFile();
390
391 m_parser->getOption<std::string>(m_option_help, m_help);
392 m_parser->getOption<double>(m_option_truncationErrorPrecisionShape, m_truncationErrorPrecisionShape);
393 m_parser->getOption<double>(m_option_truncationErrorPrecisionScale, m_truncationErrorPrecisionScale);
394 m_parser->getOption<double>(m_option_emulatorPrecisionShape, m_emulatorPrecisionShape);
395 m_parser->getOption<double>(m_option_emulatorPrecisionScale, m_emulatorPrecisionScale);
396 m_parser->getOption<bool> (m_option_calibrateObservationalPrecision, m_calibrateObservationalPrecision);
397 m_parser->getOption<double>(m_option_observationalPrecisionShape, m_observationalPrecisionShape);
398 m_parser->getOption<double>(m_option_observationalPrecisionScale, m_observationalPrecisionScale);
399 m_parser->getOption<double>(m_option_emulatorCorrelationStrengthAlpha, m_emulatorCorrelationStrengthAlpha);
400 m_parser->getOption<double>(m_option_emulatorCorrelationStrengthBeta, m_emulatorCorrelationStrengthBeta);
401 m_parser->getOption<double>(m_option_discrepancyPrecisionShape, m_discrepancyPrecisionShape);
402 m_parser->getOption<double>(m_option_discrepancyPrecisionScale, m_discrepancyPrecisionScale);
403 m_parser->getOption<double>(m_option_discrepancyCorrelationStrengthAlpha, m_discrepancyCorrelationStrengthAlpha);
404 m_parser->getOption<double>(m_option_discrepancyCorrelationStrengthBeta, m_discrepancyCorrelationStrengthBeta);
405 m_parser->getOption<double>(m_option_emulatorDataPrecisionShape, m_emulatorDataPrecisionShape);
406 m_parser->getOption<double>(m_option_emulatorDataPrecisionScale, m_emulatorDataPrecisionScale);
407 m_parser->getOption<double>(m_option_observationalPrecisionRidge, m_observationalPrecisionRidge);
408 m_parser->getOption<double>(m_option_observationalCovarianceRidge, m_observationalCovarianceRidge);
409 m_parser->getOption<bool> (m_option_autoscaleMinMaxAll, m_autoscaleMinMaxAll);
410 m_parser->getOption<bool> (m_option_autoscaleMeanVarAll, m_autoscaleMeanVarAll);
411 m_parser->getOption<int> (m_option_maxEmulatorBasisVectors, m_maxEmulatorBasisVectors);
412 m_parser->getOption<std::vector<double> >(m_option_gaussianDiscrepancyDistanceX, m_gaussianDiscrepancyDistanceX);
413 m_parser->getOption<std::vector<double> >(m_option_gaussianDiscrepancyDistanceY, m_gaussianDiscrepancyDistanceY);
414 m_parser->getOption<std::vector<double> >(m_option_gaussianDiscrepancyDistanceZ, m_gaussianDiscrepancyDistanceZ);
415 m_parser->getOption<std::vector<double> >(m_option_gaussianDiscrepancyDistanceT, m_gaussianDiscrepancyDistanceT);
416
417 // I cannot get boost::program_options to play nicely with
418 // vector<bool>, and we don't yet implement periodic discrepancy
419 // options, so let's just ignore them here for now. By the time we
420 // implement them we'll have removed the deprecated boost::po code
421 // anyway.
422 /*
423 m_parser->getOption<std::vector<bool> > (m_option_gaussianDiscrepancyPeriodicX, m_gaussianDiscrepancyPeriodicX);
424 m_parser->getOption<std::vector<bool> > (m_option_gaussianDiscrepancyPeriodicY, m_gaussianDiscrepancyPeriodicY);
425 m_parser->getOption<std::vector<bool> > (m_option_gaussianDiscrepancyPeriodicZ, m_gaussianDiscrepancyPeriodicZ);
426 m_parser->getOption<std::vector<bool> > (m_option_gaussianDiscrepancyPeriodicT, m_gaussianDiscrepancyPeriodicT);
427 */
428 m_parser->getOption<double>(m_option_gaussianDiscrepancySupportThreshold, m_gaussianDiscrepancySupportThreshold);
429 #else
430 m_help = env.input()(m_option_help, UQ_GPMSA_HELP);
431
432 m_truncationErrorPrecisionShape =
433 env.input()(m_option_truncationErrorPrecisionShape,
434 m_truncationErrorPrecisionShape);
435 m_truncationErrorPrecisionScale =
436 env.input()(m_option_truncationErrorPrecisionScale,
437 m_truncationErrorPrecisionScale);
438
439 m_emulatorPrecisionShape =
440 env.input()(m_option_emulatorPrecisionShape,
441 m_emulatorPrecisionShape);
442 m_emulatorPrecisionScale =
443 env.input()(m_option_emulatorPrecisionScale,
444 m_emulatorPrecisionScale);
445
446 m_calibrateObservationalPrecision =
447 env.input()(m_option_calibrateObservationalPrecision,
448 m_calibrateObservationalPrecision);
449 m_observationalPrecisionShape =
450 env.input()(m_option_observationalPrecisionShape,
451 m_observationalPrecisionShape);
452 m_observationalPrecisionScale =
453 env.input()(m_option_observationalPrecisionScale,
454 m_observationalPrecisionScale);
455
456 m_emulatorCorrelationStrengthAlpha =
457 env.input()(m_option_emulatorCorrelationStrengthAlpha,
458 m_emulatorCorrelationStrengthAlpha);
459 m_emulatorCorrelationStrengthBeta =
460 env.input()(m_option_emulatorCorrelationStrengthBeta,
461 m_emulatorCorrelationStrengthBeta);
462
463 m_discrepancyPrecisionShape =
464 env.input()(m_option_discrepancyPrecisionShape,
465 m_discrepancyPrecisionShape);
466 m_discrepancyPrecisionScale =
467 env.input()(m_option_discrepancyPrecisionScale,
468 m_discrepancyPrecisionScale);
469
470 m_discrepancyCorrelationStrengthAlpha =
471 env.input()(m_option_discrepancyCorrelationStrengthAlpha,
472 m_discrepancyCorrelationStrengthAlpha);
473 m_discrepancyCorrelationStrengthBeta =
474 env.input()(m_option_discrepancyCorrelationStrengthBeta,
475 m_discrepancyCorrelationStrengthBeta);
476
477 m_emulatorDataPrecisionShape =
478 env.input()(m_option_emulatorDataPrecisionShape,
479 m_emulatorDataPrecisionShape);
480 m_emulatorDataPrecisionScale =
481 env.input()(m_option_emulatorDataPrecisionScale,
482 m_emulatorDataPrecisionScale);
483
484 m_observationalPrecisionRidge =
485 env.input()(m_option_observationalPrecisionRidge,
486 m_observationalPrecisionRidge);
487 m_observationalCovarianceRidge =
488 env.input()(m_option_observationalCovarianceRidge,
489 m_observationalCovarianceRidge);
490
491 m_autoscaleMinMaxAll =
492 env.input()(m_option_autoscaleMinMaxAll,
493 m_autoscaleMinMaxAll);
494 m_autoscaleMeanVarAll =
495 env.input()(m_option_autoscaleMeanVarAll,
496 m_autoscaleMeanVarAll);
497
498 for (unsigned int i = 0,
499 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyDistanceX);
500 i != size; ++i)
501 {
502 if (m_gaussianDiscrepancyDistanceX.size() <= i)
503 m_gaussianDiscrepancyDistanceX.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
504 m_gaussianDiscrepancyDistanceX[i] =
505 env.input()(m_option_gaussianDiscrepancyDistanceX,
506 m_gaussianDiscrepancyDistanceX[i]);
507 }
508
509 for (unsigned int i = 0,
510 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyDistanceY);
511 i != size; ++i)
512 {
513 if (m_gaussianDiscrepancyDistanceY.size() <= i)
514 m_gaussianDiscrepancyDistanceY.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
515 m_gaussianDiscrepancyDistanceY[i] =
516 env.input()(m_option_gaussianDiscrepancyDistanceY,
517 m_gaussianDiscrepancyDistanceY[i]);
518 }
519
520 for (unsigned int i = 0,
521 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyDistanceZ);
522 i != size; ++i)
523 {
524 if (m_gaussianDiscrepancyDistanceZ.size() <= i)
525 m_gaussianDiscrepancyDistanceZ.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
526 m_gaussianDiscrepancyDistanceZ[i] =
527 env.input()(m_option_gaussianDiscrepancyDistanceZ,
528 m_gaussianDiscrepancyDistanceZ[i]);
529 }
530
531 for (unsigned int i = 0,
532 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyDistanceT);
533 i != size; ++i)
534 {
535 if (m_gaussianDiscrepancyDistanceT.size() <= i)
536 m_gaussianDiscrepancyDistanceT.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
537 m_gaussianDiscrepancyDistanceT[i] =
538 env.input()(m_option_gaussianDiscrepancyDistanceT,
539 m_gaussianDiscrepancyDistanceT[i]);
540 }
541
542
543 for (unsigned int i = 0,
544 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyPeriodicX);
545 i != size; ++i)
546 {
547 if (m_gaussianDiscrepancyPeriodicX.size() <= i)
548 m_gaussianDiscrepancyPeriodicX.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
549 m_gaussianDiscrepancyPeriodicX[i] =
550 env.input()(m_option_gaussianDiscrepancyPeriodicX,
551 bool(m_gaussianDiscrepancyPeriodicX[i]));
552 }
553
554 for (unsigned int i = 0,
555 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyPeriodicY);
556 i != size; ++i)
557 {
558 if (m_gaussianDiscrepancyPeriodicY.size() <= i)
559 m_gaussianDiscrepancyPeriodicY.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
560 m_gaussianDiscrepancyPeriodicY[i] =
561 env.input()(m_option_gaussianDiscrepancyPeriodicY,
562 bool(m_gaussianDiscrepancyPeriodicY[i]));
563 }
564
565 for (unsigned int i = 0,
566 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyPeriodicZ);
567 i != size; ++i)
568 {
569 if (m_gaussianDiscrepancyPeriodicZ.size() <= i)
570 m_gaussianDiscrepancyPeriodicZ.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
571 m_gaussianDiscrepancyPeriodicZ[i] =
572 env.input()(m_option_gaussianDiscrepancyPeriodicZ,
573 bool(m_gaussianDiscrepancyPeriodicZ[i]));
574 }
575
576 for (unsigned int i = 0,
577 size = env.input().vector_variable_size(m_option_gaussianDiscrepancyPeriodicT);
578 i != size; ++i)
579 {
580 if (m_gaussianDiscrepancyPeriodicT.size() <= i)
581 m_gaussianDiscrepancyPeriodicT.push_back(UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
582 m_gaussianDiscrepancyPeriodicT[i] =
583 env.input()(m_option_gaussianDiscrepancyPeriodicT,
584 bool(m_gaussianDiscrepancyPeriodicT[i]));
585 }
586
587 m_gaussianDiscrepancySupportThreshold =
588 env.input()(m_option_gaussianDiscrepancySupportThreshold,
589 m_gaussianDiscrepancySupportThreshold);
590
591 m_maxEmulatorBasisVectors =
592 env.input()(m_option_maxEmulatorBasisVectors,
593 m_maxEmulatorBasisVectors);
594 #endif // QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
595
596 checkOptions();
597 }
598
~GPMSAOptions()599 GPMSAOptions::~GPMSAOptions()
600 {
601 }
602
603
604 void
set_autoscale_minmax()605 GPMSAOptions::set_autoscale_minmax()
606 {
607 queso_require(!options_have_been_used);
608
609 this->m_autoscaleMinMaxAll = true;
610 }
611
612
613 void
set_autoscale_minmax_uncertain_parameter(unsigned int i)614 GPMSAOptions::set_autoscale_minmax_uncertain_parameter(unsigned int i)
615 {
616 queso_require(!options_have_been_used);
617
618 m_autoscaleMinMaxUncertain.insert(i);
619 }
620
621
622 void
set_autoscale_minmax_scenario_parameter(unsigned int i)623 GPMSAOptions::set_autoscale_minmax_scenario_parameter(unsigned int i)
624 {
625 queso_require(!options_have_been_used);
626
627 m_autoscaleMinMaxScenario.insert(i);
628 }
629
630
631 void
set_autoscale_minmax_output(unsigned int i)632 GPMSAOptions::set_autoscale_minmax_output(unsigned int i)
633 {
634 m_autoscaleMinMaxOutput.insert(i);
635 }
636
637
638 void
set_autoscale_meanvar()639 GPMSAOptions::set_autoscale_meanvar()
640 {
641 queso_require(!options_have_been_used);
642
643 this->m_autoscaleMeanVarAll = true;
644 }
645
646
647 void
set_autoscale_meanvar_uncertain_parameter(unsigned int i)648 GPMSAOptions::set_autoscale_meanvar_uncertain_parameter(unsigned int i)
649 {
650 queso_require(!options_have_been_used);
651
652 m_autoscaleMeanVarUncertain.insert(i);
653 }
654
655
656 void
set_autoscale_meanvar_scenario_parameter(unsigned int i)657 GPMSAOptions::set_autoscale_meanvar_scenario_parameter(unsigned int i)
658 {
659 queso_require(!options_have_been_used);
660
661 m_autoscaleMeanVarScenario.insert(i);
662 }
663
664
665 void
set_autoscale_meanvar_output(unsigned int i)666 GPMSAOptions::set_autoscale_meanvar_output(unsigned int i)
667 {
668 m_autoscaleMeanVarOutput.insert(i);
669 }
670
671
672 void
set_uncertain_parameter_scaling(unsigned int i,double range_min,double range_max)673 GPMSAOptions::set_uncertain_parameter_scaling(unsigned int i,
674 double range_min,
675 double range_max)
676 {
677 queso_require(!options_have_been_used);
678
679 if (i >= m_uncertainScaleMin.size())
680 {
681 m_uncertainScaleMin.resize(i+1, 0);
682 m_uncertainScaleRange.resize(i+1, 1);
683 }
684 m_uncertainScaleMin[i] = range_min;
685 m_uncertainScaleRange[i] = range_max - range_min;
686 }
687
688
689 void
set_scenario_parameter_scaling(unsigned int i,double range_min,double range_max)690 GPMSAOptions::set_scenario_parameter_scaling(unsigned int i,
691 double range_min,
692 double range_max)
693 {
694 queso_require(!options_have_been_used);
695
696 if (i >= m_scenarioScaleMin.size())
697 {
698 m_scenarioScaleMin.resize(i+1, 0);
699 m_scenarioScaleRange.resize(i+1, 1);
700 }
701 m_scenarioScaleMin[i] = range_min;
702 m_scenarioScaleRange[i] = range_max - range_min;
703 }
704
705
706 void
set_output_scaling(unsigned int i,double range_min,double range_max)707 GPMSAOptions::set_output_scaling(unsigned int i,
708 double range_min,
709 double range_max)
710 {
711 if (i >= m_outputScaleMin.size())
712 {
713 m_outputScaleMin.resize(i+1, 0);
714 m_outputScaleRange.resize(i+1, 1);
715 }
716 m_outputScaleMin[i] = range_min;
717 m_outputScaleRange[i] = range_max - range_min;
718 }
719
720
721 template <typename V>
722 void
set_final_scaling(const std::vector<typename SharedPtr<V>::Type> & m_simulationScenarios,const std::vector<typename SharedPtr<V>::Type> & m_simulationParameters,const std::vector<typename SharedPtr<V>::Type> & m_simulationOutputs,const std::vector<typename SharedPtr<V>::Type> & m_experimentScenarios,const std::vector<typename SharedPtr<V>::Type> & m_experimentOutputs,const std::vector<typename SharedPtr<SimulationOutputMesh<V>>::Type> & m_simulationMeshes)723 GPMSAOptions::set_final_scaling
724 (const std::vector<typename SharedPtr<V>::Type> & m_simulationScenarios,
725 const std::vector<typename SharedPtr<V>::Type> & m_simulationParameters,
726 const std::vector<typename SharedPtr<V>::Type> & m_simulationOutputs,
727 const std::vector<typename SharedPtr<V>::Type> & m_experimentScenarios,
728 const std::vector<typename SharedPtr<V>::Type> & m_experimentOutputs,
729 const std::vector<typename SharedPtr<SimulationOutputMesh<V> >::Type> & m_simulationMeshes)
730 {
731 if ((m_autoscaleMinMaxAll && m_autoscaleMeanVarAll) ||
732 ((m_autoscaleMinMaxAll || m_autoscaleMeanVarAll) &&
733 (!m_autoscaleMinMaxUncertain.empty() ||
734 !m_autoscaleMeanVarUncertain.empty() ||
735 !m_autoscaleMinMaxScenario.empty() ||
736 !m_autoscaleMeanVarScenario.empty() ||
737 !m_autoscaleMinMaxOutput.empty() ||
738 !m_autoscaleMeanVarOutput.empty())))
739 queso_error_msg("Cannot autoscale based on incompatible criteria");
740
741 unsigned int dimScenario = m_simulationScenarios.size() ?
742 m_simulationScenarios[0]->sizeGlobal() : 0;
743
744 unsigned int dimParameter = m_simulationParameters.size() ?
745 m_simulationParameters[0]->sizeGlobal() : 0;
746
747 unsigned int dimOutput = m_simulationOutputs.size() ?
748 m_simulationOutputs[0]->sizeGlobal() : 0;
749
750 queso_require_msg(dimOutput, "GPMSA needs *some* output data");
751
752 if (m_autoscaleMinMaxAll || !m_autoscaleMinMaxScenario.empty())
753 {
754 V maxScenario(*m_simulationScenarios[0]);
755
756 V minScenario(*m_simulationScenarios[0]);
757
758 // Only use simulation data for scaling.
759 for (unsigned int i=1; i < m_simulationScenarios.size(); ++i)
760 min_max_update(minScenario, maxScenario,
761 (*m_simulationScenarios[i]));
762
763 // But check the experimental data, and yell at the user if it
764 // falls outside the simulation data range, because they
765 // probably don't want to be extrapolating with their GP.
766 for (unsigned int i=0; i < m_experimentScenarios.size(); ++i)
767 min_max_update(minScenario, maxScenario,
768 (*m_experimentScenarios[i]),
769 /*warn_on_update =*/ "scenario parameter");
770
771 for (unsigned int p=0; p != dimScenario; ++p)
772 if (m_autoscaleMinMaxAll ||
773 m_autoscaleMinMaxScenario.count(p))
774 {
775 if ((m_scenarioScaleMin.size() > p) &&
776 ((m_scenarioScaleMin[p] != 0) ||
777 (m_scenarioScaleRange[p] != 1)))
778 queso_error_msg("Cannot autoscale and manually scale the same scenario parameter");
779
780 this->set_scenario_parameter_scaling(p, minScenario[p],
781 maxScenario[p]);
782 }
783 }
784
785 if (m_autoscaleMinMaxAll || !m_autoscaleMinMaxUncertain.empty())
786 {
787 V maxUncertain(*m_simulationParameters[0]);
788
789 V minUncertain(*m_simulationParameters[0]);
790
791 for (unsigned int i=1; i < m_simulationParameters.size(); ++i)
792 min_max_update(minUncertain, maxUncertain,
793 (*m_simulationParameters[i]));
794
795 for (unsigned int p=0; p != dimParameter; ++p)
796 if (m_autoscaleMinMaxAll ||
797 m_autoscaleMinMaxUncertain.count(p))
798 {
799 if ((m_uncertainScaleMin.size() > p) &&
800 ((m_uncertainScaleMin[p] != 0) ||
801 (m_uncertainScaleRange[p] != 1)))
802 queso_error_msg("Cannot autoscale and manually scale the same uncertain parameter");
803
804 this->set_uncertain_parameter_scaling(p, minUncertain[p],
805 maxUncertain[p]);
806 }
807 }
808
809 if (m_autoscaleMinMaxAll || !m_autoscaleMinMaxOutput.empty())
810 {
811 V maxOutput(*m_simulationOutputs[0]);
812
813 V minOutput(*m_simulationOutputs[0]);
814
815 // Only use simulation data for scaling. In this case we won't
816 // warn if the experimental data falls outside the simulation
817 // data bounds, because the user can't control their outputs,
818 // just their inputs.
819 for (unsigned int i=1; i < m_simulationOutputs.size(); ++i)
820 min_max_update(minOutput, maxOutput,
821 (*m_simulationOutputs[i]));
822
823 // Scale together any functional output data, data that is
824 // described by a simulation mesh, then scale multivariate
825 // outputs independently.
826 unsigned int first_multivariate_index = 0;
827 for (unsigned int m=0; m != m_simulationMeshes.size(); ++m)
828 {
829 const SimulationOutputMesh<V> & mesh = *m_simulationMeshes[m];
830 const unsigned int mesh_n_outputs = mesh.n_outputs();
831 queso_assert_greater(mesh_n_outputs, 0);
832 queso_assert_equal_to(mesh.first_solution_index(),
833 first_multivariate_index);
834 first_multivariate_index += mesh_n_outputs;
835
836 if (m_autoscaleMinMaxAll ||
837 m_autoscaleMinMaxOutput.count(m))
838 {
839 if ((m_outputScaleMin.size() > m) &&
840 ((m_outputScaleMin[m] != 0) ||
841 (m_outputScaleRange[m] != 1)))
842 queso_error_msg("Cannot autoscale and manually scale the same output data");
843
844 // FIXME - this will break if/when we introduce
845 // non-Lagrange solution bases.
846 double full_min = minOutput[mesh.first_solution_index()],
847 full_max = maxOutput[mesh.first_solution_index()];
848 for (unsigned int i=1; i != mesh_n_outputs; ++i)
849 {
850 full_min = std::min(full_min, minOutput[mesh.first_solution_index()+i]);
851 full_max = std::max(full_max, maxOutput[mesh.first_solution_index()+i]);
852 }
853 this->set_output_scaling(m, full_min, full_max);
854 }
855 }
856
857 const unsigned int n_multivariate_indices =
858 dimOutput - first_multivariate_index;
859 const unsigned int n_variables = m_simulationMeshes.size() + n_multivariate_indices;
860 const unsigned int variable_index_to_output_index = first_multivariate_index - m_simulationMeshes.size();
861 for (unsigned int p=m_simulationMeshes.size(); p != n_variables; ++p)
862 {
863 if (m_autoscaleMinMaxAll ||
864 m_autoscaleMinMaxOutput.count(p))
865 {
866 if ((m_outputScaleMin.size() > p) &&
867 ((m_outputScaleMin[p] != 0) ||
868 (m_outputScaleRange[p] != 1)))
869 queso_error_msg("Cannot autoscale and manually scale the same output data");
870
871 this->set_output_scaling
872 (p, minOutput[p+variable_index_to_output_index],
873 maxOutput[p+variable_index_to_output_index]);
874 }
875 }
876 }
877
878
879 if (m_autoscaleMeanVarAll || !m_autoscaleMeanVarScenario.empty())
880 {
881 unsigned int n=0;
882
883 V meanScenario(*m_simulationScenarios[0]);
884
885 V varScenario(m_simulationScenarios[0]->env(),
886 m_simulationScenarios[0]->map());
887
888 // For consistency with the min-max behavior, only normalize
889 // using simulation data, not experiment data.
890 for (unsigned int i=0; i < m_simulationScenarios.size(); ++i)
891 mean_var_update(n, meanScenario, varScenario,
892 *m_simulationScenarios[i]);
893
894 varScenario /= n - 1;
895
896 for (unsigned int p=0; p != dimScenario; ++p)
897 if (m_autoscaleMeanVarAll ||
898 m_autoscaleMeanVarScenario.count(p))
899 {
900 if ((m_scenarioScaleMin.size() > p) &&
901 ((m_scenarioScaleMin[p] != 0) ||
902 (m_scenarioScaleRange[p] != 1)))
903 queso_error_msg("Cannot autoscale and manually scale the same scenario parameter");
904
905 this->set_scenario_parameter_scaling(p, meanScenario[p],
906 meanScenario[p] +
907 std::sqrt(varScenario[p]));
908 }
909 }
910
911
912 if (m_autoscaleMeanVarAll || !m_autoscaleMeanVarUncertain.empty())
913 {
914 unsigned int n=0;
915
916 V meanUncertain(*m_simulationParameters[0]);
917
918 V varUncertain(m_simulationParameters[0]->env(),
919 m_simulationParameters[0]->map());
920
921 for (unsigned int i=0; i < m_simulationParameters.size(); ++i)
922 mean_var_update(n, meanUncertain, varUncertain,
923 *m_simulationParameters[i]);
924
925 varUncertain /= n - 1;
926
927 for (unsigned int p=0; p != dimScenario; ++p)
928 if (m_autoscaleMeanVarAll ||
929 m_autoscaleMeanVarUncertain.count(p))
930 {
931 if ((m_uncertainScaleMin.size() > p) &&
932 ((m_uncertainScaleMin[p] != 0) ||
933 (m_uncertainScaleRange[p] != 1)))
934 queso_error_msg("Cannot autoscale and manually scale the same uncertain parameter");
935
936 this->set_uncertain_parameter_scaling(p, meanUncertain[p],
937 meanUncertain[p] +
938 std::sqrt(varUncertain[p]));
939 }
940 }
941
942 if (m_autoscaleMeanVarAll || !m_autoscaleMeanVarOutput.empty())
943 {
944 unsigned int n=0;
945
946 V meanOutput(*m_simulationOutputs[0]);
947
948 V varOutput(m_simulationOutputs[0]->env(),
949 m_simulationOutputs[0]->map());
950
951 // This gives us the right mean and var values for any
952 // multivariate components of the problem; we'll handle
953 // functional components separately shortly.
954 for (unsigned int i=0; i < m_simulationOutputs.size(); ++i)
955 mean_var_update(n, meanOutput, varOutput,
956 *m_simulationOutputs[i]);
957
958 varOutput /= n - 1;
959
960 // Scale together any functional output data, data that is
961 // described by a simulation mesh, then scale multivariate
962 // outputs independently.
963 unsigned int first_multivariate_index = 0;
964 for (unsigned int m=0; m != m_simulationMeshes.size(); ++m)
965 {
966 const SimulationOutputMesh<V> & mesh = *m_simulationMeshes[m];
967 const unsigned int mesh_n_outputs = mesh.n_outputs();
968 queso_assert_greater(mesh_n_outputs, 0);
969 queso_assert_equal_to(mesh.first_solution_index(),
970 first_multivariate_index);
971 first_multivariate_index += mesh_n_outputs;
972
973 if (m_autoscaleMeanVarAll ||
974 m_autoscaleMeanVarOutput.count(m))
975 {
976 if ((m_outputScaleMin.size() > m) &&
977 ((m_outputScaleMin[m] != 0) ||
978 (m_outputScaleRange[m] != 1)))
979 queso_error_msg("Cannot autoscale and manually scale the same output data");
980
981 unsigned int n = 0;
982 double full_mean = 0, full_var = 0;
983 for (unsigned int i=0; i < m_simulationOutputs.size(); ++i)
984 for (unsigned int j=0; j != mesh_n_outputs; ++j)
985 mean_var_update
986 (n, full_mean, full_var,
987 (*m_simulationOutputs[i])[mesh.first_solution_index()+i]);
988
989 full_var /= (n - 1);
990
991 this->set_output_scaling
992 (m, full_mean, full_mean + std::sqrt(full_var));
993 }
994 }
995
996 const unsigned int n_multivariate_indices =
997 dimOutput - first_multivariate_index;
998 const unsigned int n_variables = m_simulationMeshes.size() + n_multivariate_indices;
999 const unsigned int variable_index_to_output_index = first_multivariate_index - m_simulationMeshes.size();
1000 for (unsigned int p=m_simulationMeshes.size(); p != n_variables; ++p)
1001 if (m_autoscaleMeanVarAll ||
1002 m_autoscaleMeanVarOutput.count(p))
1003 {
1004 if ((m_outputScaleMin.size() > p) &&
1005 ((m_outputScaleMin[p] != 0) ||
1006 (m_outputScaleRange[p] != 1)))
1007 queso_error_msg("Cannot autoscale and manually scale the same output data");
1008
1009 this->set_output_scaling
1010 (p, meanOutput[p+variable_index_to_output_index],
1011 meanOutput[p+variable_index_to_output_index] +
1012 std::sqrt(varOutput[p+variable_index_to_output_index]));
1013 }
1014 }
1015
1016 m_output_index_to_variable_index.resize(dimOutput);
1017
1018 unsigned int next_variable_begin = 0;
1019 for (unsigned int m=0; m != m_simulationMeshes.size(); ++m)
1020 {
1021 const SimulationOutputMesh<V> & mesh = *m_simulationMeshes[m];
1022 const unsigned int mesh_n_outputs = mesh.n_outputs();
1023 next_variable_begin += mesh_n_outputs;
1024 for (unsigned int i = mesh.first_solution_index();
1025 i != next_variable_begin; ++i)
1026 m_output_index_to_variable_index[i] = m;
1027 }
1028
1029 unsigned int next_mv_index = m_simulationMeshes.size();
1030 for (unsigned int i = next_variable_begin; i != dimOutput; ++i)
1031 m_output_index_to_variable_index[i] = next_mv_index++;
1032
1033 // Make sure we have gaussian discrepancy option values available
1034 // for every mesh, even if they weren't set by the input file.
1035 const unsigned int n_meshes = m_simulationMeshes.size();
1036 m_gaussianDiscrepancyDistanceX.resize
1037 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
1038 m_gaussianDiscrepancyDistanceY.resize
1039 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
1040 m_gaussianDiscrepancyDistanceZ.resize
1041 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
1042 m_gaussianDiscrepancyDistanceT.resize
1043 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_DISTANCE);
1044 m_gaussianDiscrepancyPeriodicX.resize
1045 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
1046 m_gaussianDiscrepancyPeriodicY.resize
1047 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
1048 m_gaussianDiscrepancyPeriodicZ.resize
1049 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
1050 m_gaussianDiscrepancyPeriodicT.resize
1051 (n_meshes, UQ_GPMSA_GAUSSIAN_DISCREPANCY_PERIODIC);
1052
1053 options_have_been_used = true;
1054 }
1055
1056
1057 double
normalized_scenario_parameter(unsigned int i,double physical_param)1058 GPMSAOptions::normalized_scenario_parameter(unsigned int i,
1059 double physical_param)
1060 const
1061 {
1062 if (i < m_scenarioScaleMin.size())
1063 return (physical_param - m_scenarioScaleMin[i]) /
1064 (m_scenarioScaleRange[i] ? m_scenarioScaleRange[i] : 1);
1065 return physical_param;
1066 }
1067
1068
1069 double
normalized_uncertain_parameter(unsigned int i,double physical_param)1070 GPMSAOptions::normalized_uncertain_parameter(unsigned int i,
1071 double physical_param)
1072 const
1073 {
1074 if (i < m_uncertainScaleMin.size())
1075 return (physical_param - m_uncertainScaleMin[i]) /
1076 (m_uncertainScaleRange[i] ? m_uncertainScaleRange[i] : 1);
1077 return physical_param;
1078 }
1079
1080
1081 double
normalized_output(unsigned int i,double output_data)1082 GPMSAOptions::normalized_output(unsigned int i,
1083 double output_data)
1084 const
1085 {
1086 return this->normalized_output_variable
1087 (m_output_index_to_variable_index[i], output_data);
1088 }
1089
1090
1091
1092 double
normalized_output_variable(unsigned int i,double output_data)1093 GPMSAOptions::normalized_output_variable(unsigned int i,
1094 double output_data)
1095 const
1096 {
1097 if (i < m_outputScaleMin.size())
1098 return (output_data - m_outputScaleMin[i]) /
1099 (m_outputScaleRange[i] ? m_outputScaleRange[i] : 1);
1100 return output_data;
1101 }
1102
1103
1104
1105 double
output_scale(unsigned int i)1106 GPMSAOptions::output_scale(unsigned int i)
1107 const
1108 {
1109 return this->output_scale_variable
1110 (m_output_index_to_variable_index[i]);
1111 }
1112
1113
1114
1115 double
output_scale_variable(unsigned int i)1116 GPMSAOptions::output_scale_variable(unsigned int i)
1117 const
1118 {
1119 if (i < m_outputScaleRange.size() &&
1120 m_outputScaleRange[i])
1121 return m_outputScaleRange[i];
1122 return 1;
1123 }
1124
1125
1126
1127
1128 void
checkOptions()1129 GPMSAOptions::checkOptions()
1130 {
1131 if (m_help != "") {
1132 if (m_env && m_env->subDisplayFile()) {
1133 *m_env->subDisplayFile() << (*this) << std::endl;
1134 }
1135 }
1136 }
1137
1138 void
print(std::ostream & os)1139 GPMSAOptions::print(std::ostream& os) const
1140 {
1141 os << "\n" << m_option_truncationErrorPrecisionShape << " = " << this->m_truncationErrorPrecisionShape
1142 << "\n" << m_option_truncationErrorPrecisionScale << " = " << this->m_truncationErrorPrecisionScale
1143 << "\n" << m_option_emulatorPrecisionShape << " = " << this->m_emulatorPrecisionShape
1144 << "\n" << m_option_emulatorPrecisionScale << " = " << this->m_emulatorPrecisionScale
1145 << "\n" << m_option_calibrateObservationalPrecision << " = " << this->m_calibrateObservationalPrecision
1146 << "\n" << m_option_observationalPrecisionShape << " = " << this->m_observationalPrecisionShape
1147 << "\n" << m_option_observationalPrecisionScale << " = " << this->m_observationalPrecisionScale
1148 << "\n" << m_option_emulatorCorrelationStrengthAlpha << " = " << this->m_emulatorCorrelationStrengthAlpha
1149 << "\n" << m_option_emulatorCorrelationStrengthBeta << " = " << this->m_emulatorCorrelationStrengthBeta
1150 << "\n" << m_option_discrepancyPrecisionShape << " = " << this->m_discrepancyPrecisionShape
1151 << "\n" << m_option_discrepancyPrecisionScale << " = " << this->m_discrepancyPrecisionScale
1152 << "\n" << m_option_discrepancyCorrelationStrengthAlpha << " = " << this->m_discrepancyCorrelationStrengthAlpha
1153 << "\n" << m_option_discrepancyCorrelationStrengthBeta << " = " << this->m_discrepancyCorrelationStrengthBeta
1154 << "\n" << m_option_emulatorDataPrecisionShape << " = " << this->m_emulatorDataPrecisionShape
1155 << "\n" << m_option_emulatorDataPrecisionScale << " = " << this->m_emulatorDataPrecisionScale
1156 << "\n" << m_option_observationalPrecisionRidge << " = " << this->m_observationalPrecisionRidge
1157 << "\n" << m_option_observationalCovarianceRidge << " = " << this->m_observationalCovarianceRidge
1158 << "\n" << m_option_autoscaleMinMaxAll << " = " << this->m_autoscaleMinMaxAll
1159 << "\n" << m_option_autoscaleMeanVarAll << " = " << this->m_autoscaleMeanVarAll
1160 << "\n" << m_option_maxEmulatorBasisVectors << " = " << this->m_maxEmulatorBasisVectors;
1161
1162 os << "\n" << m_option_gaussianDiscrepancyDistanceX << " = {";
1163 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyDistanceX.size(); i != size; ++i)
1164 {
1165 os << this->m_gaussianDiscrepancyDistanceX[i];
1166 if (i+1 != size)
1167 os << ',';
1168 }
1169 os << '}';
1170 os << "\n" << m_option_gaussianDiscrepancyDistanceY << " = {";
1171 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyDistanceY.size(); i != size; ++i)
1172 {
1173 os << this->m_gaussianDiscrepancyDistanceY[i];
1174 if (i+1 != size)
1175 os << ',';
1176 }
1177 os << '}';
1178 os << "\n" << m_option_gaussianDiscrepancyDistanceZ << " = {";
1179 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyDistanceZ.size(); i != size; ++i)
1180 {
1181 os << this->m_gaussianDiscrepancyDistanceZ[i];
1182 if (i+1 != size)
1183 os << ',';
1184 }
1185 os << '}';
1186 os << "\n" << m_option_gaussianDiscrepancyDistanceT << " = {";
1187 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyDistanceT.size(); i != size; ++i)
1188 {
1189 os << this->m_gaussianDiscrepancyDistanceT[i];
1190 if (i+1 != size)
1191 os << ',';
1192 }
1193 os << '}';
1194
1195 os << "\n" << m_option_gaussianDiscrepancyPeriodicX << " = {";
1196 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyPeriodicX.size(); i != size; ++i)
1197 {
1198 os << this->m_gaussianDiscrepancyPeriodicX[i];
1199 if (i+1 != size)
1200 os << ',';
1201 }
1202 os << '}';
1203 os << "\n" << m_option_gaussianDiscrepancyPeriodicY << " = {";
1204 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyPeriodicY.size(); i != size; ++i)
1205 {
1206 os << this->m_gaussianDiscrepancyPeriodicY[i];
1207 if (i+1 != size)
1208 os << ',';
1209 }
1210 os << '}';
1211 os << "\n" << m_option_gaussianDiscrepancyPeriodicZ << " = {";
1212 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyPeriodicZ.size(); i != size; ++i)
1213 {
1214 os << this->m_gaussianDiscrepancyPeriodicZ[i];
1215 if (i+1 != size)
1216 os << ',';
1217 }
1218 os << '}';
1219 os << "\n" << m_option_gaussianDiscrepancyPeriodicT << " = {";
1220 for (unsigned int i = 0, size = this->m_gaussianDiscrepancyPeriodicT.size(); i != size; ++i)
1221 {
1222 os << this->m_gaussianDiscrepancyPeriodicT[i];
1223 if (i+1 != size)
1224 os << ',';
1225 }
1226 os << '}';
1227 os << "\n" << m_option_gaussianDiscrepancySupportThreshold << " = " << this->m_gaussianDiscrepancySupportThreshold
1228 << std::endl;
1229 }
1230
1231 std::ostream &
1232 operator<<(std::ostream& os, const GPMSAOptions & obj)
1233 {
1234 #ifndef QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
1235 os << (*(obj.m_parser)) << std::endl;
1236 #endif // QUESO_DISABLE_BOOST_PROGRAM_OPTIONS
1237 obj.print(os);
1238 return os;
1239 }
1240
1241
1242
1243 // Template instantiations
1244 template
1245 void
1246 GPMSAOptions::set_final_scaling<GslVector>
1247 (const std::vector<SharedPtr<GslVector>::Type> &,
1248 const std::vector<SharedPtr<GslVector>::Type> &,
1249 const std::vector<SharedPtr<GslVector>::Type> &,
1250 const std::vector<SharedPtr<GslVector>::Type> &,
1251 const std::vector<SharedPtr<GslVector>::Type> &,
1252 const std::vector<typename SharedPtr<SimulationOutputMesh<GslVector> >::Type> &);
1253
1254
1255 } // End namespace QUESO
1256