#include <GPMSA.h>

Inheritance diagram for QUESO::GPMSAEmulator< V, M >:

Collaboration diagram for QUESO::GPMSAEmulator< V, M >:

Public Member Functions
	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< V * > &m_simulationScenarios, const std::vector< V * > &m_simulationParameters, const std::vector< V * > &m_simulationOutputs, const std::vector< V * > &m_experimentScenarios, const std::vector< V * > &m_experimentOutputs, const std::vector< V > &m_discrepancyBases, const std::vector< M > &m_observationErrorMatrices, const M &m_experimentErrors, const ConcatenatedVectorRV< V, M > &m_totalPrior, const V &residual_in, const M &BT_Wy_B_inv_in, const M &KT_K_inv_in)

virtual	~GPMSAEmulator ()

virtual double	lnValue (const V &domainVector, const V domainDirection, V gradVector, M hessianMatrix, V hessianEffect) const
	Logarithm of the value of the scalar function. More...

virtual double	actualValue (const V &domainVector, const V domainDirection, V gradVector, M hessianMatrix, V hessianEffect) const
	Actual value of the scalar function. More...

Public Member Functions inherited from QUESO::BaseScalarFunction< V, M >
	BaseScalarFunction (const char *prefix, const VectorSet< V, M > &domainSet)
	Default constructor. More...

virtual	~BaseScalarFunction ()
	Destructor. More...

const VectorSet< V, M > &	domainSet () const
	Access to the protected attribute `m_domainSet:` domain set of the scalar function. More...

Public Attributes
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< V * > &	m_simulationScenarios

const std::vector< V * > &	m_simulationParameters

const std::vector< V * > &	m_simulationOutputs

const std::vector< V * > &	m_experimentScenarios

const std::vector< V * > &	m_experimentOutputs

std::vector< V >	m_discrepancyBases

const std::vector< M > &	m_observationErrorMatrices

unsigned int	num_svd_terms

const M &	m_experimentErrors

const ConcatenatedVectorRV< V, M > &	m_totalPrior

const V &	residual

const M &	BT_Wy_B_inv

const M &	KT_K_inv

Additional Inherited Members
Protected Attributes inherited from QUESO::BaseScalarFunction< V, M >
const BaseEnvironment &	m_env

std::string	m_prefix

const VectorSet< V, M > &	m_domainSet
	Domain set of the scalar function. More...

Detailed Description

template<class V = GslVector, class M = GslMatrix>
class QUESO::GPMSAEmulator< V, M >

Definition at line 47 of file GPMSA.h.

Constructor & Destructor Documentation

template<class V , class M >

QUESO::GPMSAEmulator< V, M >::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< V * > &	m_simulationScenarios,
		const std::vector< V * > &	m_simulationParameters,
		const std::vector< V * > &	m_simulationOutputs,
		const std::vector< V * > &	m_experimentScenarios,
		const std::vector< V * > &	m_experimentOutputs,
		const std::vector< V > &	m_discrepancyBases,
		const std::vector< M > &	m_observationErrorMatrices,
		const M &	m_experimentErrors,
		const ConcatenatedVectorRV< V, M > &	m_totalPrior,
		const V &	residual_in,
		const M &	BT_Wy_B_inv_in,
		const M &	KT_K_inv_in
	)

Definition at line 32 of file GPMSA.C.

