#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 M &m_experimentErrors, const ConcatenatedVectorRV< V, M > &m_totalPrior)

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

const M &	m_experimentErrors

const ConcatenatedVectorRV< V, M > &	m_totalPrior

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 46 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 M &	m_experimentErrors,
		const ConcatenatedVectorRV< V, M > &	m_totalPrior
	)

Definition at line 32 of file GPMSA.C.

   :
   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_experimentErrors(m_experimentErrors),
   m_totalPrior(m_totalPrior)
 {
 }

template<class V , class M >

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

virtual

Definition at line 66 of file GPMSA.C.

 {
   // Do nothing?
 }

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 232 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 73 of file GPMSA.C.

References k.

 {
   // Components of domainVector:
   // theta(1)
   // theta(2)
   // ...
   // theta(dimParameterSpace)
   // emulator_mean
   // emulator_precision
   // emulator_corr_strength(1)
   // ...
   // emulator_corr_strength(dimScenario + dimParameter)
   // discrepancy_precision
   // discrepancy_corr_strength(1)
   // ...
   // discrepancy_corr_strength(dimScenario)
   // emulator_data_precision(1)
 
   // Construct covariance matrix
   unsigned int totalDim = this->m_numExperiments + this->m_numSimulations;
   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();
 
   // This is cumbersome.  All I want is a matrix.
   VectorSpace<V, M> gpSpace(this->m_scenarioSpace.env(), "", totalDim, NULL);
   V residual(gpSpace.zeroVector());
   M covMatrix(residual);
 
   V *scenario1;
   V *scenario2;
   V *parameter1;
   V *parameter2;
 
   // This for loop is a disaster and could do with a *lot* of optimisation
   for (unsigned int i = 0; i < totalDim; i++) {
     for (unsigned int j = 0; j < totalDim; j++) {
       // 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 = new V(*((this->m_experimentScenarios)[i]));
 
         // Experiment parameter (unknown)
         parameter1 = new V(*((this->m_simulationParameters)[0]));
         for (unsigned int k = 0; k < dimParameter; k++) {
           (*parameter1)[k] = domainVector[k];
         }
       }
       else {
         scenario1 =
           new V(*((this->m_simulationScenarios)[i-this->m_numExperiments]));
         parameter1 =
           new V(*((this->m_simulationParameters)[i-this->m_numExperiments]));
       }
 
       if (j < this->m_numExperiments) {
         scenario2 = new V(*((this->m_experimentScenarios)[j]));
         parameter2 = new V(*((this->m_simulationParameters)[0]));
         for (unsigned int k = 0; k < dimParameter; k++) {
           (*parameter2)[k] = domainVector[k];
         }
       }
       else {
         scenario2 =
           new V(*((this->m_simulationScenarios)[j-this->m_numExperiments]));
         parameter2 =
           new V(*((this->m_simulationParameters)[j-this->m_numExperiments]));
       }
 
       // Emulator component
       prodScenario = 1.0;
       prodParameter = 1.0;
       unsigned int emulatorCorrStrStart = dimParameter + 2;
       for (unsigned int k = 0; k < dimScenario; k++) {
         prodScenario *= std::pow(domainVector[emulatorCorrStrStart+k],
                                  4.0 * ((*scenario1)[k] - (*scenario2)[k]) *
                                        ((*scenario1)[k] - (*scenario2)[k]));
       }
 
       for (unsigned int k = 0; k < dimParameter; k++) {
         prodParameter *= std::pow(
             domainVector[emulatorCorrStrStart+dimScenario+k],
             4.0 * ((*parameter1)[k] - (*parameter2)[k]) *
                   ((*parameter1)[k] - (*parameter2)[k]));
       }
 
       double emPrecision = domainVector[dimParameter+1];
       covMatrix(i, j) = prodScenario * prodParameter /
                         emPrecision;  // emulator precision
 
       delete scenario1;
       delete scenario2;
       delete parameter1;
       delete parameter2;
 
       // If we're in the experiment cross correlation part, need extra foo
       if (i < this->m_numExperiments && j < this->m_numExperiments) {
         scenario1 = new V(*((this->m_simulationScenarios)[i]));
         scenario2 = new V(*((this->m_simulationScenarios)[j]));
         prodDiscrepancy = 1.0;
         unsigned int discrepancyCorrStrStart = dimParameter +
                                                dimParameter +
                                                dimScenario + 3;
         for (unsigned int k = 0; k < dimScenario; k++) {
           prodDiscrepancy *= std::pow(domainVector[discrepancyCorrStrStart+k],
                                       4.0 * ((*scenario1)[k] - (*scenario2)[k]) *
                                             ((*scenario1)[k] - (*scenario2)[k]));
         }
 
         covMatrix(i, j) += prodDiscrepancy /
                            domainVector[discrepancyCorrStrStart-1];
         covMatrix(i, j) += (this->m_experimentErrors)(i, j);
 
         delete scenario1;
         delete scenario2;
       }
     }
 
     // Add small white noise component to diagonal to make stuff +ve def
     unsigned int dimSum = 4 +
                           dimParameter +
                           dimParameter +
                           dimScenario +
                           dimScenario;  // yum
     double nugget = 1.0 / domainVector[dimSum-1];
     covMatrix(i, i) += nugget;
   }
 
   // Form residual = D - mean
   for (unsigned int i = 0; i < this->m_numExperiments; i++) {
     // Scalar so ok -- will need updating for nonscalar case
     residual[i] = (*((this->m_experimentOutputs)[i]))[0];
   }
   for (unsigned int i = 0; i < this->m_numSimulations; i++) {
     // Scalar so ok -- will need updating for nonscalar case
     residual[i+this->m_numExperiments] = (*((this->m_simulationOutputs)[i]))[0];
   }
 
   // Solve covMatrix * sol = residual
   V sol(covMatrix.invertMultiply(residual));
 
   // There's no dot product function in GslVector.
   double minus_2_log_lhd = 0.0;
   for (unsigned int i = 0; i < totalDim; i++) {
     minus_2_log_lhd += sol[i] * residual[i];
   }
 
   return -0.5 * minus_2_log_lhd;
 }