queso-0.55.0
Private Attributes | List of all members
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 >:
Collaboration graph
[legend]

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 BaseEnvironmentm_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 SipOptionsValuesm_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 $ \pi()$ denote the mathematical likelihood function and $ x $ denote a vector of parameters. The likelihood function object stores the routine that computes $ \pi(x) $ and whatever data necessary by such routine. The routine in the likelihood function object can compute either the actual value $ \pi(x) $ or the value $ \ln[\pi(x)] $. 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
prefixThe prefix
alternativeOptionsValuesOptions (if no input file)
priorRvThe prior RV
likelihoodFunctionThe likelihood function
postRvThe posterior RV

Definition at line 36 of file StatisticalInverseProblem.C.

References queso_require_equal_to_msg.

42  :
43  m_env (priorRv.env()),
44  m_priorRv (priorRv),
45  m_likelihoodFunction (likelihoodFunction),
46  m_postRv (postRv),
47  m_solutionDomain (NULL),
48  m_solutionPdf (NULL),
49  m_subSolutionMdf (NULL),
50  m_subSolutionCdf (NULL),
51  m_solutionRealizer (NULL),
52  m_mhSeqGenerator (NULL),
53  m_mlSampler (NULL),
54  m_chain (NULL),
55  m_logLikelihoodValues (NULL),
56  m_logTargetValues (NULL),
57  m_optionsObj (alternativeOptionsValues),
58  m_seedWithMAPEstimator (false),
60 {
61 #ifdef QUESO_MEMORY_DEBUGGING
62  std::cout << "Entering Sip" << std::endl;
63 #endif
64  if (m_env.subDisplayFile()) {
65  *m_env.subDisplayFile() << "Entering StatisticalInverseProblem<P_V,P_M>::constructor()"
66  << ": prefix = " << prefix
67  << ", alternativeOptionsValues = " << alternativeOptionsValues
68  << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
69  << std::endl;
70  }
71 
72  // If NULL, we create one
73  if (m_optionsObj == NULL) {
74  SipOptionsValues * tempOptions = new SipOptionsValues(&m_env, prefix);
75 
76  // We did this dance because scanOptionsValues is not a const method, but
77  // m_optionsObj is a pointer to const
78  m_optionsObj = tempOptions;
79 
80  // We set this flag so we don't delete the user-created object when it
81  // comes time to deconstruct
83  }
84 
85  if (m_optionsObj->m_help != "") {
86  if (m_env.subDisplayFile()) {
87  *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
88  }
89  }
90 
91 #ifdef QUESO_MEMORY_DEBUGGING
92  std::cout << "In Sip, finished scanning options" << std::endl;
93 #endif
94 
95  queso_require_equal_to_msg(priorRv.imageSet().vectorSpace().dimLocal(), likelihoodFunction.domainSet().vectorSpace().dimLocal(), "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");
96 
97  queso_require_equal_to_msg(priorRv.imageSet().vectorSpace().dimLocal(), postRv.imageSet().vectorSpace().dimLocal(), "'priorRv' and 'postRv' are related to vector spaces of different dimensions");
98 
99  if (m_env.subDisplayFile()) {
100  *m_env.subDisplayFile() << "Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
101  << ": prefix = " << m_optionsObj->m_prefix
102  << std::endl;
103  }
104 
105  return;
106 }
BaseVectorSequence< P_V, P_M > * m_chain
BaseVectorCdf< P_V, P_M > * m_subSolutionCdf
const BaseVectorRV< P_V, P_M > & m_priorRv
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:75
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
BaseJointPdf< P_V, P_M > * m_solutionPdf
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
MetropolisHastingsSG< P_V, P_M > * m_mhSeqGenerator
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:311
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
std::string m_help
If this string is non-empty, options are print to the output file.
ScalarSequence< double > * m_logLikelihoodValues
BaseVectorRealizer< P_V, P_M > * m_solutionRealizer
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:82
BaseVectorMdf< P_V, P_M > * m_subSolutionMdf
ScalarSequence< double > * m_logTargetValues
unsigned int dimLocal() const
Definition: VectorSpace.C:170
GenericVectorRV< P_V, P_M > & m_postRv
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 109 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().

