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...
const BaseVectorSequence< P_V, P_M > &	chain () const
	Returns the MCMC chain; access to private attribute m_chain. More...
const ScalarSequence< double > &	logLikelihoodValues () const
	Returns log likelihood values; access to private attribute m_logLikelihoodValues. More...
const ScalarSequence< double > &	logTargetValues () const
	Returns log target values; access to private attribute m_logTargetValues. 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
const SipOptionsValues *	m_optionsObj
bool	m_seedWithMAPEstimator
bool	m_userDidNotProvideOptions

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 = GslVector, class P_M = GslMatrix>
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 84 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 37 of file StatisticalInverseProblem.C.

References queso_require_equal_to_msg.

  :
  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_optionsObj              (alternativeOptionsValues),
  m_seedWithMAPEstimator    (false),
  m_userDidNotProvideOptions(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 NULL, we create one
  if (m_optionsObj == NULL) {
    SipOptionsValues * tempOptions = new SipOptionsValues(&m_env, prefix);

    // We did this dance because scanOptionsValues is not a const method, but
    // m_optionsObj is a pointer to const
    m_optionsObj = tempOptions;

    // We set this flag so we don't delete the user-created object when it
    // comes time to deconstruct
    m_userDidNotProvideOptions = true;
  }

  if (m_optionsObj->m_help != "") {
    if (m_env.subDisplayFile()) {
      *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
    }
  }

#ifdef QUESO_MEMORY_DEBUGGING
  std::cout << "In Sip, finished scanning options" << std::endl;
#endif

  queso_require_equal_to_msg(priorRv.imageSet().vectorSpace().dimLocal(), likelihoodFunction.domainSet().vectorSpace().dimLocal(), "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");

  queso_require_equal_to_msg(priorRv.imageSet().vectorSpace().dimLocal(), postRv.imageSet().vectorSpace().dimLocal(), "'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 110 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_env, QUESO::SipOptionsValues::m_help, QUESO::StatisticalInverseProblem< P_V, P_M >::m_likelihoodFunction, QUESO::StatisticalInverseProblem< P_V, P_M >::m_optionsObj, QUESO::SipOptionsValues::m_prefix, QUESO::StatisticalInverseProblem< P_V, P_M >::m_priorRv, QUESO::StatisticalInverseProblem< P_V, P_M >::m_userDidNotProvideOptions, QUESO::BaseEnvironment::optionsInputFileName(), queso_require_equal_to_msg, QUESO::BaseEnvironment::subDisplayFile(), and QUESO::VectorSet< V, M >::vectorSpace().

  :
  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_optionsObj              (alternativeOptionsValues),
  m_seedWithMAPEstimator    (false),
  m_userDidNotProvideOptions(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 NULL, we create one
  if (m_optionsObj == NULL) {
    SipOptionsValues * tempOptions = new SipOptionsValues(&m_env, prefix);

    // We did this dance because scanOptionsValues is not a const method, but
    // m_optionsObj is a pointer to const
    m_optionsObj = tempOptions;

    m_userDidNotProvideOptions = true;
  }

  if (m_optionsObj->m_help != "") {
    if (m_env.subDisplayFile()) {
      *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
    }
  }

  queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), m_likelihoodFunction.domainSet().vectorSpace().dimLocal(), "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");

  queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), postRv.imageSet().vectorSpace().dimLocal(), "'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 174 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 && m_userDidNotProvideOptions) {
    delete m_optionsObj;
  }
}

Member Function Documentation

template<class P_V , class P_M >

const BaseVectorSequence< P_V, P_M > & QUESO::StatisticalInverseProblem< P_V, P_M >::chain

(

)

const

Returns the MCMC chain; access to private attribute m_chain.

Only valid after solve has been called.

Definition at line 475 of file StatisticalInverseProblem.C.

References queso_require_msg.

{
  queso_require_msg(m_chain, "m_chain is NULL");
  return *m_chain;
}

template<class P_V = GslVector, class P_M = GslMatrix>

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 512 of file StatisticalInverseProblem.C.

References queso_require_msg.

