QUESO::StatisticalInverseProblem< P_V, P_M > Class Template Reference

This templated class represents a Statistical Inverse Problem. More...

#include <StatisticalInverseProblem.h>

Collaboration diagram for QUESO::StatisticalInverseProblem< P_V, P_M >:

Public Member Functions
Constructor/Destructor methods
	StatisticalInverseProblem (const char *prefix, const SipOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > &priorRv, const BaseScalarFunction< P_V, P_M > &likelihoodFunction, GenericVectorRV< P_V, P_M > &postRv)
	Constructor. More...
	StatisticalInverseProblem (const char *prefix, const SipOptionsValues *alternativeOptionsValues, const GPMSAFactory< P_V, P_M > &gpmsaFactory, GenericVectorRV< P_V, P_M > &postRv)
	Constructor for statistical inverse problems to be sovled using GPMSA. More...
	~StatisticalInverseProblem ()
	Destructor. More...
Statistical methods
bool	computeSolutionFlag () const
	Whether or not compute the solution. More...
void	solveWithBayesMetropolisHastings (const MhOptionsValues *alternativeOptionsValues, const P_V &initialValues, const P_M *initialProposalCovMatrix)
	Solves the problem via Bayes formula and a Metropolis-Hastings algorithm. More...
void	seedWithMAPEstimator ()
	Seeds the chain with the result of a deterministic optimisation. More...
void	solveWithBayesMLSampling ()
	Solves with Bayes Multi-Level (ML) sampling. More...
const MetropolisHastingsSG < P_V, P_M > &	sequenceGenerator () const
	Return the underlying MetropolisHastingSG object. More...
const BaseVectorRV< P_V, P_M > &	priorRv () const
	Returns the Prior RV; access to private attribute m_priorRv. More...
const GenericVectorRV< P_V, P_M > &	postRv () const
	Returns the Posterior RV; access to private attribute m_postrRv. More...
double	logEvidence () const
	Returns the logarithm value of the evidence. Related to ML. More...
double	meanLogLikelihood () const
	Returns the mean of the logarithm value of the likelihood. Related to ML. More...
double	eig () const

Private Attributes
const BaseEnvironment &	m_env
const BaseVectorRV< P_V, P_M > &	m_priorRv
const BaseScalarFunction< P_V, P_M > &	m_likelihoodFunction
GenericVectorRV< P_V, P_M > &	m_postRv
VectorSet< P_V, P_M > *	m_solutionDomain
BaseJointPdf< P_V, P_M > *	m_solutionPdf
BaseVectorMdf< P_V, P_M > *	m_subSolutionMdf
BaseVectorCdf< P_V, P_M > *	m_subSolutionCdf
BaseVectorRealizer< P_V, P_M > *	m_solutionRealizer
MetropolisHastingsSG< P_V, P_M > *	m_mhSeqGenerator
MLSampling< P_V, P_M > *	m_mlSampler
BaseVectorSequence< P_V, P_M > *	m_chain
ScalarSequence< double > *	m_logLikelihoodValues
ScalarSequence< double > *	m_logTargetValues
SipOptionsValues	m_alternativeOptionsValues
StatisticalInverseProblemOptions *	m_optionsObj
bool	m_seedWithMAPEstimator

I/O methods
void	print (std::ostream &os) const
	TODO: Prints the sequence. More...
std::ostream &	operator<< (std::ostream &os, const StatisticalInverseProblem< P_V, P_M > &obj)

Detailed Description

template<class P_V, class P_M>
class QUESO::StatisticalInverseProblem< P_V, P_M >

This templated class represents a Statistical Inverse Problem.

This class is templated on the type 'P_V' of vector and type 'P_M' of matrix, where 'P_' stands for 'parameter'. Some templated classes might also use 'Q_V' and 'Q_M' when referring to a vector type and a matrix type, where 'Q_' stands for 'quantities of interest'. Conceptually, a statistical inverse problem has two input entities and one output entity.
The input entities of a statistical inverse problem are:

