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	solveWithBayesMetropolisHastings (const MhOptionsValues *alternativeOptionsValues, const P_V &initialValues)
	Solves the problem via Bayes formula and a Metropolis-Hastings algorithm. More...
void	solveWithBayesMetropolisHastings (const MhOptionsValues *alternativeOptionsValues)
	Solves the problem via Bayes formula and a Metropolis-Hastings algorithm. More...
void	solveWithBayesMetropolisHastings ()
	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
ScopedPtr< VectorSet< P_V, P_M > >::Type	m_solutionDomain
ScopedPtr< BaseJointPdf< P_V, P_M > >::Type	m_solutionPdf
ScopedPtr< BaseVectorMdf< P_V, P_M > >::Type	m_subSolutionMdf
ScopedPtr< BaseVectorCdf< P_V, P_M > >::Type	m_subSolutionCdf
ScopedPtr< BaseVectorRealizer < P_V, P_M > >::Type	m_solutionRealizer
ScopedPtr < MetropolisHastingsSG< P_V, P_M > >::Type	m_mhSeqGenerator
ScopedPtr< MLSampling< P_V, P_M > >::Type	m_mlSampler
ScopedPtr< BaseVectorSequence < P_V, P_M > >::Type	m_chain
ScopedPtr< ScalarSequence < double > >::Type	m_logLikelihoodValues
ScopedPtr< ScalarSequence < double > >::Type	m_logTargetValues
ScopedPtr< const SipOptionsValues >::Type	m_optionsObj
bool	m_seedWithMAPEstimator
ScopedPtr< ArrayOfOneDGrids < P_V, P_M > >	m_subMdfGrids
ScopedPtr< ArrayOfOneDTables < P_V, P_M > >	m_subMdfValues

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 85 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 alternativeOptionsValues is NULL and an input file is specified, the constructor reads input options that begin with the string '<prefix>ip_'.

If alternativeOptionsValues is not NULL, the input file is ignored and construction copies the object pointed to by alternativeOptionsValues to and stores the copy internally. Users may delete the object poined to by alternativeOptionsValues. Users cannot change the options object after StatisticalInverseProblem has been constructed.

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::queso_require_equal_to_msg, and QUESO::BaseEnvironment::subDisplayFile().

  :
  m_env                     (priorRv.env()),
  m_priorRv                 (priorRv),
  m_likelihoodFunction      (likelihoodFunction),
  m_postRv                  (postRv),
  m_solutionDomain          (),
  m_solutionPdf             (),
  m_subSolutionMdf          (),
  m_subSolutionCdf          (),
  m_solutionRealizer        (),
  m_mhSeqGenerator          (),
  m_mlSampler               (),
  m_chain                   (),
  m_logLikelihoodValues     (),
  m_logTargetValues         (),
  m_optionsObj              (),
  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 NULL, we create one
  if (m_optionsObj == NULL) {
    m_optionsObj.reset(new SipOptionsValues(&m_env, prefix));
  }
  else {
    // If the user passed in a legitimate options object, then copy it into a
    // ScopedPtr
    m_optionsObj.reset(new SipOptionsValues(*alternativeOptionsValues));
  }

  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 alternativeOptionsValues is NULL and an input file is specified, the constructor reads input options that begin with the string '<prefix>ip_'.

If alternativeOptionsValues is not NULL, the input file is ignored and construction copies the object pointed to by alternativeOptionsValues to and stores the copy internally. Users may delete the object poined to by alternativeOptionsValues. Users cannot change the options object after StatisticalInverseProblem has been constructed.

Definition at line 106 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::StatisticalInverseProblem< P_V, P_M >::m_likelihoodFunction, QUESO::StatisticalInverseProblem< P_V, P_M >::m_optionsObj, QUESO::StatisticalInverseProblem< P_V, P_M >::m_priorRv, QUESO::BaseEnvironment::optionsInputFileName(), QUESO::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          (),
  m_solutionPdf             (),
  m_subSolutionMdf          (),
  m_subSolutionCdf          (),
  m_solutionRealizer        (),
  m_mhSeqGenerator          (),
  m_mlSampler               (),
  m_chain                   (),
  m_logLikelihoodValues     (),
  m_logTargetValues         (),
  m_optionsObj              (),
  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 NULL, we create one
  if (m_optionsObj == NULL) {
    m_optionsObj.reset(new SipOptionsValues(&m_env, prefix));
  }
  else {
    // If the user passed in a legitimate options object, then copy it into a
    // ScopedPtr
    m_optionsObj.reset(new SipOptionsValues(*alternativeOptionsValues));
  }


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

{
}

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

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

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

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

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

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

