#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< 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)

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. Deprecated. 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 >
void	setFiniteDifferenceStepSize (double fdStepSize)
	Sets the step size for finite differencing gradients. More...

void	setFiniteDifferenceStepSize (unsigned int i, double fdStepSize)

	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...

virtual double	lnValue (const V &domainVector) const
	Returns the logarithm of the function at `domainVector`. More...

virtual double	lnValue (const V &domainVector, V &gradVector) const
	Returns the logarithm of the function and its gradient at `domainVector`. More...

virtual double	lnValue (const V &domainVector, V &gradVector, const V &domainDirection, V &hessianEffect) const

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< 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

std::vector< typename SharedPtr< V >::Type >	m_discrepancyBases

const std::vector< typename SharedPtr< M >::Type > &	m_observationErrorMatrices

SharedPtr< M >::Type	m_observationErrorMatrix

unsigned int	num_svd_terms

const ConcatenatedVectorRV< V, M > &	m_totalPrior

const V &	residual

const M &	BT_Wy_B_inv

const M &	KT_K_inv

const GPMSAOptions &	m_opts

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 48 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< 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 at line 32 of file GPMSA.C.

References QUESO::VectorSpace< V, M >::dimLocal(), QUESO::GPMSAOptions::m_maxEmulatorBasisVectors, QUESO::GPMSAEmulator< V, M >::m_opts, QUESO::GPMSAEmulator< V, M >::m_simulationOutputs, and QUESO::GPMSAEmulator< V, M >::num_svd_terms.

   :
   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_observationErrorMatrix(m_observationErrorMatrix),
   m_totalPrior(m_totalPrior),
   residual(residual_in),
   BT_Wy_B_inv(BT_Wy_B_inv_in),
   KT_K_inv(KT_K_inv_in),
   m_opts(opts)
 {
   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();
 
   this->num_svd_terms = this->m_opts.m_maxEmulatorBasisVectors ?
     std::min((unsigned int)(this->m_opts.m_maxEmulatorBasisVectors), numOutputs) :
     numOutputs;
 }

template<class V , class M >

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

virtual

Definition at line 89 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 443 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. Deprecated.

Pointers will be NULL if derivative information is not required by QUESO.

Default implementation throws an exception.

Reimplemented from QUESO::BaseScalarFunction< V, M >.

Definition at line 96 of file GPMSA.C.