1. the prior RV, an instance of class 'BaseVectorRV<P_V,P_M>', and
2. the likelihood function, an instance of class 'BaseScalarFunction<P_V,P_M>'.

Let denote the mathematical likelihood function and denote a vector of parameters. The likelihood function object stores the routine that computes and whatever data necessary by such routine. The routine in the likelihood function object can compute either the actual value or the value . See files 'basic/inc/ScalarFunction.h' and 'stats/inc/JointPdf.h' for more details.
The output entity of a statistical inverse problem is the posterior RV, another instance of class 'BaseVectorRV<P_V,P_M>', which stores the solution according to the Bayesian approach. Upon return from a solution operation, the posterior RV is available through the operation 'postRv()'. Such posterior RV is able to provide:

1. a joint pdf (up to a multiplicative constant) through the operation 'postRv().pdf()', which returns an instance of the class 'BaseJointPdf<P_V,P_M>', and
2. a vector realizer through the operation 'postRv().realizer()', which returns an instance of the class 'BaseVectorRealizer<P_V,P_M>'.

Definition at line 81 of file StatisticalInverseProblem.h.

Constructor & Destructor Documentation

template<class P_V , class P_M >

QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem	(	const char *	prefix,
		const SipOptionsValues *	alternativeOptionsValues,
		const BaseVectorRV< P_V, P_M > &	priorRv,
		const BaseScalarFunction< P_V, P_M > &	likelihoodFunction,
		GenericVectorRV< P_V, P_M > &	postRv
	)

Constructor.

Requirements: 1) the image set of the vector random variable 'priorRv', 2) the domain set of the likelihood function 'likelihoodFunction' and 3) the image set of the vector random variable 'postRv' should belong to vector spaces of equal dimensions. If the requirements are satisfied, the constructor then reads input options that begin with the string '<prefix>ip_'. If no options input file is provided, the construction assigns alternativeOptionsValues to the options of the SIP.

Parameters

prefix	The prefix
alternativeOptionsValues	Options (if no input file)
priorRv	The prior RV
likelihoodFunction	The likelihood function
postRv	The posterior RV

Definition at line 35 of file StatisticalInverseProblem.C.

References UQ_FATAL_TEST_MACRO.

  :
  m_env                     (priorRv.env()),
  m_priorRv                 (priorRv),
  m_likelihoodFunction      (likelihoodFunction),
  m_postRv                  (postRv),
  m_solutionDomain          (NULL),
  m_solutionPdf             (NULL),
  m_subSolutionMdf          (NULL),
  m_subSolutionCdf          (NULL),
  m_solutionRealizer        (NULL),
  m_mhSeqGenerator          (NULL),
  m_mlSampler               (NULL),
  m_chain                   (NULL),
  m_logLikelihoodValues     (NULL),
  m_logTargetValues         (NULL),
  m_alternativeOptionsValues(),
  m_optionsObj              (NULL),
  m_seedWithMAPEstimator    (false)
{
#ifdef QUESO_MEMORY_DEBUGGING
  std::cout << "Entering Sip" << std::endl;
#endif
  if (m_env.subDisplayFile()) {
    *m_env.subDisplayFile() << "Entering StatisticalInverseProblem<P_V,P_M>::constructor()"
                            << ": prefix = " << prefix
                            << ", alternativeOptionsValues = " << alternativeOptionsValues
                            << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
                            << std::endl;
  }

  if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
  if (m_env.optionsInputFileName() == "") {
    m_optionsObj = new StatisticalInverseProblemOptions(m_env,prefix,m_alternativeOptionsValues);
  }
  else {
    m_optionsObj = new StatisticalInverseProblemOptions(m_env,prefix);
    m_optionsObj->scanOptionsValues();
  }
#ifdef QUESO_MEMORY_DEBUGGING
  std::cout << "In Sip, finished scanning options" << std::endl;
#endif

  UQ_FATAL_TEST_MACRO(priorRv.imageSet().vectorSpace().dimLocal() != likelihoodFunction.domainSet().vectorSpace().dimLocal(),
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::constructor()",
                      "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");

  UQ_FATAL_TEST_MACRO(priorRv.imageSet().vectorSpace().dimLocal() != postRv.imageSet().vectorSpace().dimLocal(),
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::constructor()",
                      "'priorRv' and 'postRv' are related to vector spaces of different dimensions");

  if (m_env.subDisplayFile()) {
    *m_env.subDisplayFile() << "Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
                            << ": prefix = " << m_optionsObj->m_prefix
                            << std::endl;
  }

  return;
}