{
  queso_require_msg(m_mlSampler, "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 498 of file StatisticalInverseProblem.C.

References queso_require_msg.

{
  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
  return m_mlSampler->logEvidence();
}

template<class P_V , class P_M >

const ScalarSequence< double > & QUESO::StatisticalInverseProblem< P_V, P_M >::logLikelihoodValues

(

)

const

Returns log likelihood values; access to private attribute m_logLikelihoodValues.

Only valid for MH and only after solve has been called.

Definition at line 483 of file StatisticalInverseProblem.C.

References queso_require_msg.

{
  queso_require_msg(m_logLikelihoodValues, "m_logLikelihoodValues is NULL");
  return *m_logLikelihoodValues;
}

template<class P_V , class P_M >

const ScalarSequence< double > & QUESO::StatisticalInverseProblem< P_V, P_M >::logTargetValues

(

)

const

Returns log target values; access to private attribute m_logTargetValues.

Only valid for MH and only after solve has been called.

Definition at line 491 of file StatisticalInverseProblem.C.

References queso_require_msg.

{
  queso_require_msg(m_logTargetValues, "m_logTargetValues is NULL");
  return *m_logTargetValues;
}

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 505 of file StatisticalInverseProblem.C.

References queso_require_msg.

{
  queso_require_msg(m_mlSampler, "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 468 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 520 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 461 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 375 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 453 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 203 of file StatisticalInverseProblem.C.

References QUESO::InstantiateIntersection(), QUESO::GslOptimizer::minimize(), QUESO::GslOptimizer::minimizer(), queso_require_equal_to_msg, queso_require_msg, QUESO::OptimizerMonitor::set_display_output(), QUESO::GslOptimizer::setInitialPoint(), 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_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;
  }

  queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialValues.sizeLocal(), "'m_priorRv' and 'initialValues' should have equal dimensions");

  if (initialProposalCovMatrix) {
    queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialProposalCovMatrix->numRowsLocal(), "'m_priorRv' and 'initialProposalCovMatrix' should have equal dimensions");
    queso_require_equal_to_msg(initialProposalCovMatrix->numCols(), initialProposalCovMatrix->numRowsGlobal(), "'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 ||
      m_optionsObj->m_seedWithMAPEstimator) {
    // Unfortunately, I didn't have the foresight to make this an input file
    // option from the beginning, hence needing to check two flags

    // Read in optimizer options from input file
    OptimizerOptions optimizer_options(&m_env, "ip_");

    // Do optimisation before sampling
    GslOptimizer optimizer(optimizer_options, *m_solutionPdf);
    optimizer.setInitialPoint(dynamic_cast<const GslVector &>(initialValues));

    OptimizerMonitor monitor(m_env);
    monitor.set_display_output(true, true);

    // If the input file option is set, then use the monitor, otherwise don't
    if (m_optionsObj->m_useOptimizerMonitor) {
      optimizer.minimize(&monitor);
    }
    else {
      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_dataOutputFileName                       != UQ_SIP_FILENAME_FOR_NO_FILE                    ) &&
      (m_optionsObj->m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->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_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_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_dataOutputFileName+"_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
      }
      queso_require_msg((ofsvar && ofsvar->is_open()), "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_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 382 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_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 = GslVector, class P_M = GslMatrix>

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

friend

Definition at line 209 of file StatisticalInverseProblem.h.

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

Member Data Documentation

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 230 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 216 of file StatisticalInverseProblem.h.

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

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 219 of file StatisticalInverseProblem.h.

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

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 231 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 232 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 228 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 229 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

const SipOptionsValues* QUESO::StatisticalInverseProblem< P_V, P_M >::m_optionsObj

private

Definition at line 234 of file StatisticalInverseProblem.h.

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

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 220 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 218 of file StatisticalInverseProblem.h.

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

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 236 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 222 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 223 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 226 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 225 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 224 of file StatisticalInverseProblem.h.

template<class P_V = GslVector, class P_M = GslMatrix>

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

private

Definition at line 243 of file StatisticalInverseProblem.h.

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

---

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

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