{
  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 458 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 510 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 451 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 373 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 443 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 175 of file StatisticalInverseProblem.C.

References QUESO::InstantiateIntersection(), QUESO::GslOptimizer::minimize(), QUESO::GslOptimizer::minimizer(), QUESO::queso_require_equal_to_msg, QUESO::OptimizerMonitor::set_display_output(), and QUESO::GslOptimizer::setInitialPoint().

{
  //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");
  }

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

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

  m_solutionPdf.reset(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.reset(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.reset(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.reset(new MetropolisHastingsSG<P_V, P_M>(
        m_optionsObj->m_prefix.c_str(), alternativeOptionsValues, m_postRv,
        initialValues, initialProposalCovMatrix));
  }


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

  m_logTargetValues.reset(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.get(),
                                     m_logTargetValues.get());

  m_solutionRealizer.reset(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.reset(new ArrayOfOneDGrids <P_V,P_M>((m_optionsObj->m_prefix+"Mdf_").c_str(),m_postRv.imageSet().vectorSpace()));
  m_subMdfValues.reset(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.reset(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
}

template<class P_V , class P_M >

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

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

This calls solveWithBayesMetropolisHastings(const MhOptionsValues *, const P_V &, P_M *); with the prior variance as the proposal covariance matrix in the third argument.

Definition at line 338 of file StatisticalInverseProblem.C.

{
  P_M proposalCovMatrix(m_priorRv.imageSet().vectorSpace().zeroVector());

  m_priorRv.pdf().distributionVariance(proposalCovMatrix);

  this->solveWithBayesMetropolisHastings(alternativeOptionsValues,
                                         initialValues,
                                         &proposalCovMatrix);
}

template<class P_V , class P_M >

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

(

const MhOptionsValues *

alternativeOptionsValues

)

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

This calls solveWithBayesMetropolisHastings(const MhOptionsValues *, const P_V &); with the prior mean as the initial point in the second argument.

Definition at line 353 of file StatisticalInverseProblem.C.

{
  P_V initialValues(m_priorRv.imageSet().vectorSpace().zeroVector());

  m_priorRv.pdf().distributionMean(initialValues);

  this->solveWithBayesMetropolisHastings(alternativeOptionsValues,
                                         initialValues);
}

template<class P_V , class P_M >

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

(

)

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

This calls solveWithBayesMetropolisHastings(const MhOptionsValues *); with NULL as the options object. Options are read from the input file.

Definition at line 366 of file StatisticalInverseProblem.C.

{
  this->solveWithBayesMetropolisHastings(NULL);
}

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 380 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;
  }

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

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

  m_solutionPdf.reset(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.reset(new SequenceOfVectors<P_V,P_M>(m_postRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"chain"));
  m_mlSampler.reset(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.reset(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 251 of file StatisticalInverseProblem.h.

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

Member Data Documentation

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

ScopedPtr<BaseVectorSequence <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_chain

private

Definition at line 272 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 258 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 261 of file StatisticalInverseProblem.h.

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

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

ScopedPtr<ScalarSequence<double> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_logLikelihoodValues

private

Definition at line 273 of file StatisticalInverseProblem.h.

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

ScopedPtr<ScalarSequence<double> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_logTargetValues

private

Definition at line 274 of file StatisticalInverseProblem.h.

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

ScopedPtr<MetropolisHastingsSG<P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_mhSeqGenerator

private

Definition at line 270 of file StatisticalInverseProblem.h.

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

ScopedPtr<MLSampling <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_mlSampler

private

Definition at line 271 of file StatisticalInverseProblem.h.

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

ScopedPtr<const SipOptionsValues>::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_optionsObj

private

Definition at line 276 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 262 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 260 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 278 of file StatisticalInverseProblem.h.

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

ScopedPtr<VectorSet <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_solutionDomain

private

Definition at line 264 of file StatisticalInverseProblem.h.

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

ScopedPtr<BaseJointPdf <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_solutionPdf

private

Definition at line 265 of file StatisticalInverseProblem.h.

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

ScopedPtr<BaseVectorRealizer <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_solutionRealizer

private

Definition at line 268 of file StatisticalInverseProblem.h.

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

ScopedPtr<ArrayOfOneDGrids <P_V,P_M> > QUESO::StatisticalInverseProblem< P_V, P_M >::m_subMdfGrids

private

Definition at line 281 of file StatisticalInverseProblem.h.

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

ScopedPtr<ArrayOfOneDTables <P_V,P_M> > QUESO::StatisticalInverseProblem< P_V, P_M >::m_subMdfValues

private

Definition at line 282 of file StatisticalInverseProblem.h.

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

ScopedPtr<BaseVectorCdf <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_subSolutionCdf

private

Definition at line 267 of file StatisticalInverseProblem.h.

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

ScopedPtr<BaseVectorMdf <P_V,P_M> >::Type QUESO::StatisticalInverseProblem< P_V, P_M >::m_subSolutionMdf

private

Definition at line 266 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