template<class P_V , class P_M >

QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem	(	const char *	prefix,
		const SipOptionsValues *	alternativeOptionsValues,
		const GPMSAFactory< P_V, P_M > &	gpmsaFactory,
		GenericVectorRV< P_V, P_M > &	postRv
	)

Constructor for statistical inverse problems to be sovled using GPMSA.

Requirements: 1) the factory for the GPMSA object 2) the image set of the vector random variable postRv (obtainable via GaussianProcessFactory::prior method) should be equal to the full prior image set (including all the hyperparameters)

If the requirements are satisfied, the constructor then reads input options that begin with the string '<prefix>ip_'. If no options input file is provided, the construction assigns alternativeOptionsValues to the options of the statistical inverse problem.

Definition at line 103 of file StatisticalInverseProblem.C.

References QUESO::VectorSpace< V, M >::dimLocal(), QUESO::BaseScalarFunction< V, M >::domainSet(), QUESO::BaseVectorRV< V, M >::imageSet(), QUESO::StatisticalInverseProblem< P_V, P_M >::m_alternativeOptionsValues, QUESO::StatisticalInverseProblem< P_V, P_M >::m_env, QUESO::StatisticalInverseProblem< P_V, P_M >::m_likelihoodFunction, QUESO::StatisticalInverseProblem< P_V, P_M >::m_optionsObj, QUESO::StatisticalInverseProblemOptions::m_prefix, QUESO::StatisticalInverseProblem< P_V, P_M >::m_priorRv, QUESO::BaseEnvironment::optionsInputFileName(), QUESO::StatisticalInverseProblemOptions::scanOptionsValues(), QUESO::BaseEnvironment::subDisplayFile(), UQ_FATAL_TEST_MACRO, QUESO::VectorSet< V, M >::vectorSpace(), and QUESO::BaseEnvironment::worldRank().

  :
  m_env                     (gpmsaFactory.m_totalPrior->env()),
  m_priorRv                 (*(gpmsaFactory.m_totalPrior)),
  m_likelihoodFunction      (gpmsaFactory.getGPMSAEmulator()),
  m_postRv                  (postRv),
  m_solutionDomain          (NULL),
  m_solutionPdf             (NULL),
  m_subSolutionMdf          (NULL),
  m_subSolutionCdf          (NULL),
  m_solutionRealizer        (NULL),
  m_mhSeqGenerator          (NULL),
  m_mlSampler               (NULL),
  m_chain                   (NULL),
  m_logLikelihoodValues     (NULL),
  m_logTargetValues         (NULL),
  m_alternativeOptionsValues(),
  m_optionsObj              (NULL),
  m_seedWithMAPEstimator    (false)
{
  if (m_env.subDisplayFile()) {
    *m_env.subDisplayFile() << "Entering StatisticalInverseProblem<P_V,P_M>::constructor()"
                            << ": prefix = " << prefix
                            << ", alternativeOptionsValues = " << alternativeOptionsValues
                            << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
                            << std::endl;
  }

  if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
  if (m_env.optionsInputFileName() == "") {
    m_optionsObj = new StatisticalInverseProblemOptions(m_env,prefix,m_alternativeOptionsValues);
  }
  else {
    m_optionsObj = new StatisticalInverseProblemOptions(m_env,prefix);
    m_optionsObj->scanOptionsValues();
  }

  UQ_FATAL_TEST_MACRO(m_priorRv.imageSet().vectorSpace().dimLocal() != m_likelihoodFunction.domainSet().vectorSpace().dimLocal(),
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::constructor()",
                      "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");

  UQ_FATAL_TEST_MACRO(m_priorRv.imageSet().vectorSpace().dimLocal() != postRv.imageSet().vectorSpace().dimLocal(),
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::constructor()",
                      "'priorRv' and 'postRv' are related to vector spaces of different dimensions");

  if (m_env.subDisplayFile()) {
    *m_env.subDisplayFile() << "Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
                            << ": prefix = " << m_optionsObj->m_prefix
                            << std::endl;
  }

  return;
}