114  :
115  m_env (gpmsaFactory.m_totalPrior->env()),
116  m_priorRv (*(gpmsaFactory.m_totalPrior)),
117  m_likelihoodFunction (gpmsaFactory.getGPMSAEmulator()),
118  m_postRv (postRv),
119  m_solutionDomain (NULL),
120  m_solutionPdf (NULL),
121  m_subSolutionMdf (NULL),
122  m_subSolutionCdf (NULL),
123  m_solutionRealizer (NULL),
124  m_mhSeqGenerator (NULL),
125  m_mlSampler (NULL),
126  m_chain (NULL),
127  m_logLikelihoodValues (NULL),
128  m_logTargetValues (NULL),
129  m_optionsObj (alternativeOptionsValues),
130  m_seedWithMAPEstimator (false),
132 {
133  if (m_env.subDisplayFile()) {
134  *m_env.subDisplayFile() << "Entering StatisticalInverseProblem<P_V,P_M>::constructor()"
135  << ": prefix = " << prefix
136  << ", alternativeOptionsValues = " << alternativeOptionsValues
137  << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
138  << std::endl;
139  }
140 
141  // If NULL, we create one
142  if (m_optionsObj == NULL) {
143  SipOptionsValues * tempOptions = new SipOptionsValues(&m_env, prefix);
144 
145  // We did this dance because scanOptionsValues is not a const method, but
146  // m_optionsObj is a pointer to const
147  m_optionsObj = tempOptions;
148 
150  }
151 
152  if (m_optionsObj->m_help != "") {
153  if (m_env.subDisplayFile()) {
154  *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
155  }
156  }
157 
158  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");
159 
160  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");
161 
162  if (m_env.subDisplayFile()) {
163  *m_env.subDisplayFile() << "Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
164  << ": prefix = " << m_optionsObj->m_prefix
165  << std::endl;
166  }
167 
168  return;
169 }
BaseVectorSequence< P_V, P_M > * m_chain
BaseVectorCdf< P_V, P_M > * m_subSolutionCdf
const BaseVectorRV< P_V, P_M > & m_priorRv
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
BaseJointPdf< P_V, P_M > * m_solutionPdf
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
MetropolisHastingsSG< P_V, P_M > * m_mhSeqGenerator
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:311
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
std::string m_help
If this string is non-empty, options are print to the output file.
ScalarSequence< double > * m_logLikelihoodValues
BaseVectorRealizer< P_V, P_M > * m_solutionRealizer
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:82
BaseVectorMdf< P_V, P_M > * m_subSolutionMdf
ScalarSequence< double > * m_logTargetValues
unsigned int dimLocal() const
Definition: VectorSpace.C:170
GenericVectorRV< P_V, P_M > & m_postRv
template<class P_V , class P_M >
QUESO::StatisticalInverseProblem< P_V, P_M >::~StatisticalInverseProblem ( )

Destructor.

Definition at line 173 of file StatisticalInverseProblem.C.

174 {
175  if (m_chain) {
176  m_chain->clear();
177  delete m_chain;
178  }
179  if (m_logLikelihoodValues) {
181  delete m_logLikelihoodValues;
182  }
183  if (m_logTargetValues) {
185  delete m_logTargetValues;
186  }
187  if (m_mlSampler ) delete m_mlSampler;
188  if (m_mhSeqGenerator ) delete m_mhSeqGenerator;
190  if (m_subSolutionCdf ) delete m_subSolutionCdf;
191  if (m_subSolutionMdf ) delete m_subSolutionMdf;
192  if (m_solutionPdf ) delete m_solutionPdf;
193  if (m_solutionDomain ) delete m_solutionDomain;
194 
196  delete m_optionsObj;
197  }
198 }
BaseVectorSequence< P_V, P_M > * m_chain
BaseVectorCdf< P_V, P_M > * m_subSolutionCdf
void clear()
Clears the sequence of scalars.
BaseJointPdf< P_V, P_M > * m_solutionPdf
void clear()
Reset the values and the size of the sequence of vectors.
MetropolisHastingsSG< P_V, P_M > * m_mhSeqGenerator
ScalarSequence< double > * m_logLikelihoodValues
BaseVectorRealizer< P_V, P_M > * m_solutionRealizer
BaseVectorMdf< P_V, P_M > * m_subSolutionMdf
ScalarSequence< double > * m_logTargetValues

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

References queso_require_msg.

472 {
473  queso_require_msg(m_chain, "m_chain is NULL");
474  return *m_chain;
475 }
BaseVectorSequence< P_V, P_M > * m_chain
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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 508 of file StatisticalInverseProblem.C.