References QUESO::MpiComm::Comm(), k, and QUESO::queso_isnan().

 {
   // Components of domainVector:
   // theta(1)                     // = "theta", "t" in Higdon et. al. 2008
   // theta(2)
   // ...
   // theta(dimParameter)
   // truncation_error_precision   // = "lambda_eta", only exists in
   //                                   vector case
   // emulator_precision(1)        // = "lambda_{w i}" in vector case
   // ...                          // = "lambda_eta" in scalar case,
   //                                   which is unrelated to the
   //                                   vector case lambda_eta
   // emulator_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_{v i}" in vector
   // ...
   // discrepancy_precision(F)     // FIXME: F = 1, |G_1| = p_delta, for now.
   // discrepancy_corr_strength(1) // = "rho_{delta k}" in scalar case,
   // ...                          //   "rho_{v i}" in vector
   // discrepancy_corr_strength(dimScenario)
   // emulator_data_precision(1)   // = "small white noise", "small ridge"
   // observation_error_precision  // = "lambda_y", iff user-requested
 
   // Other variables:
   // m_numSimulations             // = "m"
   // m_numExperiments             // = "n"
   // m_simulationScenarios        // = "eta"
   // m_experimentScenarios        // = "y"
   // m_experimentErrors           // = "Sigma_y"; now obsoleted by "W_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 = 2 +
                         (this->num_svd_terms < numOutputs) +
                         (numOutputs > 1) +
                         m_opts.m_calibrateObservationalPrecision +
                         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);
 
   typename SharedPtr<V>::Type domainVectorParameter
     (new V(*(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++) {
 
      // Scenario and uncertain input variables, *not* normalized
      typename SharedPtr<V>::Type scenario1;
      typename SharedPtr<V>::Type 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++) {
 
        // Scenario and uncertain input variables, *not* normalized
        typename SharedPtr<V>::Type scenario2;
        typename SharedPtr<V>::Type 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 + (this->num_svd_terms < numOutputs) + num_svd_terms;
       for (unsigned int k = 0; k < dimScenario; k++) {
         const double & emulator_corr_strength =
           domainVector[emulatorCorrStrStart+k];
         double scenario_param1 =
           m_opts.normalized_scenario_parameter(k, (*scenario1)[k]);
         double scenario_param2 =
           m_opts.normalized_scenario_parameter(k, (*scenario2)[k]);
         prodScenario *= std::pow(emulator_corr_strength,
                                  4.0 * (scenario_param1 - scenario_param2) *
                                        (scenario_param1 - scenario_param2));
       }
 
       queso_assert (!queso_isnan(prodScenario));
 
       // = second term in (1)
       prodParameter = 1.0;
       for (unsigned int k = 0; k < dimParameter; k++) {
         queso_assert (!queso_isnan(domainVector[emulatorCorrStrStart+dimScenario+k]));
         queso_assert (!queso_isnan((*parameter1)[k]));
         queso_assert (!queso_isnan((*parameter2)[k]));
         const double & emulator_corr_strength =
           domainVector[emulatorCorrStrStart+dimScenario+k];
         double uncertain_param1 =
           m_opts.normalized_uncertain_parameter(k, (*parameter1)[k]);
         double uncertain_param2 =
           m_opts.normalized_uncertain_parameter(k, (*parameter2)[k]);
         prodParameter *= std::pow(
             emulator_corr_strength,
             4.0 * (uncertain_param1 - uncertain_param2) *
                   (uncertain_param1 - uncertain_param2));
       }
 
       queso_assert (!queso_isnan(prodParameter));
 
       // Sigma_eta in scalar case,
       // [Sigma_u, Sigma_uw; Sigma_uw^T, Sigma_w] in vector case
       for (unsigned int basis = 0; basis != num_svd_terms; ++basis)
         {
           // 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 +
                          (this->num_svd_terms<numOutputs) + (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) {
         typename SharedPtr<V>::Type cross_scenario1 = (this->m_experimentScenarios)[i];
         typename SharedPtr<V>::Type cross_scenario2 = (this->m_experimentScenarios)[j];
         prodDiscrepancy = 1.0;
         unsigned int discrepancyCorrStrStart =
           dimParameter + num_svd_terms + dimParameter + dimScenario + 1 +
           (this->num_svd_terms<numOutputs) + (numOutputs > 1);
         for (unsigned int k = 0; k < dimScenario; k++) {
           const double & discrepancy_corr_strength =
             domainVector[discrepancyCorrStrStart+k];
           double cross_scenario_param1 =
             m_opts.normalized_scenario_parameter(k, (*cross_scenario1)[k]);
           double cross_scenario_param2 =
             m_opts.normalized_scenario_parameter(k, (*cross_scenario2)[k]);
           prodDiscrepancy *=
             std::pow(discrepancy_corr_strength, 4.0 *
                      (cross_scenario_param1 - cross_scenario_param2) *
                      (cross_scenario_param1 - cross_scenario_param2));
         }
 
         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;
 
         if (numOutputs == 1)
           {
             // Experimental error comes in via W_y now.
 /*
             // Sigma_y term from below (3)
             const double experimentalError =
               this->m_experimentErrors(i,j);
 */
             const double experimentalError =
               (*this->m_observationErrorMatrix)(i,j);
 
             queso_assert_greater_equal (experimentalError, 0);
 
             const double lambda_y =
               m_opts.m_calibrateObservationalPrecision ?
               domainVector[dimSum-1] : 1.0;
 
             covMatrix(i,j) += lambda_y * 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 =
         m_opts.m_calibrateObservationalPrecision ?
         domainVector[dimSum-1] : 1.0;
       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;
 }