template<class P_V , class P_M >

QUESO::StatisticalInverseProblem< P_V, P_M >::~StatisticalInverseProblem

(

)

Destructor.

Definition at line 165 of file StatisticalInverseProblem.C.

{
  if (m_chain) {
    m_chain->clear();
    delete m_chain;
  }
  if (m_logLikelihoodValues) {
    m_logLikelihoodValues->clear();
    delete m_logLikelihoodValues;
  }
  if (m_logTargetValues) {
    m_logTargetValues->clear();
    delete m_logTargetValues;
  }
  if (m_mlSampler       ) delete m_mlSampler;
  if (m_mhSeqGenerator  ) delete m_mhSeqGenerator;
  if (m_solutionRealizer) delete m_solutionRealizer;
  if (m_subSolutionCdf  ) delete m_subSolutionCdf;
  if (m_subSolutionMdf  ) delete m_subSolutionMdf;
  if (m_solutionPdf     ) delete m_solutionPdf;
  if (m_solutionDomain  ) delete m_solutionDomain;
  if (m_optionsObj      ) delete m_optionsObj;
}

Member Function Documentation

template<class P_V, class P_M>

bool QUESO::StatisticalInverseProblem< P_V, P_M >::computeSolutionFlag

(

)

const

Whether or not compute the solution.

template<class P_V , class P_M >

double QUESO::StatisticalInverseProblem< P_V, P_M >::eig

(

)

const

Definition at line 477 of file StatisticalInverseProblem.C.

References UQ_FATAL_TEST_MACRO.

{
  UQ_FATAL_TEST_MACRO(m_mlSampler == NULL,
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::eig()",
                      "m_mlSampler is NULL");
  return m_mlSampler->eig();
}

template<class P_V , class P_M >

double QUESO::StatisticalInverseProblem< P_V, P_M >::logEvidence

(

)

const

Returns the logarithm value of the evidence. Related to ML.

Definition at line 457 of file StatisticalInverseProblem.C.

References UQ_FATAL_TEST_MACRO.

{
  UQ_FATAL_TEST_MACRO(m_mlSampler == NULL,
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::logEvidence()",
                      "m_mlSampler is NULL");
  return m_mlSampler->logEvidence();
}

template<class P_V , class P_M >

double QUESO::StatisticalInverseProblem< P_V, P_M >::meanLogLikelihood

(

)

const

Returns the mean of the logarithm value of the likelihood. Related to ML.

Definition at line 467 of file StatisticalInverseProblem.C.

References UQ_FATAL_TEST_MACRO.

{
  UQ_FATAL_TEST_MACRO(m_mlSampler == NULL,
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::meanLogLikelihood()",
                      "m_mlSampler is NULL");
  return m_mlSampler->meanLogLikelihood();
}

template<class P_V , class P_M >

const GenericVectorRV< P_V, P_M > & QUESO::StatisticalInverseProblem< P_V, P_M >::postRv

(

)

const

Returns the Posterior RV; access to private attribute m_postrRv.

The Posterior RV contains the solution of the Bayes problem.

Definition at line 451 of file StatisticalInverseProblem.C.

{
  return m_postRv;
}

template<class P_V , class P_M >

void QUESO::StatisticalInverseProblem< P_V, P_M >::print

(

std::ostream &

os

)

const

TODO: Prints the sequence.

