queso-0.53.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 35 of file StatisticalInverseProblem.C.

References queso_require_equal_to_msg.

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

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

Destructor.

Definition at line 172 of file StatisticalInverseProblem.C.

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

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

References queso_require_msg.

457 {
458  queso_require_msg(m_chain, "m_chain is NULL");
459  return *m_chain;
460 }
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
BaseVectorSequence< P_V, P_M > * 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 493 of file StatisticalInverseProblem.C.

References queso_require_msg.

494 {
495  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
496  return m_mlSampler->eig();
497 }
#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 479 of file StatisticalInverseProblem.C.

References queso_require_msg.

480 {
481  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
482  return m_mlSampler->logEvidence();
483 }
#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 464 of file StatisticalInverseProblem.C.

References queso_require_msg.

465 {
466  queso_require_msg(m_logLikelihoodValues, "m_logLikelihoodValues is NULL");
467  return *m_logLikelihoodValues;
468 }
#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 472 of file StatisticalInverseProblem.C.

References queso_require_msg.

473 {
474  queso_require_msg(m_logTargetValues, "m_logTargetValues is NULL");
475  return *m_logTargetValues;
476 }
ScalarSequence< double > * m_logTargetValues
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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 486 of file StatisticalInverseProblem.C.

References queso_require_msg.

487 {
488  queso_require_msg(m_mlSampler, "m_mlSampler is NULL");
489  return m_mlSampler->meanLogLikelihood();
490 }
#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 449 of file StatisticalInverseProblem.C.

450 {
451  return m_postRv;
452 }
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 501 of file StatisticalInverseProblem.C.

502 {
503  return;
504 }
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 442 of file StatisticalInverseProblem.C.

443 {
444  return m_priorRv;
445 }
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 356 of file StatisticalInverseProblem.C.

357 {
358  this->m_seedWithMAPEstimator = true;
359 }
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 434 of file StatisticalInverseProblem.C.

435 {
436  return *m_mhSeqGenerator;
437 }
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 201 of file StatisticalInverseProblem.C.

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

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

References QUESO::InstantiateIntersection().

364 {
365  m_env.fullComm().Barrier();
366  m_env.fullComm().syncPrintDebugMsg("Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
367 
368  if (m_optionsObj->m_computeSolution == false) {
369  if ((m_env.subDisplayFile())) {
370  *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
371  << ": avoiding solution, as requested by user"
372  << std::endl;
373  }
374  return;
375  }
376  if ((m_env.subDisplayFile())) {
377  *m_env.subDisplayFile() << "In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
378  << ": computing solution, as requested by user"
379  << std::endl;
380  }
381 
382  if (m_mlSampler ) delete m_mlSampler;
383  if (m_mhSeqGenerator ) delete m_mhSeqGenerator;
385  if (m_subSolutionCdf ) delete m_subSolutionCdf;
386  if (m_subSolutionMdf ) delete m_subSolutionMdf;
387  if (m_solutionPdf ) delete m_solutionPdf;
388  if (m_solutionDomain ) delete m_solutionDomain;
389 
390  P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
391  numEvaluationPointsVec.cwSet(250.);
392 
393  // Compute output pdf up to a multiplicative constant: Bayesian approach
395 
396  m_solutionPdf = new BayesianJointPdf<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
397  m_priorRv.pdf(),
399  1.,
401 
402  m_postRv.setPdf(*m_solutionPdf);
403 
404  // Compute output realizer: ML approach
405  m_chain = new SequenceOfVectors<P_V,P_M>(m_postRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"chain");
406  m_mlSampler = new MLSampling<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
407  //m_postRv,
408  m_priorRv,
410  // initialValues,
411  // initialProposalCovMatrix);
412 
413  m_mlSampler->generateSequence(*m_chain,
414  NULL,
415  NULL);
416 
417  m_solutionRealizer = new SequentialVectorRealizer<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
418  *m_chain);
419 
421 
422  if (m_env.subDisplayFile()) {
423  *m_env.subDisplayFile() << std::endl;
424  }
425 
426  m_env.fullComm().syncPrintDebugMsg("Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
427  m_env.fullComm().Barrier();
428 
429  return;
430 }
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
void Barrier() const
Pause every process in *this communicator until all the processes reach this point.
Definition: MpiComm.C:122
const BaseVectorRV< P_V, P_M > & m_priorRv
const MpiComm & fullComm() const
Access function for MpiComm full communicator.
Definition: Environment.C:228
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:189
BaseJointPdf< P_V, P_M > * m_solutionPdf
BaseVectorSequence< P_V, P_M > * m_chain
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:274
BaseVectorMdf< P_V, P_M > * m_subSolutionMdf
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:86
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
MetropolisHastingsSG< P_V, P_M > * m_mhSeqGenerator
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:79
void setRealizer(BaseVectorRealizer< V, M > &realizer)
Sets the realizer of this vector RV to realizer.
BaseVectorCdf< P_V, P_M > * m_subSolutionCdf
GenericVectorRV< P_V, P_M > & m_postRv
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
void syncPrintDebugMsg(const char *msg, unsigned int msgVerbosity, unsigned int numUSecs) const
Synchronizes all the processes and print debug message.
Definition: MpiComm.C:195
BaseVectorRealizer< P_V, P_M > * m_solutionRealizer

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 Thu Jun 11 2015 13:52:36 for queso-0.53.0 by  doxygen 1.8.5