References queso_require_msg.

509 {
510  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
511  return m_mlSampler->eig();
512 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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 494 of file StatisticalInverseProblem.C.

References queso_require_msg.

495 {
496  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
497  return m_mlSampler->logEvidence();
498 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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 479 of file StatisticalInverseProblem.C.

References queso_require_msg.

480 {
481  queso_require_msg(m_logLikelihoodValues, "m_logLikelihoodValues is NULL");
482  return *m_logLikelihoodValues;
483 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
ScalarSequence< double > * 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 487 of file StatisticalInverseProblem.C.

References queso_require_msg.

488 {
489  queso_require_msg(m_logTargetValues, "m_logTargetValues is NULL");
490  return *m_logTargetValues;
491 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
ScalarSequence< double > * 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 501 of file StatisticalInverseProblem.C.

References queso_require_msg.

502 {
503  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
504  return m_mlSampler->meanLogLikelihood();
505 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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 464 of file StatisticalInverseProblem.C.

465 {
466  return m_postRv;
467 }
GenericVectorRV< P_V, P_M > & 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 516 of file StatisticalInverseProblem.C.

517 {
518  return;
519 }
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 457 of file StatisticalInverseProblem.C.

458 {
459  return m_priorRv;
460 }
const BaseVectorRV< P_V, P_M > & 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 371 of file StatisticalInverseProblem.C.

372 {
373  this->m_seedWithMAPEstimator = true;
374 }
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 449 of file StatisticalInverseProblem.C.

450 {
451  return *m_mhSeqGenerator;
452 }
MetropolisHastingsSG< P_V, P_M > * 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 202 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.

206 {
207  //grvy_timer_begin("BayesMetropolisHastings"); TODO: revisit timing output
208  //std::cout << "proc " << m_env.fullRank() << ", HERE sip 000" << std::endl;
209  m_env.fullComm().Barrier();
210  //std::cout << "proc " << m_env.fullRank() << ", HERE sip 001" << std::endl;
211  m_env.fullComm().syncPrintDebugMsg("Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
212 
213  if (m_optionsObj->m_computeSolution == false) {
214  if ((m_env.subDisplayFile())) {
215  *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
216  << ": avoiding solution, as requested by user"
217  << std::endl;
218  }
219  return;
220  }
221  if ((m_env.subDisplayFile())) {
222  *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
223  << ": computing solution, as requested by user"
224  << std::endl;
225  }
226 
227  queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialValues.sizeLocal(), "'m_priorRv' and 'initialValues' should have equal dimensions");
228 
229  if (initialProposalCovMatrix) {
230  queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialProposalCovMatrix->numRowsLocal(), "'m_priorRv' and 'initialProposalCovMatrix' should have equal dimensions");
231  queso_require_equal_to_msg(initialProposalCovMatrix->numCols(), initialProposalCovMatrix->numRowsGlobal(), "'initialProposalCovMatrix' should be a square matrix");
232  }
233 
234  if (m_mlSampler ) delete m_mlSampler;
235  if (m_mhSeqGenerator ) delete m_mhSeqGenerator;
237  if (m_subSolutionCdf ) delete m_subSolutionCdf;
238  if (m_subSolutionMdf ) delete m_subSolutionMdf;
239  if (m_solutionPdf ) delete m_solutionPdf;
240  if (m_solutionDomain ) delete m_solutionDomain;
241 
242  P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
243  numEvaluationPointsVec.cwSet(250.);
244 
245  // Compute output pdf up to a multiplicative constant: Bayesian approach
247 
248  m_solutionPdf = new BayesianJointPdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
249  m_priorRv.pdf(),
251  1.,
253 
254  m_postRv.setPdf(*m_solutionPdf);
255  m_chain = new SequenceOfVectors<P_V,P_M>(m_postRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"chain");
256 
257  // Decide whether or not to create a MetropolisHastingsSG instance from the
258  // user-provided initial seed, or use the user-provided seed for a
259  // deterministic optimisation instead and seed the chain with the result of
260  // the optimisation
261  if (this->m_seedWithMAPEstimator ||
263  // Unfortunately, I didn't have the foresight to make this an input file
264  // option from the beginning, hence needing to check two flags
265 
266  // Do optimisation before sampling
267  GslOptimizer optimizer(*m_solutionPdf);
268  optimizer.setInitialPoint(dynamic_cast<const GslVector &>(initialValues));
269 
270  OptimizerMonitor monitor(m_env);
271  monitor.set_display_output(true, true);
272 
273  // If the input file option is set, then use the monitor, otherwise don't
275  optimizer.minimize(&monitor);
276  }
277  else {
278  optimizer.minimize();
279  }
280 
281  // Compute output realizer: Metropolis-Hastings approach
282  m_mhSeqGenerator = new MetropolisHastingsSG<P_V, P_M>(
283  m_optionsObj->m_prefix.c_str(), alternativeOptionsValues,
284  m_postRv, optimizer.minimizer(), initialProposalCovMatrix);
285  }
286  else {
287  // Compute output realizer: Metropolis-Hastings approach
288  m_mhSeqGenerator = new MetropolisHastingsSG<P_V, P_M>(
289  m_optionsObj->m_prefix.c_str(), alternativeOptionsValues, m_postRv,
290  initialValues, initialProposalCovMatrix);
291  }
292 
293 
294  m_logLikelihoodValues = new ScalarSequence<double>(m_env, 0,
296  "logLike");
297 
298  m_logTargetValues = new ScalarSequence<double>(m_env, 0,
300  "logTarget");
301 
302  // m_logLikelihoodValues and m_logTargetValues may be NULL
303  m_mhSeqGenerator->generateSequence(*m_chain, m_logLikelihoodValues,
304  m_logTargetValues);
305 
306  m_solutionRealizer = new SequentialVectorRealizer<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
307  *m_chain);
308 
309  m_postRv.setRealizer(*m_solutionRealizer);
310 
311  m_env.fullComm().syncPrintDebugMsg("In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings(), code place 1",3,3000000);
312  //m_env.fullComm().Barrier();
313 
314  // Compute output mdf: uniform sampling approach
315 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
316  m_subMdfGrids = new ArrayOfOneDGrids <P_V,P_M>((m_optionsObj->m_prefix+"Mdf_").c_str(),m_postRv.imageSet().vectorSpace());
317  m_subMdfValues = new ArrayOfOneDTables<P_V,P_M>((m_optionsObj->m_prefix+"Mdf_").c_str(),m_postRv.imageSet().vectorSpace());
318  m_chain->subUniformlySampledMdf(numEvaluationPointsVec, // input
319  *m_subMdfGrids, // output
320  *m_subMdfValues); // output
321  m_subSolutionMdf = new SampledVectorMdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
322  *m_subMdfGrids,
323  *m_subMdfValues);
324  m_postRv.setMdf(*m_subSolutionMdf);
325 
327  (m_optionsObj->m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_dataOutputAllowedSet.end())) {
328  if (m_env.subRank() == 0) {
329  // Write data output file
330  if (m_env.subDisplayFile()) {
331  *m_env.subDisplayFile() << "Opening data output file '" << m_optionsObj->m_dataOutputFileName
332  << "' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
333  << std::endl;
334  }
335 
336  // Open file
337  // Always write at the end of an eventual pre-existing file
338  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);
339  if ((ofsvar == NULL ) ||
340  (ofsvar->is_open() == false)) {
341  delete ofsvar;
342  ofsvar = new std::ofstream((m_optionsObj->m_dataOutputFileName+"_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
343  }
344  queso_require_msg((ofsvar && ofsvar->is_open()), "failed to open file");
345 
346  m_postRv.mdf().print(*ofsvar);
347 
348  // Close file
349  //ofsvar->close();
350  delete ofsvar;
351  if (m_env.subDisplayFile()) {
352  *m_env.subDisplayFile() << "Closed data output file '" << m_optionsObj->m_dataOutputFileName
353  << "' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
354  << std::endl;
355  }
356  }
357  }
358 #endif
359  if (m_env.subDisplayFile()) {
360  *m_env.subDisplayFile() << std::endl;
361  }
362 
363  m_env.fullComm().syncPrintDebugMsg("Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
364  m_env.fullComm().Barrier();
365  // grvy_timer_end("BayesMetropolisHastings"); TODO: revisit timers
366  return;
367 }
BaseVectorSequence< P_V, P_M > * m_chain
void syncPrintDebugMsg(const char *msg, unsigned int msgVerbosity, unsigned int numUSecs) const
Synchronizes all the processes and print debug message.
Definition: MpiComm.C:323
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
BaseVectorCdf< P_V, P_M > * m_subSolutionCdf
const BaseVectorRV< P_V, P_M > & m_priorRv
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
const MpiComm & fullComm() const
Access function for MpiComm full communicator.
Definition: Environment.C:232
BaseJointPdf< P_V, P_M > * m_solutionPdf
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:164
MetropolisHastingsSG< P_V, P_M > * m_mhSeqGenerator
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:89
ScalarSequence< double > * m_logLikelihoodValues
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:189
BaseVectorRealizer< P_V, P_M > * m_solutionRealizer
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:82
#define UQ_SIP_FILENAME_FOR_NO_FILE
std::set< unsigned int > m_dataOutputAllowedSet
BaseVectorMdf< P_V, P_M > * m_subSolutionMdf
ScalarSequence< double > * m_logTargetValues
unsigned int dimLocal() const
Definition: VectorSpace.C:170
GenericVectorRV< P_V, P_M > & m_postRv
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 378 of file StatisticalInverseProblem.C.

References QUESO::InstantiateIntersection().

379 {
380  m_env.fullComm().Barrier();
381  m_env.fullComm().syncPrintDebugMsg("Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
382 
383  if (m_optionsObj->m_computeSolution == false) {
384  if ((m_env.subDisplayFile())) {
385  *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
386  << ": avoiding solution, as requested by user"
387  << std::endl;
388  }
389  return;
390  }
391  if ((m_env.subDisplayFile())) {
392  *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
393  << ": computing solution, as requested by user"
394  << std::endl;
395  }
396 
397  if (m_mlSampler ) delete m_mlSampler;
398  if (m_mhSeqGenerator ) delete m_mhSeqGenerator;
400  if (m_subSolutionCdf ) delete m_subSolutionCdf;
401  if (m_subSolutionMdf ) delete m_subSolutionMdf;
402  if (m_solutionPdf ) delete m_solutionPdf;
403  if (m_solutionDomain ) delete m_solutionDomain;
404 
405  P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
406  numEvaluationPointsVec.cwSet(250.);
407 
408  // Compute output pdf up to a multiplicative constant: Bayesian approach
410 
411  m_solutionPdf = new BayesianJointPdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
412  m_priorRv.pdf(),
414  1.,
416 
417  m_postRv.setPdf(*m_solutionPdf);
418 
419  // Compute output realizer: ML approach
420  m_chain = new SequenceOfVectors<P_V,P_M>(m_postRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"chain");
421  m_mlSampler = new MLSampling<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
422  //m_postRv,
423  m_priorRv,
425  // initialValues,
426  // initialProposalCovMatrix);
427 
428  m_mlSampler->generateSequence(*m_chain,
429  NULL,
430  NULL);
431 
432  m_solutionRealizer = new SequentialVectorRealizer<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
433  *m_chain);
434 
436 
437  if (m_env.subDisplayFile()) {
438  *m_env.subDisplayFile() << std::endl;
439  }
440 
441  m_env.fullComm().syncPrintDebugMsg("Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
442  m_env.fullComm().Barrier();
443 
444  return;
445 }
BaseVectorSequence< P_V, P_M > * m_chain
void syncPrintDebugMsg(const char *msg, unsigned int msgVerbosity, unsigned int numUSecs) const
Synchronizes all the processes and print debug message.
Definition: MpiComm.C:323
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
BaseVectorCdf< P_V, P_M > * m_subSolutionCdf
const BaseVectorRV< P_V, P_M > & m_priorRv
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
void setRealizer(BaseVectorRealizer< V, M > &realizer)
Sets the realizer of this vector RV to realizer.
const MpiComm & fullComm() const
Access function for MpiComm full communicator.
Definition: Environment.C:232
BaseJointPdf< P_V, P_M > * m_solutionPdf
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:164
MetropolisHastingsSG< P_V, P_M > * m_mhSeqGenerator
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:89
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:189
BaseVectorRealizer< P_V, P_M > * m_solutionRealizer
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:82
BaseVectorMdf< P_V, P_M > * m_subSolutionMdf
GenericVectorRV< P_V, P_M > & m_postRv

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.

210  {
211  obj.print(os);
212  return os;
213  }

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
template<class P_V = GslVector, class P_M = GslMatrix>
const BaseScalarFunction<P_V,P_M>& QUESO::StatisticalInverseProblem< P_V, P_M >::m_likelihoodFunction
private
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
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
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

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

Generated on Fri Jun 17 2016 14:17:46 for queso-0.55.0 by  doxygen 1.8.5