Todo:

: implement me!

Definition at line 488 of file StatisticalInverseProblem.C.

{
  return;
}

template<class P_V , class P_M >

const BaseVectorRV< P_V, P_M > & QUESO::StatisticalInverseProblem< P_V, P_M >::priorRv

(

)

const

Returns the Prior RV; access to private attribute m_priorRv.

Definition at line 444 of file StatisticalInverseProblem.C.

{
  return m_priorRv;
}

template<class P_V , class P_M >

void QUESO::StatisticalInverseProblem< P_V, P_M >::seedWithMAPEstimator

(

)

Seeds the chain with the result of a deterministic optimisation.

This only works for Metropolis-Hastings right now. Multi-level is not currently supported.

Definition at line 358 of file StatisticalInverseProblem.C.

{
  this->m_seedWithMAPEstimator = true;
}

template<class P_V , class P_M >

const MetropolisHastingsSG< P_V, P_M > & QUESO::StatisticalInverseProblem< P_V, P_M >::sequenceGenerator

(

)

const

Return the underlying MetropolisHastingSG object.

Definition at line 436 of file StatisticalInverseProblem.C.

{
  return *m_mhSeqGenerator;
}

template<class P_V , class P_M >

void QUESO::StatisticalInverseProblem< P_V, P_M >::solveWithBayesMetropolisHastings	(	const MhOptionsValues *	alternativeOptionsValues,
		const P_V &	initialValues,
		const P_M *	initialProposalCovMatrix
	)

Solves the problem via Bayes formula and a Metropolis-Hastings algorithm.

Requirements:

1. 'initialValues' should have the same number of components as member variable 'm_priorRv'
2. if 'initialProposalCovMatrix' is not NULL, it should be square and its size should be equal to the size of 'initialValues'.

If the requirements are satisfied, this methods checks the member flag 'm_computeSolution' (one of the options read from the input file during construction). If the flag is 'false', the operation returns immediately, computing nothing; otherwise, the operation sets the member variable 'm_postRv' accordingly. The operation:

1. sets the pdf of 'm_postRv' equal to an instance of BayesianJointPdf<P_V,P_M>,
2. instantiates SequenceOfVectors<P_V,P_M> (the chain),
3. instantiates MetropolisHastingsSG<P_V,P_M> (the Metropolis-Hastings algorithm),
4. populates the chain with the Metropolis-Hastings algorithm, and
5. sets the realizer of 'm_postRv' with the contents of the chain.

If seedWithMapEstimator() is called before this method, then initialValues is used as the initial condition for an optimization procedure, the result of which will be used as the seed for the chain.

Definition at line 191 of file StatisticalInverseProblem.C.

References QUESO::InstantiateIntersection(), QUESO::GslOptimizer::minimize(), QUESO::GslOptimizer::minimizer(), QUESO::GslOptimizer::setInitialPoint(), UQ_FATAL_TEST_MACRO, and UQ_SIP_FILENAME_FOR_NO_FILE.

