queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
GPMSA.C
Go to the documentation of this file.
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/GPMSA.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 template <class V, class M>
33  const VectorSet<V, M> & domain,
34  const VectorSpace<V, M> & m_scenarioSpace,
35  const VectorSpace<V, M> & m_parameterSpace,
36  const VectorSpace<V, M> & m_simulationOutputSpace,
37  const VectorSpace<V, M> & m_experimentOutputSpace,
38  const unsigned int m_numSimulations,
39  const unsigned int m_numExperiments,
40  const std::vector<typename SharedPtr<V>::Type> & m_simulationScenarios,
41  const std::vector<typename SharedPtr<V>::Type> & m_simulationParameters,
42  const std::vector<typename SharedPtr<V>::Type> & m_simulationOutputs,
43  const std::vector<typename SharedPtr<V>::Type> & m_experimentScenarios,
44  const std::vector<typename SharedPtr<V>::Type> & m_experimentOutputs,
45  const std::vector<typename SharedPtr<V>::Type> & m_discrepancyBases,
46  const std::vector<typename SharedPtr<M>::Type> & m_observationErrorMatrices,
47  const typename SharedPtr<M>::Type & m_observationErrorMatrix,
48  const ConcatenatedVectorRV<V, M> & m_totalPrior,
49  const V & residual_in,
50  const M & BT_Wy_B_inv_in,
51  const M & KT_K_inv_in,
52  const GPMSAOptions & opts)
53  :
54  BaseScalarFunction<V, M>("", m_totalPrior.imageSet()),
55  m_scenarioSpace(m_scenarioSpace),
56  m_parameterSpace(m_parameterSpace),
57  m_simulationOutputSpace(m_simulationOutputSpace),
58  m_experimentOutputSpace(m_experimentOutputSpace),
59  m_numSimulations(m_numSimulations),
60  m_numExperiments(m_numExperiments),
61  m_simulationScenarios(m_simulationScenarios),
62  m_simulationParameters(m_simulationParameters),
63  m_simulationOutputs(m_simulationOutputs),
64  m_experimentScenarios(m_experimentScenarios),
65  m_experimentOutputs(m_experimentOutputs),
66  m_discrepancyBases(m_discrepancyBases),
67  m_observationErrorMatrices(m_observationErrorMatrices),
68  m_observationErrorMatrix(m_observationErrorMatrix),
69  m_totalPrior(m_totalPrior),
70  residual(residual_in),
71  BT_Wy_B_inv(BT_Wy_B_inv_in),
72  KT_K_inv(KT_K_inv_in),
73  m_opts(opts)
74 {
75  queso_assert_greater(m_numSimulations, 0);
76 
77  queso_assert_equal_to
78  (m_simulationOutputs[0]->map().Comm().NumProc(), 1);
79 
80  const unsigned int numOutputs =
81  this->m_experimentOutputSpace.dimLocal();
82 
84  std::min((unsigned int)(this->m_opts.m_maxEmulatorBasisVectors), numOutputs) :
85  numOutputs;
86 }
87 
88 template <class V, class M>
90 {
91  // heap allocations get deleted via ScopedPtr
92 }
93 
94 template <class V, class M>
95 double
96 GPMSAEmulator<V, M>::lnValue(const V & domainVector,
97  const V * /* domainDirection */,
98  V * /* gradVector */,
99  M * /* hessianMatrix */,
100  V * /* hessianEffect */) const
101 {
102  // Components of domainVector:
103  // theta(1) // = "theta", "t" in Higdon et. al. 2008
104  // theta(2)
105  // ...
106  // theta(dimParameter)
107  // truncation_error_precision // = "lambda_eta", only exists in
108  // vector case
109  // emulator_precision(1) // = "lambda_{w i}" in vector case
110  // ... // = "lambda_eta" in scalar case,
111  // which is unrelated to the
112  // vector case lambda_eta
113  // emulator_precision(num_svd_terms)
114  // emulator_corr_strength(1) // = "rho_{eta k}"
115  // ... // dimScenario = "p_x", dimParameter = "p_t"
116  // emulator_corr_strength(dimScenario + dimParameter)
117  // discrepancy_precision(1) // = "lambda_delta" in scalar case,
118  // // "lambda_{v i}" in vector
119  // ...
120  // discrepancy_precision(F) // FIXME: F = 1, |G_1| = p_delta, for now.
121  // discrepancy_corr_strength(1) // = "rho_{delta k}" in scalar case,
122  // ... // "rho_{v i}" in vector
123  // discrepancy_corr_strength(dimScenario)
124  // emulator_data_precision(1) // = "small white noise", "small ridge"
125  // observation_error_precision // = "lambda_y", iff user-requested
126 
127  // Other variables:
128  // m_numSimulations // = "m"
129  // m_numExperiments // = "n"
130  // m_simulationScenarios // = "eta"
131  // m_experimentScenarios // = "y"
132  // m_experimentErrors // = "Sigma_y"; now obsoleted by "W_y"
133  // numOutputs // = "n_eta"
134  // // "n_y" := sum(n_y_i)
135  // // (== n*n_eta for us for now)
136  // dimScenario // = "p_x"
137  // num_svd_terms // = "p_eta"
138  // num_discrepancy_bases // = "p_delta"
139  // m_TruncatedSVD_simulationOutputs // = "K_eta"
140  // covMatrix // = "Sigma_D" in scalar case,
141  // "Sigma_zhat" in vector
142  // m_discrepancyMatrices // = "D_i"
143  // m_observationErrorMatrices // = "W_i"
144  // m_observationErrorMatrix // = "W_y"
145  // m_BMatrix // = "B"
146  // num_discrepancy_groups // = "F"
147  // m_emulatorPrecisionShapeVec // = "a_eta"
148  // 1.0/m_emulatorPrecisionScaleVec // = "b_eta"
149  // m_observationalPrecisionShapeVec // = "a_y"
150  // 1.0/m_observationalPrecisionScaleVec // = "b_y"
151 
152  // Construct covariance matrix
153  const unsigned int totalRuns = this->m_numExperiments + this->m_numSimulations;
154  const unsigned int numOutputs = this->m_experimentOutputSpace.dimLocal();
155  const unsigned int totalOutputs = totalRuns * numOutputs;
156  const unsigned int num_discrepancy_bases = m_discrepancyBases.size();
157  const unsigned int residualSize = (numOutputs == 1) ? totalOutputs :
158  totalRuns * num_svd_terms + m_numExperiments * num_discrepancy_bases;
159 
160  double prodScenario = 1.0;
161  double prodParameter = 1.0;
162  double prodDiscrepancy = 1.0;
163  unsigned int dimScenario = (this->m_scenarioSpace).dimLocal();
164  unsigned int dimParameter = (this->m_parameterSpace).dimLocal();
165 
166  // Length of prior+hyperprior inputs
167  unsigned int dimSum = 2 +
168  (this->num_svd_terms < numOutputs) +
169  (numOutputs > 1) +
170  m_opts.m_calibrateObservationalPrecision +
171  num_svd_terms +
172  dimParameter +
173  dimParameter +
174  dimScenario +
175  dimScenario; // yum
176 
177  // Offset for Sigma_eta equivalent in vector case
178  const unsigned int offset1 = (numOutputs == 1) ?
179  0 : m_numExperiments * num_discrepancy_bases;
180 
181  // Offset for Sigma_w in vector case
182  const unsigned int offset1b = offset1 +
183  m_numExperiments * num_svd_terms;
184 
185  // Offset for lambda_eta term in zhat covariance in vector case
186  const unsigned int offset2 = (numOutputs == 1) ?
187  0 : m_numExperiments * (num_discrepancy_bases + num_svd_terms);
188 
189  // This is cumbersome. All I want is a matrix.
190  const MpiComm & comm = domainVector.map().Comm();
191  Map z_map(residualSize, 0, comm);
192  M covMatrix(this->m_env, z_map, residualSize);
193 
194  typename SharedPtr<V>::Type domainVectorParameter
195  (new V(*(this->m_simulationParameters[0])));
196  for (unsigned int k = 0; k < dimParameter; k++) {
197  queso_assert (!queso_isnan(domainVector[k]));
198  (*domainVectorParameter)[k] = domainVector[k];
199  }
200 
201  // This for loop is a disaster and could do with a *lot* of optimisation
202  for (unsigned int i = 0; i < totalRuns; i++) {
203 
204  // Scenario and uncertain input variables, *not* normalized
205  typename SharedPtr<V>::Type scenario1;
206  typename SharedPtr<V>::Type parameter1;
207 
208  // Decide whether to do experiment part of the covariance matrix
209  // Get i-th simulation-parameter pair
210  if (i < this->m_numExperiments) {
211  // Experiment scenario (known)
212  scenario1 = (this->m_experimentScenarios)[i];
213 
214  // Experiment parameter (unknown)
215  parameter1 = domainVectorParameter;
216  }
217  else {
218  scenario1 =
219  (this->m_simulationScenarios)[i-this->m_numExperiments];
220  parameter1 =
221  (this->m_simulationParameters)[i-this->m_numExperiments];
222  }
223 
224  for (unsigned int j = 0; j < totalRuns; j++) {
225 
226  // Scenario and uncertain input variables, *not* normalized
227  typename SharedPtr<V>::Type scenario2;
228  typename SharedPtr<V>::Type parameter2;
229 
230  if (j < this->m_numExperiments) {
231  scenario2 = (this->m_experimentScenarios)[j];
232  parameter2 = domainVectorParameter;
233  }
234  else {
235  scenario2 =
236  (this->m_simulationScenarios)[j-this->m_numExperiments];
237  parameter2 =
238  (this->m_simulationParameters)[j-this->m_numExperiments];
239  }
240 
241  // Emulator component // = first term in (1)
242  prodScenario = 1.0;
243  unsigned int emulatorCorrStrStart =
244  dimParameter + (this->num_svd_terms < numOutputs) + num_svd_terms;
245  for (unsigned int k = 0; k < dimScenario; k++) {
246  const double & emulator_corr_strength =
247  domainVector[emulatorCorrStrStart+k];
248  double scenario_param1 =
249  m_opts.normalized_scenario_parameter(k, (*scenario1)[k]);
250  double scenario_param2 =
251  m_opts.normalized_scenario_parameter(k, (*scenario2)[k]);
252  prodScenario *= std::pow(emulator_corr_strength,
253  4.0 * (scenario_param1 - scenario_param2) *
254  (scenario_param1 - scenario_param2));
255  }
256 
257  queso_assert (!queso_isnan(prodScenario));
258 
259  // = second term in (1)
260  prodParameter = 1.0;
261  for (unsigned int k = 0; k < dimParameter; k++) {
262  queso_assert (!queso_isnan(domainVector[emulatorCorrStrStart+dimScenario+k]));
263  queso_assert (!queso_isnan((*parameter1)[k]));
264  queso_assert (!queso_isnan((*parameter2)[k]));
265  const double & emulator_corr_strength =
266  domainVector[emulatorCorrStrStart+dimScenario+k];
267  double uncertain_param1 =
268  m_opts.normalized_uncertain_parameter(k, (*parameter1)[k]);
269  double uncertain_param2 =
270  m_opts.normalized_uncertain_parameter(k, (*parameter2)[k]);
271  prodParameter *= std::pow(
272  emulator_corr_strength,
273  4.0 * (uncertain_param1 - uncertain_param2) *
274  (uncertain_param1 - uncertain_param2));
275  }
276 
277  queso_assert (!queso_isnan(prodParameter));
278 
279  // Sigma_eta in scalar case,
280  // [Sigma_u, Sigma_uw; Sigma_uw^T, Sigma_w] in vector case
281  for (unsigned int basis = 0; basis != num_svd_terms; ++basis)
282  {
283  // coefficient in (1)
284  // The relevant precision for Sigma_eta is lambda_eta; for
285  // Sigma_uw etc. it's lambda_wi and we skip lambda_eta
286  const double relevant_precision =
287  domainVector[dimParameter + basis +
288  (this->num_svd_terms<numOutputs) + (numOutputs>1)];
289  queso_assert_greater(relevant_precision, 0.0);
290 
291  const unsigned int stridei =
292  (i < this->m_numExperiments) ? m_numExperiments : m_numSimulations;
293  const unsigned int offseti =
294  (i < this->m_numExperiments) ? offset1 : offset1b - m_numExperiments;
295  const unsigned int stridej =
296  (j < this->m_numExperiments) ? m_numExperiments : m_numSimulations;
297  const unsigned int offsetj =
298  (j < this->m_numExperiments) ? offset1 : offset1b - m_numExperiments;
299 
300  covMatrix(offseti+basis*stridei+i,
301  offsetj+basis*stridej+j) =
302  prodScenario * prodParameter / relevant_precision;
303  }
304 
305  // If we're in the experiment cross correlation part, need extra
306  // foo: Sigma_delta/Sigma_v and Sigma_y
307  if (i < this->m_numExperiments && j < this->m_numExperiments) {
308  typename SharedPtr<V>::Type cross_scenario1 = (this->m_experimentScenarios)[i];
309  typename SharedPtr<V>::Type cross_scenario2 = (this->m_experimentScenarios)[j];
310  prodDiscrepancy = 1.0;
311  unsigned int discrepancyCorrStrStart =
312  dimParameter + num_svd_terms + dimParameter + dimScenario + 1 +
313  (this->num_svd_terms<numOutputs) + (numOutputs > 1);
314  for (unsigned int k = 0; k < dimScenario; k++) {
315  const double & discrepancy_corr_strength =
316  domainVector[discrepancyCorrStrStart+k];
317  double cross_scenario_param1 =
318  m_opts.normalized_scenario_parameter(k, (*cross_scenario1)[k]);
319  double cross_scenario_param2 =
320  m_opts.normalized_scenario_parameter(k, (*cross_scenario2)[k]);
321  prodDiscrepancy *=
322  std::pow(discrepancy_corr_strength, 4.0 *
323  (cross_scenario_param1 - cross_scenario_param2) *
324  (cross_scenario_param1 - cross_scenario_param2));
325  }
326 
327  const double discrepancy_precision =
328  domainVector[discrepancyCorrStrStart-1];
329  queso_assert_greater(discrepancy_precision, 0);
330  queso_assert (!queso_isnan(prodDiscrepancy));
331 
332  // Sigma_delta term from below (3) in univariate case
333  // Sigma_v term from p. 576 in multivariate case
334  const double R_v = prodDiscrepancy / discrepancy_precision;
335  for (unsigned int disc = 0; disc != num_discrepancy_bases;
336  ++disc)
337  covMatrix(disc*m_numExperiments+i,
338  disc*m_numExperiments+j) += R_v;
339 
340  if (numOutputs == 1)
341  {
342  // Experimental error comes in via W_y now.
343 /*
344  // Sigma_y term from below (3)
345  const double experimentalError =
346  this->m_experimentErrors(i,j);
347 */
348  const double experimentalError =
349  (*this->m_observationErrorMatrix)(i,j);
350 
351  queso_assert_greater_equal (experimentalError, 0);
352 
353  const double lambda_y =
354  m_opts.m_calibrateObservationalPrecision ?
355  domainVector[dimSum-1] : 1.0;
356 
357  covMatrix(i,j) += lambda_y * experimentalError;
358  }
359  }
360  }
361 
362  // Add small white noise component to diagonal to make stuff +ve def
363  // = "small ridge"
364  const double emulator_data_precision = domainVector[dimSum-1-(numOutputs>1)];
365  queso_assert_greater(emulator_data_precision, 0);
366  double nugget = 1.0 / emulator_data_precision;
367 
368  for (unsigned int disc = 0; disc != num_discrepancy_bases;
369  ++disc)
370  covMatrix(disc*m_numExperiments+i,
371  disc*m_numExperiments+i) += nugget;
372  }
373 
374  // If we're in the multivariate case, we've built the full Sigma_z
375  // matrix; now add the remaining Sigma_zhat terms
376  if (numOutputs > 1)
377  {
378  const double lambda_y =
379  m_opts.m_calibrateObservationalPrecision ?
380  domainVector[dimSum-1] : 1.0;
381  const double inv_lambda_y = 1.0/lambda_y;
382 
383  unsigned int BT_Wy_B_size = BT_Wy_B_inv.numCols();
384  for (unsigned int i=0; i != BT_Wy_B_size; ++i)
385  for (unsigned int j=0; j != BT_Wy_B_size; ++j)
386  covMatrix(i,j) += BT_Wy_B_inv(i,j) * inv_lambda_y;
387 
388  const double emulator_precision =
389  domainVector[dimParameter+1];
390  const double inv_emulator_precision = 1.0/emulator_precision;
391 
392  unsigned int KT_K_size = KT_K_inv.numCols();
393  for (unsigned int i=0; i != KT_K_size; ++i)
394  for (unsigned int j=0; j != KT_K_size; ++j)
395  covMatrix(i+offset2,j+offset2) +=
396  KT_K_inv(i,j) * inv_emulator_precision;
397  }
398 
399 
400  // Solve covMatrix * sol = residual
401  // = Sigma_D^-1 * (D - mu 1) from (3)
402  V sol(covMatrix.invertMultiply(residual));
403 
404  // Premultiply by residual^T as in (3)
405  double minus_2_log_lhd = 0.0;
406  // There's no dot product function in GslVector.
407  for (unsigned int i = 0; i < residualSize; i++) {
408  minus_2_log_lhd += sol[i] * residual[i];
409  }
410 
411 if (queso_isnan(minus_2_log_lhd))
412  for (unsigned int i = 0; i < residualSize; i++) {
413  if (queso_isnan(sol[i]))
414  std::cout << "NaN sol[" << i << ']' << std::endl;
415  if (queso_isnan(residual[i]))
416  std::cout << "NaN residual[" << i << ']' << std::endl;
417 
418  std::cout << "Covariance Matrix:" << std::endl;
419  covMatrix.print(std::cout);
420  }
421 
422 // std::cout << "minus_2_log_lhd = " << minus_2_log_lhd << std::endl;
423 
424  queso_assert_greater(minus_2_log_lhd, 0);
425 
426  double cov_det = covMatrix.determinant();
427 
428  if (cov_det <= 0)
429  {
430  std::cout << "Non-positive determinant for covMatrix = " << std::endl;
431  covMatrix.print(std::cout);
432  queso_error();
433  }
434 
435  minus_2_log_lhd += std::log(covMatrix.determinant());
436 
437  // Multiply by -1/2 coefficient from (3)
438  return -0.5 * minus_2_log_lhd;
439 }
440 
441 template <class V, class M>
442 double
443 GPMSAEmulator<V, M>::actualValue(const V & /* domainVector */,
444  const V * /* domainDirection */,
445  V * /* gradVector */,
446  M * /* hessianMatrix */,
447  V * /* hessianEffect */) const
448 {
449  // Do nothing?
450  return 1.0;
451 }
452 
453 template <class V, class M>
455  const BaseEnvironment & env,
456  GPMSAOptions * opts,
457  const BaseVectorRV<V, M> & parameterPrior,
458  const VectorSpace<V, M> & scenarioSpace,
459  const VectorSpace<V, M> & parameterSpace,
460  const VectorSpace<V, M> & simulationOutputSpace,
461  const VectorSpace<V, M> & experimentOutputSpace,
462  unsigned int numSimulations,
463  unsigned int numExperiments)
464  :
465  m_env(env),
466  m_parameterPrior(parameterPrior),
467  m_scenarioSpace(scenarioSpace),
468  m_parameterSpace(parameterSpace),
469  m_simulationOutputSpace(simulationOutputSpace),
470  m_experimentOutputSpace(experimentOutputSpace),
471  m_numSimulations(numSimulations),
472  m_numExperiments(numExperiments),
473  m_simulationScenarios(numSimulations),
474  m_simulationParameters(numSimulations),
475  m_simulationOutputs(numSimulations),
476  m_experimentScenarios(numExperiments),
477  m_experimentOutputs(numExperiments),
478  m_numSimulationAdds(0),
479  m_numExperimentAdds(0),
480  num_svd_terms(0),
481  priors(),
482  m_constructedGP(false)
483 {
484  // We should have the same number of outputs from both simulations
485  // and experiments
486  queso_assert_equal_to(simulationOutputSpace.dimGlobal(),
487  experimentOutputSpace.dimGlobal());
488 
489  {
490  const unsigned int numOutputs =
491  this->m_experimentOutputSpace.dimLocal();
492  const Map & output_map = experimentOutputSpace.map();
493 
494  // Set up the default discrepancy basis:
495  typename SharedPtr<V>::Type all_ones_basis(new V(env, output_map));
496  for (unsigned int i=0; i != numOutputs; ++i)
497  (*all_ones_basis)[i] = 1;
498  m_discrepancyBases.push_back(all_ones_basis);
499 
500  // Set up the default observation error covariance matrix:
501  for (unsigned int i = 0; i != numExperiments; ++i)
502  {
503  typename SharedPtr<M>::Type identity_matrix(new M(env, output_map, 1.0));
504  m_observationErrorMatrices.push_back(identity_matrix);
505  }
506  }
507 
508  // DM: Not sure if the logic in these 3 if-blocks is correct
509  if ((opts == NULL) && (this->m_env.optionsInputFileName() == "")) {
510  queso_error_msg("Must options object or an input file");
511  }
512 
513  if (opts != NULL) {
514  allocated_m_opts = false;
515  this->m_opts = opts;
516  }
517  else {
518  // Create a default one
519  allocated_m_opts = true;
520  this->m_opts = new GPMSAOptions(this->m_env, "");
521  }
522 
523  // FIXME: WTF? - RHS
524  // this->m_constructedGP = false;
525 }
526 
527 template <class V, class M>
529 {
530  if (this->allocated_m_opts)
531  delete this->m_opts;
532 }
533 
534 template <class V, class M>
535 const GPMSAOptions &
537 {
538  return *this->m_opts;
539 }
540 
541 
542 template <class V, class M>
543 GPMSAOptions &
545 {
546  return *this->m_opts;
547 }
548 
549 
550 template <class V, class M>
551 unsigned int
553 {
554  return this->m_numSimulations;
555 }
556 
557 template <class V, class M>
558 unsigned int
560 {
561  return this->m_numExperiments;
562 }
563 
564 template <class V, class M>
565 const VectorSpace<V, M> &
567 {
568  return this->m_scenarioSpace;
569 }
570 
571 template <class V, class M>
572 const VectorSpace<V, M> &
574 {
575  return this->m_parameterSpace;
576 }
577 
578 template <class V, class M>
579 const VectorSpace<V, M> &
581 {
582  return this->m_simulationOutputSpace;
583 }
584 
585 template <class V, class M>
586 const VectorSpace<V, M> &
588 {
589  return this->m_experimentOutputSpace;
590 }
591 
592 template <class V, class M>
593 const V &
595  unsigned int simulationId) const
596 {
597  queso_require_less_msg(simulationId, m_simulationScenarios.size(), "simulationId is too large");
598 
599  queso_require_msg(m_simulationScenarios[simulationId], "vector is NULL");
600 
601  return *(this->m_simulationScenarios[simulationId]);
602 }
603 
604 template <class V, class M>
605 const std::vector<typename SharedPtr<V>::Type> &
607 {
608  return this->m_simulationScenarios;
609 }
610 
611 template <class V, class M>
612 const V &
614  unsigned int simulationId) const
615 {
616  queso_require_less_msg(simulationId, m_simulationParameters.size(), "simulationId is too large");
617 
618  queso_require_msg(m_simulationParameters[simulationId], "vector is NULL");
619 
620  return *(this->m_simulationParameters[simulationId]);
621 }
622 
623 template <class V, class M>
624 const std::vector<typename SharedPtr<V>::Type> &
626 {
627  return this->m_simulationParameters;
628 }
629 
630 template <class V, class M>
631 const V &
633  unsigned int simulationId) const
634 {
635  queso_require_less_msg(simulationId, m_simulationOutputs.size(), "simulationId is too large");
636 
637  queso_require_msg(m_simulationOutputs[simulationId], "vector is NULL");
638 
639  return *(this->m_simulationOutputs[simulationId]);
640 }
641 
642 template <class V, class M>
643 const std::vector<typename SharedPtr<V>::Type> &
645 {
646  return this->m_simulationOutputs;
647 }
648 
649 template <class V, class M>
650 const V &
652  unsigned int experimentId) const
653 {
654  queso_require_less_msg(experimentId, (this->m_experimentScenarios).size(), "experimentId is too large");
655 
656  queso_require_msg(this->m_experimentScenarios[experimentId], "vector is NULL");
657 
658  return *(this->m_experimentScenarios[experimentId]);
659 }
660 
661 template <class V, class M>
662 const std::vector<typename SharedPtr<V>::Type> &
664 {
665  return this->m_experimentScenarios;
666 }
667 
668 template <class V, class M>
669 const V &
671  unsigned int experimentId) const
672 {
673  queso_require_less_msg(experimentId, (this->m_experimentOutputs).size(), "experimentId is too large");
674 
675  queso_require_msg(this->m_experimentOutputs[experimentId], "vector is NULL");
676 
677  return *(this->m_experimentOutputs[experimentId]);
678 }
679 
680 template <class V, class M>
681 const std::vector<typename SharedPtr<V>::Type> &
683 {
684  return this->m_experimentOutputs;
685 }
686 
687 template <class V, class M>
688 const BaseEnvironment &
690 {
691  return this->m_env;
692 }
693 
694 template <class V, class M>
695 const GPMSAEmulator<V, M> &
697 {
698  return *(this->gpmsaEmulator);
699 }
700 
701 template <class V, class M>
702 void
704  typename SharedPtr<V>::Type simulationParameter,
705  typename SharedPtr<V>::Type simulationOutput)
706 {
707  queso_require_less_msg(this->m_numSimulationAdds, this->m_numSimulations, "too many simulation adds...");
708 
709  this->m_simulationScenarios[this->m_numSimulationAdds] = simulationScenario;
710  this->m_simulationParameters[this->m_numSimulationAdds] = simulationParameter;
711  this->m_simulationOutputs[this->m_numSimulationAdds] = simulationOutput;
712  this->m_numSimulationAdds++;
713 
714  if ((this->m_numSimulationAdds == this->m_numSimulations) &&
715  (this->m_numExperimentAdds == this->m_numExperiments) &&
716  (this->m_constructedGP == false)) {
717  this->setUpEmulator();
718  }
719 }
720 
721 template <class V, class M>
722 void
724  const std::vector<typename SharedPtr<V>::Type> & simulationScenarios,
725  const std::vector<typename SharedPtr<V>::Type> & simulationParameters,
726  const std::vector<typename SharedPtr<V>::Type> & simulationOutputs)
727 {
728  for (unsigned int i = 0; i < this->m_numSimulations; i++) {
729  this->addSimulation(simulationScenarios[i],
730  simulationParameters[i],
731  simulationOutputs[i]);
732  }
733 }
734 
735 
736 
737 template <class V, class M>
738 void
740 {
741  this->m_opts->template set_final_scaling<V>
742  (m_simulationScenarios,
743  m_simulationParameters,
744  m_simulationOutputs,
745  m_experimentScenarios,
746  m_experimentOutputs);
747 
748  const unsigned int numOutputs =
749  this->m_experimentOutputSpace.dimLocal();
750 
751  this->num_svd_terms = this->m_opts->m_maxEmulatorBasisVectors ?
752  std::min((unsigned int)(this->m_opts->m_maxEmulatorBasisVectors), numOutputs) :
753  numOutputs;
754 
755  const Map & output_map = m_simulationOutputs[0]->map();
756 
757  const MpiComm & comm = output_map.Comm();
758 
759  Map serial_map(m_numSimulations, 0, comm);
760 
761  const BaseEnvironment &env = m_simulationOutputs[0]->env();
762 
763  simulationOutputMeans.reset
764  (new V (env, output_map));
765 
766  for (unsigned int i=0; i != m_numSimulations; ++i)
767  for (unsigned int j=0; j != numOutputs; ++j)
768  (*simulationOutputMeans)[j] += (*m_simulationOutputs[i])[j];
769 
770  for (unsigned int j=0; j != numOutputs; ++j)
771  (*simulationOutputMeans)[j] /= m_numSimulations;
772 
773  M simulation_matrix(env, serial_map, numOutputs);
774 
775  for (unsigned int i=0; i != m_numSimulations; ++i)
776  for (unsigned int j=0; j != numOutputs; ++j)
777  simulation_matrix(i,j) =
778  (*m_simulationOutputs[i])[j] - (*simulationOutputMeans)[j];
779 
780  // GSL only finds left singular vectors if n_rows>=n_columns, so we need to
781  // calculate them indirectly from the eigenvalues of M^T*M
782 
783  M S_trans(simulation_matrix.transpose());
784 
785  M SM_squared(S_trans*simulation_matrix);
786 
787  M SM_singularVectors(env, SM_squared.map(), numOutputs);
788  V SM_singularValues(env, SM_squared.map());
789 
790  SM_squared.eigen(SM_singularValues, &SM_singularVectors);
791 
792  // Copy only those vectors we want into K_eta
793  m_TruncatedSVD_simulationOutputs.reset
794  (new M(env, output_map, num_svd_terms));
795 
796  for (unsigned int i=0; i != numOutputs; ++i)
797  for (unsigned int k = 0; k != num_svd_terms; ++k)
798  (*m_TruncatedSVD_simulationOutputs)(i,k) = SM_singularVectors(i,k);
799 
800  Map copied_map(numOutputs * m_numSimulations, 0, comm);
801 
802  K.reset
803  (new M(env, copied_map, m_numSimulations * num_svd_terms));
804  for (unsigned int k=0; k != num_svd_terms; ++k)
805  for (unsigned int i1=0; i1 != m_numSimulations; ++i1)
806  for (unsigned int i2=0; i2 != numOutputs; ++i2)
807  {
808  const unsigned int i = i1 * numOutputs + i2;
809  const unsigned int j = k * m_numSimulations + i1;
810  (*K)(i,j) = SM_singularVectors(i2,k);
811  }
812 
813  KT_K_inv.reset
814  (new M((K->transpose() * *K).inverse()));
815 
816  Map serial_output_map(numOutputs, 0, comm);
817 
818  for (unsigned int i = 0; i != m_numExperiments; ++i)
819  {
820  M D_i(env, serial_output_map,
821  (unsigned int)(m_discrepancyBases.size()));
822 
823  for (unsigned int j=0; j != numOutputs; ++j)
824  for (unsigned int k=0; k != m_discrepancyBases.size(); ++k)
825  D_i(j,k) = (*m_discrepancyBases[k])[j];
826 
827  m_discrepancyMatrices.push_back(D_i);
828  }
829 
830  // Create the giant matrices!
831 
832  // We build B from two parts, the diag(D_i)*P_D^T block and
833  // the diag(K_i)*P_K^T block, but we build simultaneously so we only
834  // need one loop over experiments.
835  // Since P_D and P_K are permutation matrices we'll just apply the
836  // permutations as we insert.
837 
838  // W_y is simple block diagonal.
839 
840  // Neither of these ought to be dense matrices but we're sacrificing
841  // efficiency for clarity for now.
842 
843  const unsigned int num_discrepancy_bases = m_discrepancyBases.size();
844  const unsigned int Brows = m_numExperiments * numOutputs;
845  const unsigned int Bcols =
846  m_numExperiments * (num_discrepancy_bases + num_svd_terms);
847 
848  const Map B_row_map(Brows, 0, comm);
849 
850  m_BMatrix.reset
851  (new M(env, B_row_map, Bcols));
852 
853  const unsigned int Wyrows = m_numExperiments * numOutputs;
854 
855  const Map Wy_row_map(Wyrows, 0, comm);
856 
857  m_observationErrorMatrix.reset
858  (new M(env, Wy_row_map, Wyrows));
859 
860  M& B = *m_BMatrix;
861  M& Wy = *m_observationErrorMatrix;
862 
863  for (unsigned int ex = 0; ex != m_numExperiments; ++ex)
864  {
865  const M & D_i = m_discrepancyMatrices[ex];
866 
867  // Each error covariance matrix had better be SPD
868  const M & W_i = (*m_observationErrorMatrices[ex]);
869 #ifndef NDEBUG
870  M W_i_copy(W_i);
871  int rv = W_i_copy.chol();
872  queso_assert_msg(!rv, "Observation error matrix W_" << ex <<
873  " was not SPD!");
874 #endif
875 
876  // For the multivariate case, the bases K_eta computed from
877  // simulator outputs are the same as the bases K_i which apply
878  // to quantities of interest, because simulator outputs are QoIs
879  // alone.
880  //
881  // FIXME - we need to interpolate K_i in the functional case.
882 
883  for (unsigned int outi = 0; outi != numOutputs; ++outi)
884  {
885  unsigned int i = ex*numOutputs+outi;
886  for (unsigned int outj = 0; outj != num_discrepancy_bases; ++outj)
887  {
888  unsigned int j = ex + m_numExperiments * outj;
889 
890  B(i,j) = D_i(outi,outj);
891  }
892 
893  for (unsigned int outj = 0; outj != num_svd_terms; ++outj)
894  {
895  unsigned int j = ex +
896  m_numExperiments * (num_discrepancy_bases + outj);
897 
898  B(i,j) = (*m_TruncatedSVD_simulationOutputs)(outi,outj);
899  }
900 
901  for (unsigned int outj = 0; outj != numOutputs; ++outj)
902  {
903  // No fancy perturbation here
904  unsigned int j = ex*numOutputs+outj;
905 
906  Wy(i,j) = W_i(outi,outj);
907  }
908  }
909  }
910 
911  M BT_Wy_B (B.transpose() * Wy * B);
912 
913  // Adding a "small ridge" to make sure this is invertible, as on
914  // p.577 - using 1e-4 from discussion notes.
915  for (unsigned int i=0; i != Brows; ++i)
916  BT_Wy_B(i,i) += 1.e-4;
917 
918  BT_Wy_B_inv.reset(new M(BT_Wy_B.inverse()));
919 
920  this->setUpHyperpriors();
921 
922  this->m_constructedGP = true;
923  this->gpmsaEmulator.reset
924  (new GPMSAEmulator<V, M>(
925  this->prior().imageSet(),
926  this->m_scenarioSpace,
927  this->m_parameterSpace,
928  this->m_simulationOutputSpace,
929  this->m_experimentOutputSpace,
930  this->m_numSimulations,
931  this->m_numExperiments,
932  this->m_simulationScenarios,
933  this->m_simulationParameters,
934  this->m_simulationOutputs,
935  this->m_experimentScenarios,
936  this->m_experimentOutputs,
937  this->m_discrepancyBases,
938  this->m_observationErrorMatrices,
939  this->m_observationErrorMatrix,
940  *(this->m_totalPrior),
941  *this->residual,
942  *this->BT_Wy_B_inv,
943  *this->KT_K_inv,
944  *this->m_opts));
945 }
946 
947 
948 
949 template <class V, class M>
950 void
952  const std::vector<typename SharedPtr<V>::Type> & experimentScenarios,
953  const std::vector<typename SharedPtr<V>::Type> & experimentOutputs,
954  const typename SharedPtr<M>::Type experimentErrors)
955 {
956  queso_require_less_equal_msg(experimentScenarios.size(), this->m_numExperiments, "too many experiments...");
957 
958  unsigned int offset = 0;
959  for (unsigned int i = 0; i < this->m_experimentScenarios.size(); i++) {
960  this->m_experimentScenarios[i] = experimentScenarios[i];
961  this->m_experimentOutputs[i] = experimentOutputs[i];
962 
963  const unsigned int outsize =
964  this->m_experimentOutputs[i]->sizeGlobal();
965  for (unsigned int outi = 0; outi != outsize; ++outi)
966  for (unsigned int outj = 0; outj != outsize; ++outj)
967  (*this->m_observationErrorMatrices[i])(outi,outj) =
968  (*experimentErrors)(offset+outi, offset+outj);
969 
970  offset += outsize;
971  }
972 
973  this->m_numExperimentAdds += experimentScenarios.size();
974 
975  if ((this->m_numSimulationAdds == this->m_numSimulations) &&
976  (this->m_numExperimentAdds == this->m_numExperiments) &&
977  (this->m_constructedGP == false)) {
978  this->setUpEmulator();
979  }
980 }
981 
982 
983 template <class V, class M>
984 void
986  const std::vector<typename SharedPtr<V>::Type> & experimentScenarios,
987  const std::vector<typename SharedPtr<V>::Type> & experimentOutputs,
988  const std::vector<typename SharedPtr<M>::Type> & experimentErrors)
989 {
990  queso_require_less_equal_msg(experimentScenarios.size(), this->m_numExperiments, "too many experiments...");
991  queso_require_equal_to(experimentScenarios.size(),
992  experimentOutputs.size());
993  queso_require_equal_to(experimentScenarios.size(),
994  experimentErrors.size());
995 
996  for (unsigned int i = 0; i < this->m_experimentScenarios.size(); i++) {
997  this->m_experimentScenarios[i] = experimentScenarios[i];
998  this->m_experimentOutputs[i] = experimentOutputs[i];
999  this->m_observationErrorMatrices[i] = experimentErrors[i];
1000  }
1001  this->m_numExperimentAdds += experimentScenarios.size();
1002 
1003  if ((this->m_numSimulationAdds == this->m_numSimulations) &&
1004  (this->m_numExperimentAdds == this->m_numExperiments) &&
1005  (this->m_constructedGP == false)) {
1006  this->setUpEmulator();
1007  }
1008 }
1009 
1010 
1011 template <class V, class M>
1012 void
1014  const std::vector<typename SharedPtr<V>::Type> & discrepancyBases)
1015 {
1016  m_discrepancyBases = discrepancyBases;
1017 
1018  // We should not yet have constructed the underlying GP model
1019  queso_assert_equal_to(this->m_constructedGP, false);
1020 }
1021 
1022 
1023 
1024 template <class V, class M>
1025 M &
1027  (unsigned int simulationNumber)
1028 {
1029  queso_assert_less(simulationNumber, m_numSimulations);
1030  queso_assert_equal_to(m_observationErrorMatrices.size(), m_numSimulations);
1031 
1032  return *m_observationErrorMatrices[simulationNumber];
1033 }
1034 
1035 
1036 template <class V, class M>
1037 const M &
1039  (unsigned int simulationNumber) const
1040 {
1041  queso_assert_less(simulationNumber, m_numSimulations);
1042  queso_assert_equal_to(m_observationErrorMatrices.size(), m_numSimulations);
1043 
1044  return *m_observationErrorMatrices[simulationNumber];
1045 }
1046 
1047 
1048 
1049 template <class V, class M>
1052 {
1053  return *(this->m_totalPrior);
1054 }
1055 
1056 template <class V, class M>
1057 void
1058 GPMSAFactory<V, M>::print(std::ostream& /* os */) const
1059 {
1060  // Do nothing
1061 }
1062 
1063 // Private methods follow
1064 template <class V, class M>
1065 void
1067 {
1068  const unsigned int numOutputs =
1069  this->m_experimentOutputSpace.dimLocal();
1070 
1071  const unsigned int num_discrepancy_bases = m_discrepancyBases.size();
1072 
1073  const MpiComm & comm = m_simulationOutputs[0]->map().Comm();
1074 
1075  unsigned int rank_B;
1076  if (m_BMatrix->numRowsGlobal() > m_BMatrix->numCols())
1077  rank_B = m_BMatrix->rank(0, 1.e-4);
1078  else
1079  rank_B = m_BMatrix->transpose().rank(0, 1.e-4);
1080 
1081  double emulatorPrecisionShape = this->m_opts->m_emulatorPrecisionShape;
1082  double emulatorPrecisionScale = this->m_opts->m_emulatorPrecisionScale;
1083 
1084  double observationalPrecisionShape = this->m_opts->m_observationalPrecisionShape;
1085  double observationalPrecisionScale = this->m_opts->m_observationalPrecisionScale;
1086 
1087  if (numOutputs > 1)
1088  {
1089  Map y_map(m_numExperiments * numOutputs, 0, comm);
1090  Map eta_map(m_numSimulations * numOutputs, 0, comm);
1091 
1092  const unsigned int yhat_size =
1093  m_numExperiments * (num_discrepancy_bases + num_svd_terms);
1094 
1095  const unsigned int etahat_size =
1096  m_numSimulations * num_svd_terms;
1097 
1098  Map zhat_map(yhat_size + etahat_size, 0, comm);
1099 
1100  V y(this->m_env, y_map);
1101  V eta(this->m_env, eta_map);
1102 
1103  for (unsigned int i = 0; i < this->m_numExperiments; i++) {
1104  for (unsigned int k = 0; k != numOutputs; ++k)
1105  y[i*numOutputs+k] =
1106  (*((this->m_experimentOutputs)[i]))[k] -
1107  (*simulationOutputMeans)[k];
1108  }
1109 
1110  for (unsigned int i = 0; i < this->m_numSimulations; i++) {
1111  for (unsigned int k = 0; k != numOutputs; ++k)
1112  eta[i*numOutputs+k] =
1113  (*((this->m_simulationOutputs)[i]))[k] -
1114  (*simulationOutputMeans)[k];
1115  }
1116 
1117  M& B = *m_BMatrix;
1118  M& Wy = *m_observationErrorMatrix;
1119 
1120  V yhat(*BT_Wy_B_inv * (B.transpose() * (Wy * y)));
1121 
1122  queso_assert_equal_to(yhat.sizeGlobal(), yhat_size);
1123 
1124  V etahat(*KT_K_inv * (K->transpose() * eta));
1125 
1126  residual.reset(new V(this->m_env, zhat_map));
1127  for (unsigned int i = 0; i < yhat_size; ++i)
1128  (*residual)[i] = yhat[i];
1129 
1130  for (unsigned int i = 0; i < etahat_size; ++i)
1131  (*residual)[yhat_size+i] = etahat[i];
1132 
1133  emulatorPrecisionShape +=
1134  (this->m_numSimulations * (numOutputs - num_svd_terms)) / 2.0;
1135 
1136  V eta_temp(eta);
1137  eta_temp -= *K * etahat;
1138 
1139  emulatorPrecisionScale +=
1140  scalarProduct(eta, eta_temp) / 2.0;
1141 
1142  observationalPrecisionShape +=
1143  (this->m_numExperiments * numOutputs - rank_B) / 2.0;
1144 
1145  V y_temp(Wy * y);
1146  y_temp -= Wy * B * yhat;
1147 
1148  observationalPrecisionScale +=
1149  scalarProduct(y, y_temp) / 2.0;
1150  }
1151  else
1152  {
1153  const unsigned int totalRuns = this->m_numExperiments + this->m_numSimulations;
1154  Map z_map(totalRuns, 0, comm);
1155  residual.reset(new V (this->m_env, z_map));
1156 
1157  // Form residual = D - mean // = D - mu*1 in (3)
1158  // We currently use the mean of the simulation data, not a free
1159  // hyperparameter mean
1160  for (unsigned int i = 0; i < this->m_numExperiments; i++) {
1161  (*residual)[i] = (*((this->m_experimentOutputs)[i]))[0] -
1162  (*simulationOutputMeans)[0];
1163  }
1164  for (unsigned int i = 0; i < this->m_numSimulations; i++) {
1165  (*residual)[i+this->m_numExperiments] =
1166  (*((this->m_simulationOutputs)[i]))[0] -
1167  (*simulationOutputMeans)[0];
1168  }
1169  }
1170 
1171 
1172  double emulatorCorrelationStrengthAlpha = this->m_opts->m_emulatorCorrelationStrengthAlpha;
1173  double emulatorCorrelationStrengthBeta = this->m_opts->m_emulatorCorrelationStrengthBeta;
1174  double discrepancyPrecisionShape = this->m_opts->m_discrepancyPrecisionShape;
1175  double discrepancyPrecisionScale = this->m_opts->m_discrepancyPrecisionScale;
1176  double discrepancyCorrelationStrengthAlpha = this->m_opts->m_discrepancyCorrelationStrengthAlpha;
1177  double discrepancyCorrelationStrengthBeta = this->m_opts->m_discrepancyCorrelationStrengthBeta;
1178  double emulatorDataPrecisionShape = this->m_opts->m_emulatorDataPrecisionShape;
1179  double emulatorDataPrecisionScale = this->m_opts->m_emulatorDataPrecisionScale;
1180 
1181  this->oneDSpace.reset
1182  (new VectorSpace<V, M>(this->m_env, "", 1, NULL));
1183 
1184  // Truncation error precision
1185  if (this->num_svd_terms < numOutputs)
1186  {
1187  this->truncationErrorPrecisionMin.reset(new V(this->oneDSpace->zeroVector()));
1188  this->truncationErrorPrecisionMax.reset(new V(this->oneDSpace->zeroVector()));
1189  this->truncationErrorPrecisionMin->cwSet(-INFINITY);
1190  this->truncationErrorPrecisionMax->cwSet(INFINITY);
1191 
1192  this->truncationErrorPrecisionDomain.reset
1193  (new BoxSubset<V, M>
1194  ("",
1195  *(this->oneDSpace),
1196  *(this->truncationErrorPrecisionMin),
1197  *(this->truncationErrorPrecisionMax)));
1198 
1199  this->m_truncationErrorPrecision.reset
1201  ("",
1202  *(this->truncationErrorPrecisionDomain)));
1203  }
1204 
1205  // Emulator precision
1206  this->emulatorPrecisionSpace.reset
1207  (new VectorSpace<V, M>
1208  (this->m_env,
1209  "",
1210  num_svd_terms + (numOutputs > 1),
1211  NULL));
1212 
1213  this->emulatorPrecisionMin.reset
1214  (new V(this->emulatorPrecisionSpace->zeroVector()));
1215  this->emulatorPrecisionMax.reset
1216  (new V(this->emulatorPrecisionSpace->zeroVector()));
1217  this->m_emulatorPrecisionShapeVec.reset
1218  (new V(this->emulatorPrecisionSpace->zeroVector()));
1219  this->m_emulatorPrecisionScaleVec.reset
1220  (new V(this->emulatorPrecisionSpace->zeroVector()));
1221  this->emulatorPrecisionMin->cwSet(0.3);
1222  this->emulatorPrecisionMax->cwSet(INFINITY);
1223  this->m_emulatorPrecisionShapeVec->cwSet(emulatorPrecisionShape);
1224  this->m_emulatorPrecisionScaleVec->cwSet(emulatorPrecisionScale);
1225 
1226  this->emulatorPrecisionDomain.reset
1227  (new BoxSubset<V, M>
1228  ("",
1229  *(this->emulatorPrecisionSpace),
1230  *(this->emulatorPrecisionMin),
1231  *(this->emulatorPrecisionMax)));
1232 
1233  this->m_emulatorPrecision.reset
1234  (new GammaVectorRV<V, M>
1235  ("",
1236  *(this->emulatorPrecisionDomain),
1237  *(this->m_emulatorPrecisionShapeVec),
1238  *(this->m_emulatorPrecisionScaleVec)));
1239 
1240  // Emulator correlation strength
1241  unsigned int dimScenario = (this->scenarioSpace()).dimLocal();
1242  unsigned int dimParameter = (this->parameterSpace()).dimLocal();
1243  this->emulatorCorrelationSpace.reset
1244  (new VectorSpace<V, M>
1245  (this->m_env,
1246  "",
1247  dimScenario + dimParameter,
1248  NULL));
1249 
1250  this->emulatorCorrelationMin.reset
1251  (new V(this->emulatorCorrelationSpace->zeroVector()));
1252  this->emulatorCorrelationMax.reset
1253  (new V(this->emulatorCorrelationSpace->zeroVector()));
1254  this->m_emulatorCorrelationStrengthAlphaVec.reset
1255  (new V(this->emulatorCorrelationSpace->zeroVector()));
1256  this->m_emulatorCorrelationStrengthBetaVec.reset
1257  (new V(this->emulatorCorrelationSpace->zeroVector()));
1258  this->emulatorCorrelationMin->cwSet(0);
1259  this->emulatorCorrelationMax->cwSet(1);
1260  this->m_emulatorCorrelationStrengthAlphaVec->cwSet(emulatorCorrelationStrengthAlpha);
1261  this->m_emulatorCorrelationStrengthBetaVec->cwSet(emulatorCorrelationStrengthBeta);
1262 
1263  this->emulatorCorrelationDomain.reset
1264  (new BoxSubset<V, M>
1265  ("",
1266  *(this->emulatorCorrelationSpace),
1267  *(this->emulatorCorrelationMin),
1268  *(this->emulatorCorrelationMax)));
1269 
1270  this->m_emulatorCorrelationStrength.reset
1271  (new BetaVectorRV<V, M>
1272  ("",
1273  *(this->emulatorCorrelationDomain),
1274  *(this->m_emulatorCorrelationStrengthAlphaVec),
1275  *(this->m_emulatorCorrelationStrengthBetaVec)));
1276 
1277  // Observation precision
1278  this->observationalPrecisionSpace.reset
1279  (new VectorSpace<V, M>
1280  (this->m_env,
1281  "",
1282  1,
1283  NULL));
1284 
1285  this->observationalPrecisionMin.reset
1286  (new V(this->observationalPrecisionSpace->zeroVector()));
1287  this->observationalPrecisionMax.reset
1288  (new V(this->observationalPrecisionSpace->zeroVector()));
1289  this->m_observationalPrecisionShapeVec.reset
1290  (new V(this->observationalPrecisionSpace->zeroVector()));
1291  this->m_observationalPrecisionScaleVec.reset
1292  (new V(this->observationalPrecisionSpace->zeroVector()));
1293  this->observationalPrecisionMin->cwSet(0.3);
1294  this->observationalPrecisionMax->cwSet(INFINITY);
1295  this->m_observationalPrecisionShapeVec->cwSet
1296  (this->m_opts->m_observationalPrecisionShape);
1297  this->m_observationalPrecisionScaleVec->cwSet
1298  (this->m_opts->m_observationalPrecisionScale);
1299 
1300  this->observationalPrecisionDomain.reset
1301  (new BoxSubset<V, M>
1302  ("",
1303  *(this->observationalPrecisionSpace),
1304  *(this->observationalPrecisionMin),
1305  *(this->observationalPrecisionMax)));
1306 
1307  this->m_observationalPrecision.reset
1308  (new GammaVectorRV<V, M>
1309  ("",
1310  *(this->observationalPrecisionDomain),
1311  *(this->m_observationalPrecisionShapeVec),
1312  *(this->m_observationalPrecisionScaleVec)));
1313 
1314  // Discrepancy precision
1315  this->discrepancyPrecisionMin.reset
1316  (new V(this->oneDSpace->zeroVector()));
1317  this->discrepancyPrecisionMax.reset
1318  (new V(this->oneDSpace->zeroVector()));
1319  this->m_discrepancyPrecisionShapeVec.reset
1320  (new V(this->oneDSpace->zeroVector()));
1321  this->m_discrepancyPrecisionScaleVec.reset
1322  (new V(this->oneDSpace->zeroVector()));
1323  this->discrepancyPrecisionMin->cwSet(0);
1324  this->discrepancyPrecisionMax->cwSet(INFINITY);
1325  this->m_discrepancyPrecisionShapeVec->cwSet(discrepancyPrecisionShape);
1326  this->m_discrepancyPrecisionScaleVec->cwSet(discrepancyPrecisionScale);
1327 
1328  this->discrepancyPrecisionDomain.reset
1329  (new BoxSubset<V, M>
1330  ("",
1331  *(this->oneDSpace),
1332  *(this->discrepancyPrecisionMin),
1333  *(this->discrepancyPrecisionMax)));
1334 
1335  this->m_discrepancyPrecision.reset
1336  (new GammaVectorRV<V, M>
1337  ("",
1338  *(this->discrepancyPrecisionDomain),
1339  *(this->m_discrepancyPrecisionShapeVec),
1340  *(this->m_discrepancyPrecisionScaleVec)));
1341 
1342  // Discrepancy correlation strength
1343  this->discrepancyCorrelationSpace.reset
1344  (new VectorSpace<V, M>
1345  (this->m_env,
1346  "",
1347  dimScenario,
1348  NULL));
1349 
1350  this->discrepancyCorrelationMin.reset
1351  (new V(this->discrepancyCorrelationSpace->zeroVector()));
1352  this->discrepancyCorrelationMax.reset
1353  (new V(this->discrepancyCorrelationSpace->zeroVector()));
1354  this->m_discrepancyCorrelationStrengthAlphaVec.reset
1355  (new V(this->discrepancyCorrelationSpace->zeroVector()));
1356  this->m_discrepancyCorrelationStrengthBetaVec.reset
1357  (new V(this->discrepancyCorrelationSpace->zeroVector()));
1358  this->discrepancyCorrelationMin->cwSet(0);
1359  this->discrepancyCorrelationMax->cwSet(1);
1360  this->m_discrepancyCorrelationStrengthAlphaVec->cwSet(discrepancyCorrelationStrengthAlpha);
1361  this->m_discrepancyCorrelationStrengthBetaVec->cwSet(discrepancyCorrelationStrengthBeta);
1362 
1363  this->discrepancyCorrelationDomain.reset
1364  (new BoxSubset<V, M>
1365  ("",
1366  *(this->discrepancyCorrelationSpace),
1367  *(this->discrepancyCorrelationMin),
1368  *(this->discrepancyCorrelationMax)));
1369 
1370  this->m_discrepancyCorrelationStrength.reset
1371  (new BetaVectorRV<V, M>
1372  ("",
1373  *(this->discrepancyCorrelationDomain),
1374  *(this->m_discrepancyCorrelationStrengthAlphaVec),
1375  *(this->m_discrepancyCorrelationStrengthBetaVec)));
1376 
1377  // Emulator data precision
1378  this->emulatorDataPrecisionMin.reset
1379  (new V(this->oneDSpace->zeroVector()));
1380  this->emulatorDataPrecisionMax.reset
1381  (new V(this->oneDSpace->zeroVector()));
1382  this->m_emulatorDataPrecisionShapeVec.reset
1383  (new V(this->oneDSpace->zeroVector()));
1384  this->m_emulatorDataPrecisionScaleVec.reset
1385  (new V(this->oneDSpace->zeroVector()));
1386  this->emulatorDataPrecisionMin->cwSet(60.0);
1387  this->emulatorDataPrecisionMax->cwSet(1e5);
1388  this->m_emulatorDataPrecisionShapeVec->cwSet(emulatorDataPrecisionShape);
1389  this->m_emulatorDataPrecisionScaleVec->cwSet(emulatorDataPrecisionScale);
1390 
1391  this->emulatorDataPrecisionDomain.reset
1392  (new BoxSubset<V, M>
1393  ("",
1394  *(this->oneDSpace),
1395  *(this->emulatorDataPrecisionMin),
1396  *(this->emulatorDataPrecisionMax)));
1397 
1398  this->m_emulatorDataPrecision.reset
1399  (new GammaVectorRV<V, M>
1400  ("",
1401  *(this->emulatorDataPrecisionDomain),
1402  *(this->m_emulatorDataPrecisionShapeVec),
1403  *(this->m_emulatorDataPrecisionScaleVec)));
1404 
1405  // Now form full prior
1406  unsigned int dimSum = 2 +
1407  (this->num_svd_terms < numOutputs) +
1408  (numOutputs > 1) +
1409  this->m_opts->m_calibrateObservationalPrecision +
1410  num_svd_terms +
1411  dimParameter +
1412  dimParameter +
1413  dimScenario +
1414  dimScenario; // yum
1415 
1416  this->totalSpace.reset
1417  (new VectorSpace<V, M>
1418  (this->m_env,
1419  "",
1420  dimSum,
1421  NULL));
1422  this->totalMins.reset(new V(this->totalSpace->zeroVector()));
1423  this->totalMaxs.reset(new V(this->totalSpace->zeroVector()));
1424 
1425  // Hackety hack McHackington. There's no better way to do this unfortunately
1426  this->totalMins->cwSet(0);
1427  this->totalMaxs->cwSet(1);
1428 
1429  // Min emulator precision
1430  (*(this->totalMins))[dimParameter] = 0.3;
1431  // Max emulator precision
1432  (*(this->totalMaxs))[dimParameter] = INFINITY;
1433 
1434  if (numOutputs > 1)
1435  for (unsigned int basis = 0; basis != num_svd_terms; ++basis)
1436  {
1437  // Min weights precision
1438  (*(this->totalMins))[dimParameter+1+basis] = 0.3;
1439  // Max weights precision
1440  (*(this->totalMaxs))[dimParameter+1+basis] = INFINITY;
1441  }
1442 
1443  // FIXME: F = 1 for now
1444  // Min discrepancy precision
1445  (*(this->totalMins))[dimParameter+(numOutputs>1)+num_svd_terms+dimScenario+dimParameter] = 0;
1446  // Max discrepancy precision
1447  (*(this->totalMaxs))[dimParameter+(numOutputs>1)+num_svd_terms+dimScenario+dimParameter] = INFINITY;
1448 
1449  const int emulator_data_precision_index =
1450  dimSum - 1 - this->m_opts->m_calibrateObservationalPrecision;
1451  (*(this->totalMins))[emulator_data_precision_index] = 60.0; // Min emulator data precision
1452  (*(this->totalMaxs))[emulator_data_precision_index] = 1e5; // Max emulator data precision
1453 
1454  if (this->m_opts->m_calibrateObservationalPrecision) {
1455  (*(this->totalMins))[dimSum-1] = 0.3; // Min observation error precision
1456  (*(this->totalMaxs))[dimSum-1] = INFINITY; // Max observation error precision
1457  }
1458 
1459  this->totalDomain.reset
1460  (new BoxSubset<V, M>
1461  ("",
1462  *(this->totalSpace),
1463  *(this->totalMins),
1464  *(this->totalMaxs)));
1465 
1466  this->priors.push_back(&(this->m_parameterPrior));
1467 
1468  if (this->num_svd_terms < numOutputs)
1469  this->priors.push_back(this->m_truncationErrorPrecision.get());
1470 
1471  this->priors.push_back(this->m_emulatorPrecision.get());
1472  this->priors.push_back(this->m_emulatorCorrelationStrength.get());
1473  this->priors.push_back(this->m_discrepancyPrecision.get());
1474  this->priors.push_back(this->m_discrepancyCorrelationStrength.get());
1475  this->priors.push_back(this->m_emulatorDataPrecision.get());
1476  if (this->m_opts->m_calibrateObservationalPrecision)
1477  this->priors.push_back(this->m_observationalPrecision.get());
1478 
1479  // Finally
1480  this->m_totalPrior.reset
1482  ("",
1483  this->priors,
1484  *(this->totalDomain)));
1485 }
1486 
1487 } // End namespace QUESO
1488 
std::shared_ptr< T > Type
Definition: SharedPtr.h:69
void setUpHyperpriors()
Definition: GPMSA.C:1066
unsigned int dimGlobal() const
Definition: VectorSpace.C:161
A templated (base) class for handling scalar functions.
std::vector< typename SharedPtr< M >::Type > m_observationErrorMatrices
Definition: GPMSA.h:357
const GPMSAOptions & options() const
Return GPMSAOptions structure.
Definition: GPMSA.C:536
A class for partitioning vectors and matrices.
Definition: Map.h:49
std::vector< typename SharedPtr< V >::Type > m_discrepancyBases
Definition: GPMSA.h:355
unsigned int numSimulations() const
Return number of simulations.
Definition: GPMSA.C:552
const V & simulationParameter(unsigned int simulationId) const
Return the point in parameterSpace for simulation simulationId.
Definition: GPMSA.C:613
A class representing a vector space.
Definition: VectorSet.h:49
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
const BaseEnvironment & m_env
Definition: GPMSA.h:334
M & getObservationErrorCovariance(unsigned int simulationNumber)
Definition: GPMSA.C:1027
const BaseEnvironment & env() const
Return the QUESO environment.
Definition: GPMSA.C:689
const VectorSpace< V, M > & parameterSpace() const
Return the vector space in which parameters live.
Definition: GPMSA.C:573
A class representing a uniform vector RV.
const VectorSpace< V, M > & simulationOutputSpace() const
Return the vector space in which simulations live.
Definition: GPMSA.C:580
const GPMSAEmulator< V, M > & getGPMSAEmulator() const
Return the GPMSAEmulator likelihood object.
Definition: GPMSA.C:696
const VectorSpace< V, M > & experimentOutputSpace() const
Return the vector space in which experiments live.
Definition: GPMSA.C:587
const std::vector< typename SharedPtr< V >::Type > & experimentOutputs() const
Return all points in experimentOutputSpace for all experiments.
Definition: GPMSA.C:682
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:354
A class representing a vector RV constructed via Beta distribution.
Definition: BetaVectorRV.h:62
GPMSAOptions * m_opts
Definition: GPMSA.h:482
virtual double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the scalar function. Deprecated.
Definition: GPMSA.C:96
int m_maxEmulatorBasisVectors
The maximum number of basis vectors to use for approximating.
Definition: GPMSAOptions.h:80
bool queso_isnan(T arg)
Definition: math_macros.h:39
A templated class for handling sets.
Definition: VectorSet.h:52
unsigned int numExperiments() const
Return number of experiments.
Definition: GPMSA.C:559
const std::vector< typename SharedPtr< V >::Type > & simulationOutputs() const
Return all points in simulationOutputSpace for all simulations.
Definition: GPMSA.C:644
double scalarProduct(const GslVector &x, const GslVector &y)
Definition: GslVector.C:1137
unsigned int num_svd_terms
Definition: GPMSA.h:110
A class representing concatenated vector RVs.
const V & experimentOutput(unsigned int experimentId) const
Return the experiment output for experiment experimentId.
Definition: GPMSA.C:670
void addSimulation(typename SharedPtr< V >::Type simulationScenario, typename SharedPtr< V >::Type simulationParameter, typename SharedPtr< V >::Type simulationOutput)
Add a simulation to this.
Definition: GPMSA.C:703
~GPMSAFactory()
Destructor.
Definition: GPMSA.C:528
Class representing a subset of a vector space shaped like a hypercube.
Definition: BoxSubset.h:44
const GPMSAOptions & m_opts
Definition: GPMSA.h:123
void print(std::ostream &os) const
Definition: GPMSA.C:1058
void addExperiments(const std::vector< typename SharedPtr< V >::Type > &experimentScenarios, const std::vector< typename SharedPtr< V >::Type > &experimentOutputs, const typename SharedPtr< M >::Type experimentErrors)
Add all experiments to this.
Definition: GPMSA.C:951
GPMSAEmulator(const VectorSet< V, M > &domain, const VectorSpace< V, M > &m_scenarioSpace, const VectorSpace< V, M > &m_parameterSpace, const VectorSpace< V, M > &m_simulationOutputSpace, const VectorSpace< V, M > &m_experimentOutputSpace, const unsigned int m_numSimulations, const unsigned int m_numExperiments, 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< V >::Type > &m_discrepancyBases, const std::vector< typename SharedPtr< M >::Type > &m_observationErrorMatrices, const typename SharedPtr< M >::Type &m_observationErrorMatrix, const ConcatenatedVectorRV< V, M > &m_totalPrior, const V &residual_in, const M &BT_Wy_B_inv_in, const M &KT_K_inv_in, const GPMSAOptions &opts)
Definition: GPMSA.C:32
virtual double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the scalar function.
Definition: GPMSA.C:443
GPMSAFactory(const BaseEnvironment &env, GPMSAOptions *opts, const BaseVectorRV< V, M > &parameterPrior, const VectorSpace< V, M > &scenarioSpace, const VectorSpace< V, M > &parameterSpace, const VectorSpace< V, M > &simulationOutputSpace, const VectorSpace< V, M > &experimentOutputSpace, unsigned int numSimulations, unsigned int numExperiments)
Constructor.
Definition: GPMSA.C:454
const ConcatenatedVectorRV< V, M > & prior() const
Definition: GPMSA.C:1051
unsigned int dimLocal() const
Definition: VectorSpace.C:155
void addSimulations(const std::vector< typename SharedPtr< V >::Type > &simulationScenarios, const std::vector< typename SharedPtr< V >::Type > &simulationParameters, const std::vector< typename SharedPtr< V >::Type > &simulationOutputs)
Adds multiple simulations to this.
Definition: GPMSA.C:723
bool allocated_m_opts
Definition: GPMSA.h:481
A templated base class for handling vector RV.
Definition: VectorRV.h:53
const Map & map() const
Map.
Definition: VectorSpace.C:142
void setUpEmulator()
Definition: GPMSA.C:739
void setDiscrepancyBases(const std::vector< typename SharedPtr< V >::Type > &discrepancyBases)
Add all discrepancy bases to this.
Definition: GPMSA.C:1013
The QUESO MPI Communicator Class.
Definition: MpiComm.h:203
const std::vector< typename SharedPtr< V >::Type > & experimentScenarios() const
Return all points in scenarioSpace for all experiments.
Definition: GPMSA.C:663
const std::vector< typename SharedPtr< V >::Type > & m_simulationOutputs
Definition: GPMSA.h:96
A class representing a vector RV constructed via Gamma distribution.
Definition: GammaVectorRV.h:64
const std::vector< typename SharedPtr< V >::Type > & simulationScenarios() const
Return all points in scenarioSpace for all simulations.
Definition: GPMSA.C:606
const std::vector< typename SharedPtr< V >::Type > & simulationParameters() const
Return all points in parameterSpace for all simulations.
Definition: GPMSA.C:625
const MpiComm & Comm() const
Access function for MpiComm communicator.
Definition: Map.C:131
const VectorSpace< V, M > & scenarioSpace() const
Return the vector space in which scenarios live.
Definition: GPMSA.C:566
const V & simulationOutput(unsigned int simulationId) const
Return the simulation output for simulation simulationId.
Definition: GPMSA.C:632
const V & experimentScenario(unsigned int experimentId) const
Return the point in scenarioSpace for experiment experimentId.
Definition: GPMSA.C:651
This class defines the options that specify the behaviour of the Gaussian process emulator...
Definition: GPMSAOptions.h:45
const VectorSpace< V, M > & m_experimentOutputSpace
Definition: GPMSA.h:341
virtual ~GPMSAEmulator()
Definition: GPMSA.C:89
const V & simulationScenario(unsigned int simulationId) const
Return the point in scenarioSpace for simulation simulationId.
Definition: GPMSA.C:594
RawType_MPI_Comm Comm() const
Extract MPI Communicator from a MpiComm object.
Definition: MpiComm.C:111

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5