25 #include <queso/StatisticalInverseProblem.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 #include <queso/GPMSA.h>
29 #include <queso/GslOptimizer.h>
30 #include <queso/OptimizerMonitor.h>
35 template <
class P_V,
class P_M>
43 m_env (priorRv.env()),
45 m_likelihoodFunction (likelihoodFunction),
47 m_solutionDomain (NULL),
49 m_subSolutionMdf (NULL),
50 m_subSolutionCdf (NULL),
51 m_solutionRealizer (NULL),
52 m_mhSeqGenerator (NULL),
55 m_logLikelihoodValues (NULL),
56 m_logTargetValues (NULL),
57 m_optionsObj (alternativeOptionsValues),
58 m_seedWithMAPEstimator (
false),
59 m_userDidNotProvideOptions(
false)
61 #ifdef QUESO_MEMORY_DEBUGGING
62 std::cout <<
"Entering Sip" << std::endl;
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()
73 if (m_optionsObj == NULL) {
78 m_optionsObj = tempOptions;
82 m_userDidNotProvideOptions =
true;
85 if (m_optionsObj->m_help !=
"") {
86 if (m_env.subDisplayFile()) {
87 *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
91 #ifdef QUESO_MEMORY_DEBUGGING
92 std::cout <<
"In Sip, finished scanning options" << std::endl;
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");
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");
99 if (m_env.subDisplayFile()) {
100 *m_env.subDisplayFile() <<
"Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
101 <<
": prefix = " << m_optionsObj->m_prefix
108 template <
class P_V,
class P_M>
115 m_env (gpmsaFactory.m_totalPrior->env()),
116 m_priorRv (*(gpmsaFactory.m_totalPrior)),
117 m_likelihoodFunction (gpmsaFactory.getGPMSAEmulator()),
119 m_solutionDomain (NULL),
120 m_solutionPdf (NULL),
121 m_subSolutionMdf (NULL),
122 m_subSolutionCdf (NULL),
123 m_solutionRealizer (NULL),
124 m_mhSeqGenerator (NULL),
127 m_logLikelihoodValues (NULL),
128 m_logTargetValues (NULL),
129 m_optionsObj (alternativeOptionsValues),
130 m_seedWithMAPEstimator (false),
131 m_userDidNotProvideOptions(false)
135 <<
": prefix = " << prefix
136 <<
", alternativeOptionsValues = " << alternativeOptionsValues
172 template <
class P_V,
class P_M>
179 if (m_logLikelihoodValues) {
180 m_logLikelihoodValues->clear();
181 delete m_logLikelihoodValues;
183 if (m_logTargetValues) {
184 m_logTargetValues->clear();
185 delete m_logTargetValues;
187 if (m_mlSampler )
delete m_mlSampler;
188 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
189 if (m_solutionRealizer)
delete m_solutionRealizer;
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;
195 if (m_optionsObj && m_userDidNotProvideOptions) {
200 template <
class P_V,
class P_M>
204 const P_V& initialValues,
205 const P_M* initialProposalCovMatrix)
209 m_env.fullComm().Barrier();
211 m_env.fullComm().syncPrintDebugMsg(
"Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
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"
221 if ((m_env.subDisplayFile())) {
222 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
223 <<
": computing solution, as requested by user"
227 queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialValues.sizeLocal(),
"'m_priorRv' and 'initialValues' should have equal dimensions");
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");
234 if (m_mlSampler )
delete m_mlSampler;
235 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
236 if (m_solutionRealizer)
delete m_solutionRealizer;
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;
242 P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
243 numEvaluationPointsVec.cwSet(250.);
250 m_likelihoodFunction,
254 m_postRv.setPdf(*m_solutionPdf);
261 if (this->m_seedWithMAPEstimator ||
262 m_optionsObj->m_seedWithMAPEstimator) {
268 optimizer.
setInitialPoint(dynamic_cast<const GslVector &>(initialValues));
274 if (m_optionsObj->m_useOptimizerMonitor) {
283 m_optionsObj->m_prefix.c_str(), alternativeOptionsValues,
284 m_postRv, optimizer.
minimizer(), initialProposalCovMatrix);
289 m_optionsObj->m_prefix.c_str(), alternativeOptionsValues, m_postRv,
290 initialValues, initialProposalCovMatrix);
295 m_optionsObj->m_prefix +
299 m_optionsObj->m_prefix +
303 m_mhSeqGenerator->generateSequence(*m_chain, m_logLikelihoodValues,
309 m_postRv.setRealizer(*m_solutionRealizer);
311 m_env.fullComm().syncPrintDebugMsg(
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings(), code place 1",3,3000000);
315 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
318 m_chain->subUniformlySampledMdf(numEvaluationPointsVec,
324 m_postRv.setMdf(*m_subSolutionMdf);
327 (m_optionsObj->m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_dataOutputAllowedSet.end())) {
328 if (m_env.subRank() == 0) {
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
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)) {
342 ofsvar =
new std::ofstream((m_optionsObj->m_dataOutputFileName+
"_sub"+m_env.subIdString()+
".m").c_str(), std::ofstream::out | std::ofstream::trunc);
346 m_postRv.mdf().print(*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
359 if (m_env.subDisplayFile()) {
360 *m_env.subDisplayFile() << std::endl;
363 m_env.fullComm().syncPrintDebugMsg(
"Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
364 m_env.fullComm().Barrier();
369 template <
class P_V,
class P_M>
373 this->m_seedWithMAPEstimator =
true;
376 template <
class P_V,
class P_M>
380 m_env.fullComm().Barrier();
381 m_env.fullComm().syncPrintDebugMsg(
"Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
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"
391 if ((m_env.subDisplayFile())) {
392 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
393 <<
": computing solution, as requested by user"
397 if (m_mlSampler )
delete m_mlSampler;
398 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
399 if (m_solutionRealizer)
delete m_solutionRealizer;
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;
405 P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
406 numEvaluationPointsVec.cwSet(250.);
413 m_likelihoodFunction,
417 m_postRv.setPdf(*m_solutionPdf);
424 m_likelihoodFunction);
428 m_mlSampler->generateSequence(*m_chain,
435 m_postRv.setRealizer(*m_solutionRealizer);
437 if (m_env.subDisplayFile()) {
438 *m_env.subDisplayFile() << std::endl;
441 m_env.fullComm().syncPrintDebugMsg(
"Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
442 m_env.fullComm().Barrier();
447 template <
class P_V,
class P_M>
451 return *m_mhSeqGenerator;
455 template <
class P_V,
class P_M>
462 template <
class P_V,
class P_M>
469 template <
class P_V,
class P_M>
477 template <
class P_V,
class P_M>
482 return *m_logLikelihoodValues;
485 template <
class P_V,
class P_M>
490 return *m_logTargetValues;
493 template <
class P_V,
class P_M>
497 return m_mlSampler->logEvidence();
500 template <
class P_V,
class P_M>
504 return m_mlSampler->meanLogLikelihood();
507 template <
class P_V,
class P_M>
511 return m_mlSampler->eig();
514 template <
class P_V,
class P_M>
Class to accommodate arrays of one-dimensional grid.
void solveWithBayesMetropolisHastings(const MhOptionsValues *alternativeOptionsValues, const P_V &initialValues, const P_M *initialProposalCovMatrix)
Solves the problem via Bayes formula and a Metropolis-Hastings algorithm.
This templated class represents a Statistical Inverse Problem.
void seedWithMAPEstimator()
Seeds the chain with the result of a deterministic optimisation.
A class for handling sampled vector MDFs.
This class provides options for the Metropolis-Hastings generator of samples if no input file is avai...
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.
const SipOptionsValues * m_optionsObj
A base class for handling optimisation of scalar functions.
const BaseVectorRV< P_V, P_M > & m_priorRv
void setInitialPoint(const GslVector &intialPoint)
Set the point at which the optimization starts.
const BaseVectorSequence< P_V, P_M > & chain() const
Returns the MCMC chain; access to private attribute m_chain.
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
const ScalarSequence< double > & logTargetValues() const
Returns log target values; access to private attribute m_logTargetValues.
const BaseEnvironment & m_env
const BaseVectorRV< P_V, P_M > & priorRv() const
Returns the Prior RV; access to private attribute m_priorRv.
const ScalarSequence< double > & logLikelihoodValues() const
Returns log likelihood values; access to private attribute m_logLikelihoodValues. ...
virtual void minimize(OptimizerMonitor *monitor=NULL)
Minimize the objective function, starting at m_initialPoint.
#define queso_require_msg(asserted, msg)
#define queso_require_equal_to_msg(expr1, expr2, msg)
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
A class for handling sequential draws (sampling) from probability density distributions.
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).
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
const MetropolisHastingsSG< P_V, P_M > & sequenceGenerator() const
Return the underlying MetropolisHastingSG object.
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
std::string m_help
If this string is non-empty, options are print to the output file.
This class provides options for a Statistical Inverse Problem if no input file is available...
void print(std::ostream &os) const
TODO: Prints the sequence.
A templated class that represents a Metropolis-Hastings generator of samples.
Class to accommodate arrays of one-dimensional tables.
const GslVector & minimizer() const
Return the point that minimizes the objective function.
~StatisticalInverseProblem()
Destructor.
void set_display_output(bool enable_output, bool print_xmin)
Enabling output to std::cout everytime append is called.
bool m_userDidNotProvideOptions
const GenericVectorRV< P_V, P_M > & postRv() const
Returns the Posterior RV; access to private attribute m_postrRv.
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
void solveWithBayesMLSampling()
Solves with Bayes Multi-Level (ML) sampling.
double meanLogLikelihood() const
Returns the mean of the logarithm value of the likelihood. Related to ML.
#define UQ_SIP_FILENAME_FOR_NO_FILE
double logEvidence() const
Returns the logarithm value of the evidence. Related to ML.
unsigned int dimLocal() const
Object to monitor convergence of optimizers.
A templated class that represents a Multilevel generator of samples.