{
  //grvy_timer_begin("BayesMetropolisHastings"); TODO: revisit timing output
  //std::cout << "proc " << m_env.fullRank() << ", HERE sip 000" << std::endl;
  m_env.fullComm().Barrier();
  //std::cout << "proc " << m_env.fullRank() << ", HERE sip 001" << std::endl;
  m_env.fullComm().syncPrintDebugMsg("Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);

  if (m_optionsObj->m_ov.m_computeSolution == false) {
    if ((m_env.subDisplayFile())) {
      *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
                              << ": avoiding solution, as requested by user"
                              << std::endl;
    }
    return;
  }
  if ((m_env.subDisplayFile())) {
    *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
                            << ": computing solution, as requested by user"
                            << std::endl;
  }

  UQ_FATAL_TEST_MACRO(m_priorRv.imageSet().vectorSpace().dimLocal() != initialValues.sizeLocal(),
                      m_env.worldRank(),
                      "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
                      "'m_priorRv' and 'initialValues' should have equal dimensions");

  if (initialProposalCovMatrix) {
    UQ_FATAL_TEST_MACRO(m_priorRv.imageSet().vectorSpace().dimLocal() != initialProposalCovMatrix->numRowsLocal(),
                        m_env.worldRank(),
                        "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
                        "'m_priorRv' and 'initialProposalCovMatrix' should have equal dimensions");
    UQ_FATAL_TEST_MACRO(initialProposalCovMatrix->numCols() != initialProposalCovMatrix->numRowsGlobal(),
                        m_env.worldRank(),
                        "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
                        "'initialProposalCovMatrix' should be a square matrix");
  }

  if (m_mlSampler       ) delete m_mlSampler;
  if (m_mhSeqGenerator  ) delete m_mhSeqGenerator;
  if (m_solutionRealizer) delete m_solutionRealizer;
  if (m_subSolutionCdf  ) delete m_subSolutionCdf;
  if (m_subSolutionMdf  ) delete m_subSolutionMdf;
  if (m_solutionPdf     ) delete m_solutionPdf;
  if (m_solutionDomain  ) delete m_solutionDomain;

  P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
  numEvaluationPointsVec.cwSet(250.);

  // Compute output pdf up to a multiplicative constant: Bayesian approach
  m_solutionDomain = InstantiateIntersection(m_priorRv.pdf().domainSet(),m_likelihoodFunction.domainSet());

  m_solutionPdf = new BayesianJointPdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
                                                       m_priorRv.pdf(),
                                                       m_likelihoodFunction,
                                                       1.,
                                                       *m_solutionDomain);

  m_postRv.setPdf(*m_solutionPdf);
  m_chain = new SequenceOfVectors<P_V,P_M>(m_postRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"chain");

  // Decide whether or not to create a MetropolisHastingsSG instance from the
  // user-provided initial seed, or use the user-provided seed for a
  // deterministic optimisation instead and seed the chain with the result of
  // the optimisation
  if (this->m_seedWithMAPEstimator) {
    // Do optimisation before sampling
    GslOptimizer optimizer(*m_solutionPdf);
    optimizer.setInitialPoint(dynamic_cast<const GslVector &>(initialValues));
    optimizer.minimize();

    // Compute output realizer: Metropolis-Hastings approach
    m_mhSeqGenerator = new MetropolisHastingsSG<P_V, P_M>(
        m_optionsObj->m_prefix.c_str(), alternativeOptionsValues,
        m_postRv, optimizer.minimizer(), initialProposalCovMatrix);
  }
  else {
    // Compute output realizer: Metropolis-Hastings approach
    m_mhSeqGenerator = new MetropolisHastingsSG<P_V, P_M>(
        m_optionsObj->m_prefix.c_str(), alternativeOptionsValues, m_postRv,
        initialValues, initialProposalCovMatrix);
  }


  m_logLikelihoodValues = new ScalarSequence<double>(m_env, 0,
                                                     m_optionsObj->m_prefix +
                                                     "logLike");

  m_logTargetValues = new ScalarSequence<double>(m_env, 0,
                                                 m_optionsObj->m_prefix +
                                                 "logTarget");

  // m_logLikelihoodValues and m_logTargetValues may be NULL
  m_mhSeqGenerator->generateSequence(*m_chain, m_logLikelihoodValues,
                                     m_logTargetValues);

  m_solutionRealizer = new SequentialVectorRealizer<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
                                                                    *m_chain);

  m_postRv.setRealizer(*m_solutionRealizer);

  m_env.fullComm().syncPrintDebugMsg("In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings(), code place 1",3,3000000);
  //m_env.fullComm().Barrier();

  // Compute output mdf: uniform sampling approach
#ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
  m_subMdfGrids  = new ArrayOfOneDGrids <P_V,P_M>((m_optionsObj->m_prefix+"Mdf_").c_str(),m_postRv.imageSet().vectorSpace());
  m_subMdfValues = new ArrayOfOneDTables<P_V,P_M>((m_optionsObj->m_prefix+"Mdf_").c_str(),m_postRv.imageSet().vectorSpace());
  m_chain->subUniformlySampledMdf(numEvaluationPointsVec, // input
                                  *m_subMdfGrids,         // output
                                  *m_subMdfValues);       // output
  m_subSolutionMdf = new SampledVectorMdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
                                                          *m_subMdfGrids,
                                                          *m_subMdfValues);
  m_postRv.setMdf(*m_subSolutionMdf);

  if ((m_optionsObj->m_ov.m_dataOutputFileName                       != UQ_SIP_FILENAME_FOR_NO_FILE                    ) &&
      (m_optionsObj->m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_ov.m_dataOutputAllowedSet.end())) {
    if (m_env.subRank() == 0) {
      // Write data output file
      if (m_env.subDisplayFile()) {
        *m_env.subDisplayFile() << "Opening data output file '" << m_optionsObj->m_ov.m_dataOutputFileName
                                << "' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
                                << std::endl;
      }

      // Open file
      // Always write at the end of an eventual pre-existing file
      std::ofstream* ofsvar = new std::ofstream((m_optionsObj->m_ov.m_dataOutputFileName+"_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
      if ((ofsvar            == NULL ) ||
          (ofsvar->is_open() == false)) {
        delete ofsvar;
        ofsvar = new std::ofstream((m_optionsObj->m_ov.m_dataOutputFileName+"_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
      }
      UQ_FATAL_TEST_MACRO((ofsvar && ofsvar->is_open()) == false,
                          m_env.worldRank(),
                          "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
                          "failed to open file");

      m_postRv.mdf().print(*ofsvar);

      // Close file
      //ofsvar->close();
      delete ofsvar;
      if (m_env.subDisplayFile()) {
        *m_env.subDisplayFile() << "Closed data output file '" << m_optionsObj->m_ov.m_dataOutputFileName
                                << "' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
                                << std::endl;
      }
    }
  }
#endif
  if (m_env.subDisplayFile()) {
    *m_env.subDisplayFile() << std::endl;
  }

  m_env.fullComm().syncPrintDebugMsg("Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
  m_env.fullComm().Barrier();
  // grvy_timer_end("BayesMetropolisHastings"); TODO: revisit timers
  return;
}

template<class P_V , class P_M >

void QUESO::StatisticalInverseProblem< P_V, P_M >::solveWithBayesMLSampling

(

)

Solves with Bayes Multi-Level (ML) sampling.

Definition at line 365 of file StatisticalInverseProblem.C.

References QUESO::InstantiateIntersection().

{
  m_env.fullComm().Barrier();
  m_env.fullComm().syncPrintDebugMsg("Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);

  if (m_optionsObj->m_ov.m_computeSolution == false) {
    if ((m_env.subDisplayFile())) {
      *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
                              << ": avoiding solution, as requested by user"
                              << std::endl;
    }
    return;
  }
  if ((m_env.subDisplayFile())) {
    *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
                            << ": computing solution, as requested by user"
                            << std::endl;
  }

  if (m_mlSampler       ) delete m_mlSampler;
  if (m_mhSeqGenerator  ) delete m_mhSeqGenerator;
  if (m_solutionRealizer) delete m_solutionRealizer;
  if (m_subSolutionCdf  ) delete m_subSolutionCdf;
  if (m_subSolutionMdf  ) delete m_subSolutionMdf;
  if (m_solutionPdf     ) delete m_solutionPdf;
  if (m_solutionDomain  ) delete m_solutionDomain;

  P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
  numEvaluationPointsVec.cwSet(250.);

  // Compute output pdf up to a multiplicative constant: Bayesian approach
  m_solutionDomain = InstantiateIntersection(m_priorRv.pdf().domainSet(),m_likelihoodFunction.domainSet());

  m_solutionPdf = new BayesianJointPdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
                                                       m_priorRv.pdf(),
                                                       m_likelihoodFunction,
                                                       1.,
                                                       *m_solutionDomain);

  m_postRv.setPdf(*m_solutionPdf);

  // Compute output realizer: ML approach
  m_chain = new SequenceOfVectors<P_V,P_M>(m_postRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"chain");
  m_mlSampler = new MLSampling<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
                                             //m_postRv,
                                               m_priorRv,
                                               m_likelihoodFunction);
  //                                           initialValues,
  //                                           initialProposalCovMatrix);

  m_mlSampler->generateSequence(*m_chain,
                                NULL,
                                NULL);

  m_solutionRealizer = new SequentialVectorRealizer<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
                                                                    *m_chain);

  m_postRv.setRealizer(*m_solutionRealizer);

  if (m_env.subDisplayFile()) {
    *m_env.subDisplayFile() << std::endl;
  }

  m_env.fullComm().syncPrintDebugMsg("Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
  m_env.fullComm().Barrier();

  return;
}

Friends And Related Function Documentation

template<class P_V, class P_M>

std::ostream& operator<< ( std::ostream &  os, const StatisticalInverseProblem< P_V, P_M > &  obj  )

friend

Definition at line 194 of file StatisticalInverseProblem.h.

                                                     {
    obj.print(os);
    return os;
  }

Member Data Documentation

template<class P_V, class P_M>

SipOptionsValues QUESO::StatisticalInverseProblem< P_V, P_M >::m_alternativeOptionsValues

private

Definition at line 219 of file StatisticalInverseProblem.h.

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem().

template<class P_V, class P_M>

BaseVectorSequence<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_chain

private

Definition at line 215 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

const BaseEnvironment& QUESO::StatisticalInverseProblem< P_V, P_M >::m_env

private

Definition at line 201 of file StatisticalInverseProblem.h.

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem().

template<class P_V, class P_M>

const BaseScalarFunction<P_V,P_M>& QUESO::StatisticalInverseProblem< P_V, P_M >::m_likelihoodFunction

private

Definition at line 204 of file StatisticalInverseProblem.h.

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem().

template<class P_V, class P_M>

ScalarSequence<double>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_logLikelihoodValues

private

Definition at line 216 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

ScalarSequence<double>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_logTargetValues

private

Definition at line 217 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

MetropolisHastingsSG<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_mhSeqGenerator

private

Definition at line 213 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

MLSampling<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_mlSampler

private

Definition at line 214 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

StatisticalInverseProblemOptions* QUESO::StatisticalInverseProblem< P_V, P_M >::m_optionsObj

private

Definition at line 220 of file StatisticalInverseProblem.h.

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem().

template<class P_V, class P_M>

GenericVectorRV<P_V,P_M>& QUESO::StatisticalInverseProblem< P_V, P_M >::m_postRv

private

Definition at line 205 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

const BaseVectorRV<P_V,P_M>& QUESO::StatisticalInverseProblem< P_V, P_M >::m_priorRv

private

Definition at line 203 of file StatisticalInverseProblem.h.

Referenced by QUESO::StatisticalInverseProblem< P_V, P_M >::StatisticalInverseProblem().

template<class P_V, class P_M>

bool QUESO::StatisticalInverseProblem< P_V, P_M >::m_seedWithMAPEstimator

private

Definition at line 222 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

VectorSet<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_solutionDomain

private

Definition at line 207 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

BaseJointPdf<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_solutionPdf

private

Definition at line 208 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

BaseVectorRealizer<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_solutionRealizer

private

Definition at line 211 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

BaseVectorCdf<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_subSolutionCdf

private

Definition at line 210 of file StatisticalInverseProblem.h.

template<class P_V, class P_M>

BaseVectorMdf<P_V,P_M>* QUESO::StatisticalInverseProblem< P_V, P_M >::m_subSolutionMdf

private

Definition at line 209 of file StatisticalInverseProblem.h.

---

The documentation for this class was generated from the following files:

• src/stats/inc/StatisticalInverseProblem.h
• src/stats/src/StatisticalInverseProblem.C