References QUESO::VectorSpace< V, M >::dimLocal(), QUESO::GPMSAEmulator< V, M >::num_svd_terms, queso_assert_equal_to, and queso_assert_greater.

   :
   BaseScalarFunction<V, M>("", m_totalPrior.imageSet()),
   m_scenarioSpace(m_scenarioSpace),
   m_parameterSpace(m_parameterSpace),
   m_simulationOutputSpace(m_simulationOutputSpace),
   m_experimentOutputSpace(m_experimentOutputSpace),
   m_numSimulations(m_numSimulations),
   m_numExperiments(m_numExperiments),
   m_simulationScenarios(m_simulationScenarios),
   m_simulationParameters(m_simulationParameters),
   m_simulationOutputs(m_simulationOutputs),
   m_experimentScenarios(m_experimentScenarios),
   m_experimentOutputs(m_experimentOutputs),
   m_discrepancyBases(m_discrepancyBases),
   m_observationErrorMatrices(m_observationErrorMatrices),
   m_experimentErrors(m_experimentErrors),
   m_totalPrior(m_totalPrior),
   residual(residual_in),
   BT_Wy_B_inv(BT_Wy_B_inv_in),
   KT_K_inv(KT_K_inv_in)
 {
   queso_assert_greater(m_numSimulations, 0);
 
   queso_assert_equal_to
     (m_simulationOutputs[0]->map().Comm().NumProc(), 1);
 
   const unsigned int numOutputs =
     this->m_experimentOutputSpace.dimLocal();
 
   const unsigned int MAX_SVD_TERMS =
     std::min(m_numSimulations,(unsigned int)(5));
   num_svd_terms = std::min(MAX_SVD_TERMS, numOutputs);
 }

template<class V , class M >

QUESO::GPMSAEmulator< V, M >::~GPMSAEmulator ( )

virtual

Definition at line 87 of file GPMSA.C.

 {
   // heap allocations get deleted via ScopedPtr
 }

Member Function Documentation

template<class V , class M >

double QUESO::GPMSAEmulator< V, M >::actualValue	(	const V &	domainVector,
		const V *	domainDirection,
		V *	gradVector,
		M *	hessianMatrix,
		V *	hessianEffect
	)		const

virtual

Actual value of the scalar function.

Implements QUESO::BaseScalarFunction< V, M >.

Definition at line 414 of file GPMSA.C.

 {
   // Do nothing?
   return 1.0;
 }

template<class V , class M >

double QUESO::GPMSAEmulator< V, M >::lnValue	(	const V &	domainVector,
		const V *	domainDirection,
		V *	gradVector,
		M *	hessianMatrix,
		V *	hessianEffect
	)		const

virtual

Logarithm of the value of the scalar function.

Implements QUESO::BaseScalarFunction< V, M >.

Definition at line 94 of file GPMSA.C.

References QUESO::MpiComm::Comm(), k, queso_assert, queso_assert_greater, queso_assert_greater_equal, queso_error, and QUESO::queso_isnan().

 {
   // Components of domainVector:
   // theta(1)                     // = "theta", "t" in Higdon et. al. 2008
   // theta(2)
   // ...
   // theta(dimParameter)
   // emulator_mean                // = "mu"
   // emulator_precision           // = "lambda_eta"
   // weights_precision(1)         // = "lambda_{wi}" in vector case
   // ...                          // nonexistent in scalar case
   // weights_precision(num_svd_terms)
   // emulator_corr_strength(1)    // = "rho_{eta k}"
   // ...                          // dimScenario = "p_x", dimParameter = "p_t"
   // emulator_corr_strength(dimScenario + dimParameter)
   // discrepancy_precision(1)     // = "lambda_delta" in scalar case,
   //                              //   "lambda_{vi}" in vector
   // ...
   // discrepancy_precision(F)     // FIXME: F = 1, |G_1| = p_delta, for now.
   // discrepancy_corr_strength(1) // = "rho_{deltak}"
   // ...
   // discrepancy_corr_strength(dimScenario)
   // emulator_data_precision(1)   // = "small white noise", "small ridge"
   // observation_error_precision  // = "lambda_y", only in vector case
 
   // Other variables:
   // m_numSimulations             // = "m"
   // m_numExperiments             // = "n"
   // m_simulationScenarios        // = "eta"
   // m_experimentScenarios        // = "y"
   // m_experimentErrors           // = "Sigma_y"
   // numOutputs                   // = "n_eta"
   //                              // "n_y" := sum(n_y_i)
   //                              //      (== n*n_eta for us for now)
   // dimScenario                  // = "p_x"
   // num_svd_terms                // = "p_eta"
   // num_discrepancy_bases        // = "p_delta"
   // m_TruncatedSVD_simulationOutputs  // = "K_eta"
   // covMatrix                    // = "Sigma_D" in scalar case,
   //                                   "Sigma_zhat" in vector
   // m_discrepancyMatrices        // = "D_i"
   // m_observationErrorMatrices   // = "W_i"
   // m_observationErrorMatrix     // = "W_y"
   // m_BMatrix                    // = "B"
   // num_discrepancy_groups       // = "F"
   // m_emulatorPrecisionShapeVec          // = "a_eta"
   // 1.0/m_emulatorPrecisionScaleVec      // = "b_eta"
   // m_observationalPrecisionShapeVec     // = "a_y"
   // 1.0/m_observationalPrecisionScaleVec // = "b_y"
 
   // Construct covariance matrix
   const unsigned int totalRuns = this->m_numExperiments + this->m_numSimulations;
   const unsigned int numOutputs = this->m_experimentOutputSpace.dimLocal();
   const unsigned int totalOutputs = totalRuns * numOutputs;
   const unsigned int num_discrepancy_bases = m_discrepancyBases.size();
   const unsigned int residualSize = (numOutputs == 1) ?  totalOutputs :
     totalRuns * num_svd_terms + m_numExperiments * num_discrepancy_bases;
 
   double prodScenario = 1.0;
   double prodParameter = 1.0;
   double prodDiscrepancy = 1.0;
   unsigned int dimScenario = (this->m_scenarioSpace).dimLocal();
   unsigned int dimParameter = (this->m_parameterSpace).dimLocal();
 
   // Length of prior+hyperprior inputs
   unsigned int dimSum = 3 +
                         (numOutputs > 1) * 2 +
                         num_svd_terms +
                         dimParameter +
                         dimParameter +
                         dimScenario +
                         dimScenario;  // yum
 
   // Offset for Sigma_eta equivalent in vector case
   const unsigned int offset1 = (numOutputs == 1) ?
     0 : m_numExperiments * num_discrepancy_bases;
 
   // Offset for Sigma_w in vector case
   const unsigned int offset1b = offset1 +
     m_numExperiments * num_svd_terms;
 
   // Offset for lambda_eta term in zhat covariance in vector case
   const unsigned int offset2 = (numOutputs == 1) ?
     0 : m_numExperiments * (num_discrepancy_bases + num_svd_terms);
 
   // This is cumbersome.  All I want is a matrix.
   const MpiComm & comm = domainVector.map().Comm();
   Map z_map(residualSize, 0, comm);
   M covMatrix(this->m_env, z_map, residualSize);
 
   V domainVectorParameter(*(this->m_simulationParameters[0]));
   for (unsigned int k = 0; k < dimParameter; k++) {
     queso_assert (!queso_isnan(domainVector[k]));
     domainVectorParameter[k] = domainVector[k];
   }
 
   // This for loop is a disaster and could do with a *lot* of optimisation
   for (unsigned int i = 0; i < totalRuns; i++) {
 
      const V* scenario1;
      const V* parameter1;
 
     // Decide whether to do experiment part of the covariance matrix
     // Get i-th simulation-parameter pair
     if (i < this->m_numExperiments) {
       // Experiment scenario (known)
       scenario1 = (this->m_experimentScenarios)[i];
 
       // Experiment parameter (unknown)
       parameter1 = &domainVectorParameter;
     }
     else {
       scenario1 =
         (this->m_simulationScenarios)[i-this->m_numExperiments];
       parameter1 =
         (this->m_simulationParameters)[i-this->m_numExperiments];
     }
 
     for (unsigned int j = 0; j < totalRuns; j++) {
 
        const V* scenario2;
        const V* parameter2;
 
       if (j < this->m_numExperiments) {
         scenario2 = (this->m_experimentScenarios)[j];
         parameter2 = &domainVectorParameter;
       }
       else {
         scenario2 =
           (this->m_simulationScenarios)[j-this->m_numExperiments];
         parameter2 =
           (this->m_simulationParameters)[j-this->m_numExperiments];
       }
 
       // Emulator component       // = first term in (1)
       prodScenario = 1.0;
       unsigned int emulatorCorrStrStart =
         dimParameter + 1 + num_svd_terms;
       for (unsigned int k = 0; k < dimScenario; k++) {
         const double & emulator_corr_strength =
           domainVector[emulatorCorrStrStart+k];
         prodScenario *= std::pow(emulator_corr_strength,
                                  4.0 * ((*scenario1)[k] - (*scenario2)[k]) *
                                        ((*scenario1)[k] - (*scenario2)[k]));
       }
 
       queso_assert (!queso_isnan(prodScenario));
 
       // = second term in (1)
       prodParameter = 1.0;
       for (unsigned int k = 0; k < dimParameter; k++) {
         queso_assert (!queso_isnan(domainVector[emulatorCorrStrStart+dimScenario+k]));
         queso_assert (!queso_isnan((*parameter1)[k]));
         queso_assert (!queso_isnan((*parameter2)[k]));
         const double & emulator_corr_strength =
           domainVector[emulatorCorrStrStart+dimScenario+k];
         prodParameter *= std::pow(
             emulator_corr_strength,
             4.0 * ((*parameter1)[k] - (*parameter2)[k]) *
                   ((*parameter1)[k] - (*parameter2)[k]));
       }
 
       queso_assert (!queso_isnan(prodParameter));
 
       // Sigma_eta in scalar case,
       // [Sigma_u, Sigma_uw; Sigma_uw^T, Sigma_w] in vector case
       for (unsigned int basis = 0; basis != num_svd_terms; ++basis)
         {
           // coefficient in (1)
           // The relevant precision for Sigma_eta is lambda_eta; for
           // Sigma_uw etc. it's lambda_wi and we skip lambda_eta
           const double relevant_precision =
             domainVector[dimParameter+basis+1+(numOutputs>1)];
           queso_assert_greater(relevant_precision, 0.0);
 
           const unsigned int stridei =
             (i < this->m_numExperiments) ? m_numExperiments : m_numSimulations;
           const unsigned int offseti =
             (i < this->m_numExperiments) ? offset1 : offset1b - m_numExperiments;
           const unsigned int stridej =
             (j < this->m_numExperiments) ? m_numExperiments : m_numSimulations;
           const unsigned int offsetj =
             (j < this->m_numExperiments) ? offset1 : offset1b - m_numExperiments;
 
           covMatrix(offseti+basis*stridei+i,
                     offsetj+basis*stridej+j) =
             prodScenario * prodParameter / relevant_precision;
         }
 
       // If we're in the experiment cross correlation part, need extra
       // foo: Sigma_delta/Sigma_v and Sigma_y
       if (i < this->m_numExperiments && j < this->m_numExperiments) {
         V* cross_scenario1 = (this->m_simulationScenarios)[i];
         V* cross_scenario2 = (this->m_simulationScenarios)[j];
         prodDiscrepancy = 1.0;
         unsigned int discrepancyCorrStrStart = dimParameter +
                                                num_svd_terms +
                                                dimParameter +
                                                dimScenario + 2 +
                                                (numOutputs > 1);
         for (unsigned int k = 0; k < dimScenario; k++) {
           const double & discrepancy_corr_strength =
             domainVector[discrepancyCorrStrStart+k];
           prodDiscrepancy *=
             std::pow(discrepancy_corr_strength, 4.0 *
                      ((*cross_scenario1)[k] - (*cross_scenario2)[k]) *
                      ((*cross_scenario1)[k] - (*cross_scenario2)[k]));
         }
 
         const double discrepancy_precision =
           domainVector[discrepancyCorrStrStart-1];
         queso_assert_greater(discrepancy_precision, 0);
         queso_assert (!queso_isnan(prodDiscrepancy));
 
         // Sigma_delta term from below (3) in univariate case
         // Sigma_v term from p. 576 in multivariate case
         const double R_v = prodDiscrepancy / discrepancy_precision;
         for (unsigned int disc = 0; disc != num_discrepancy_bases;
              ++disc)
           covMatrix(disc*m_numExperiments+i,
                     disc*m_numExperiments+j) += R_v;
 
         // Experimental error comes in via K in the multivariate
         // case, but comes in via Sigma_y in the univariate case here
         if (numOutputs == 1)
           {
             // Sigma_y term from below (3)
             const double experimentalError =
               (this->m_experimentErrors)(i,j);
 
             queso_assert_greater_equal (experimentalError, 0);
 
             covMatrix(i,j) += experimentalError;
           }
       }
     }
 
     // Add small white noise component to diagonal to make stuff +ve def
     // = "small ridge"
     const double emulator_data_precision = domainVector[dimSum-1-(numOutputs>1)];
     queso_assert_greater(emulator_data_precision, 0);
     double nugget = 1.0 / emulator_data_precision;
 
     for (unsigned int disc = 0; disc != num_discrepancy_bases;
          ++disc)
       covMatrix(disc*m_numExperiments+i,
                 disc*m_numExperiments+i) += nugget;
   }
 
   // If we're in the multivariate case, we've built the full Sigma_z
   // matrix; now add the remaining Sigma_zhat terms
   if (numOutputs > 1)
     {
       const double lambda_y = domainVector[dimSum-1];
       const double inv_lambda_y = 1.0/lambda_y;
 
       unsigned int BT_Wy_B_size = BT_Wy_B_inv.numCols();
       for (unsigned int i=0; i != BT_Wy_B_size; ++i)
         for (unsigned int j=0; j != BT_Wy_B_size; ++j)
           covMatrix(i,j) += BT_Wy_B_inv(i,j) * inv_lambda_y;
 
       const double emulator_precision =
         domainVector[dimParameter+1];
       const double inv_emulator_precision = 1.0/emulator_precision;
 
       unsigned int KT_K_size = KT_K_inv.numCols();
       for (unsigned int i=0; i != KT_K_size; ++i)
         for (unsigned int j=0; j != KT_K_size; ++j)
           covMatrix(i+offset2,j+offset2) +=
             KT_K_inv(i,j) * inv_emulator_precision;
     }
 
 
   // Solve covMatrix * sol = residual
   // = Sigma_D^-1 * (D - mu 1) from (3)
   V sol(covMatrix.invertMultiply(residual));
 
   // Premultiply by residual^T as in (3)
   double minus_2_log_lhd = 0.0;
   // There's no dot product function in GslVector.
   for (unsigned int i = 0; i < residualSize; i++) {
     minus_2_log_lhd += sol[i] * residual[i];
   }
 
 if (queso_isnan(minus_2_log_lhd))
   for (unsigned int i = 0; i < residualSize; i++) {
     if (queso_isnan(sol[i]))
       std::cout << "NaN sol[" << i << ']' << std::endl;
     if (queso_isnan(residual[i]))
       std::cout << "NaN residual[" << i << ']' << std::endl;
 
     std::cout << "Covariance Matrix:" << std::endl;
     covMatrix.print(std::cout);
   }
 
 // std::cout << "minus_2_log_lhd = " << minus_2_log_lhd << std::endl;
 
   queso_assert_greater(minus_2_log_lhd, 0);
 
   double cov_det = covMatrix.determinant();
 
   if (cov_det <= 0)
     {
       std::cout << "Non-positive determinant for covMatrix = " << std::endl;
       covMatrix.print(std::cout);
       queso_error();
     }
 
   minus_2_log_lhd += std::log(covMatrix.determinant());
 
   // Multiply by -1/2 coefficient from (3)
   return -0.5 * minus_2_log_lhd;
 }