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>
31 #include <queso/BayesianJointPdf.h>
36 template <
class P_V,
class P_M>
44 m_env (priorRv.env()),
46 m_likelihoodFunction (likelihoodFunction),
48 m_solutionDomain (NULL),
50 m_subSolutionMdf (NULL),
51 m_subSolutionCdf (NULL),
52 m_solutionRealizer (NULL),
53 m_mhSeqGenerator (NULL),
56 m_logLikelihoodValues (NULL),
57 m_logTargetValues (NULL),
58 m_optionsObj (alternativeOptionsValues),
59 m_seedWithMAPEstimator (
false),
60 m_userDidNotProvideOptions(
false)
62 #ifdef QUESO_MEMORY_DEBUGGING
63 std::cout <<
"Entering Sip" << std::endl;
65 if (m_env.subDisplayFile()) {
66 *m_env.subDisplayFile() <<
"Entering StatisticalInverseProblem<P_V,P_M>::constructor()"
67 <<
": prefix = " << prefix
68 <<
", alternativeOptionsValues = " << alternativeOptionsValues
69 <<
", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
74 if (m_optionsObj == NULL) {
79 m_optionsObj = tempOptions;
83 m_userDidNotProvideOptions =
true;
86 if (m_optionsObj->m_help !=
"") {
87 if (m_env.subDisplayFile()) {
88 *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
92 #ifdef QUESO_MEMORY_DEBUGGING
93 std::cout <<
"In Sip, finished scanning options" << std::endl;
96 queso_require_equal_to_msg(priorRv.imageSet().vectorSpace().dimLocal(), likelihoodFunction.domainSet().vectorSpace().dimLocal(),
"'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");
98 queso_require_equal_to_msg(priorRv.imageSet().vectorSpace().dimLocal(), postRv.imageSet().vectorSpace().dimLocal(),
"'priorRv' and 'postRv' are related to vector spaces of different dimensions");
100 if (m_env.subDisplayFile()) {
101 *m_env.subDisplayFile() <<
"Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
102 <<
": prefix = " << m_optionsObj->m_prefix
109 template <
class P_V,
class P_M>
116 m_env (gpmsaFactory.m_totalPrior->env()),
117 m_priorRv (*(gpmsaFactory.m_totalPrior)),
118 m_likelihoodFunction (gpmsaFactory.getGPMSAEmulator()),
120 m_solutionDomain (NULL),
121 m_solutionPdf (NULL),
122 m_subSolutionMdf (NULL),
123 m_subSolutionCdf (NULL),
124 m_solutionRealizer (NULL),
125 m_mhSeqGenerator (NULL),
128 m_logLikelihoodValues (NULL),
129 m_logTargetValues (NULL),
130 m_optionsObj (alternativeOptionsValues),
131 m_seedWithMAPEstimator (false),
132 m_userDidNotProvideOptions(false)
136 <<
": prefix = " << prefix
137 <<
", alternativeOptionsValues = " << alternativeOptionsValues
173 template <
class P_V,
class P_M>
180 if (m_logLikelihoodValues) {
181 m_logLikelihoodValues->clear();
182 delete m_logLikelihoodValues;
184 if (m_logTargetValues) {
185 m_logTargetValues->clear();
186 delete m_logTargetValues;
188 if (m_mlSampler )
delete m_mlSampler;
189 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
190 if (m_solutionRealizer)
delete m_solutionRealizer;
191 if (m_subSolutionCdf )
delete m_subSolutionCdf;
192 if (m_subSolutionMdf )
delete m_subSolutionMdf;
193 if (m_solutionPdf )
delete m_solutionPdf;
194 if (m_solutionDomain )
delete m_solutionDomain;
196 if (m_optionsObj && m_userDidNotProvideOptions) {
201 template <
class P_V,
class P_M>
205 const P_V& initialValues,
206 const P_M* initialProposalCovMatrix)
210 m_env.fullComm().Barrier();
212 m_env.fullComm().syncPrintDebugMsg(
"Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
214 if (m_optionsObj->m_computeSolution ==
false) {
215 if ((m_env.subDisplayFile())) {
216 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
217 <<
": avoiding solution, as requested by user"
222 if ((m_env.subDisplayFile())) {
223 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
224 <<
": computing solution, as requested by user"
228 queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialValues.sizeLocal(),
"'m_priorRv' and 'initialValues' should have equal dimensions");
230 if (initialProposalCovMatrix) {
231 queso_require_equal_to_msg(m_priorRv.imageSet().vectorSpace().dimLocal(), initialProposalCovMatrix->numRowsLocal(),
"'m_priorRv' and 'initialProposalCovMatrix' should have equal dimensions");
232 queso_require_equal_to_msg(initialProposalCovMatrix->numCols(), initialProposalCovMatrix->numRowsGlobal(),
"'initialProposalCovMatrix' should be a square matrix");
235 if (m_mlSampler )
delete m_mlSampler;
236 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
237 if (m_solutionRealizer)
delete m_solutionRealizer;
238 if (m_subSolutionCdf )
delete m_subSolutionCdf;
239 if (m_subSolutionMdf )
delete m_subSolutionMdf;
240 if (m_solutionPdf )
delete m_solutionPdf;
241 if (m_solutionDomain )
delete m_solutionDomain;
243 P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
244 numEvaluationPointsVec.cwSet(250.);
251 m_likelihoodFunction,
255 m_postRv.setPdf(*m_solutionPdf);
262 if (this->m_seedWithMAPEstimator ||
263 m_optionsObj->m_seedWithMAPEstimator) {
271 GslOptimizer optimizer(optimizer_options, *m_solutionPdf);
272 optimizer.
setInitialPoint(dynamic_cast<const GslVector &>(initialValues));
278 if (m_optionsObj->m_useOptimizerMonitor) {
287 m_optionsObj->m_prefix.c_str(), alternativeOptionsValues,
288 m_postRv, optimizer.
minimizer(), initialProposalCovMatrix);
293 m_optionsObj->m_prefix.c_str(), alternativeOptionsValues, m_postRv,
294 initialValues, initialProposalCovMatrix);
299 m_optionsObj->m_prefix +
303 m_optionsObj->m_prefix +
307 m_mhSeqGenerator->generateSequence(*m_chain, m_logLikelihoodValues,
313 m_postRv.setRealizer(*m_solutionRealizer);
315 m_env.fullComm().syncPrintDebugMsg(
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings(), code place 1",3,3000000);
319 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
322 m_chain->subUniformlySampledMdf(numEvaluationPointsVec,
328 m_postRv.setMdf(*m_subSolutionMdf);
331 (m_optionsObj->m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_dataOutputAllowedSet.end())) {
332 if (m_env.subRank() == 0) {
334 if (m_env.subDisplayFile()) {
335 *m_env.subDisplayFile() <<
"Opening data output file '" << m_optionsObj->m_dataOutputFileName
336 <<
"' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
342 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);
343 if ((ofsvar == NULL ) ||
344 (ofsvar->is_open() ==
false)) {
346 ofsvar =
new std::ofstream((m_optionsObj->m_dataOutputFileName+
"_sub"+m_env.subIdString()+
".m").c_str(), std::ofstream::out | std::ofstream::trunc);
350 m_postRv.mdf().print(*ofsvar);
355 if (m_env.subDisplayFile()) {
356 *m_env.subDisplayFile() <<
"Closed data output file '" << m_optionsObj->m_dataOutputFileName
357 <<
"' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
363 if (m_env.subDisplayFile()) {
364 *m_env.subDisplayFile() << std::endl;
367 m_env.fullComm().syncPrintDebugMsg(
"Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
368 m_env.fullComm().Barrier();
373 template <
class P_V,
class P_M>
377 this->m_seedWithMAPEstimator =
true;
380 template <
class P_V,
class P_M>
384 m_env.fullComm().Barrier();
385 m_env.fullComm().syncPrintDebugMsg(
"Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
387 if (m_optionsObj->m_computeSolution ==
false) {
388 if ((m_env.subDisplayFile())) {
389 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
390 <<
": avoiding solution, as requested by user"
395 if ((m_env.subDisplayFile())) {
396 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
397 <<
": computing solution, as requested by user"
401 if (m_mlSampler )
delete m_mlSampler;
402 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
403 if (m_solutionRealizer)
delete m_solutionRealizer;
404 if (m_subSolutionCdf )
delete m_subSolutionCdf;
405 if (m_subSolutionMdf )
delete m_subSolutionMdf;
406 if (m_solutionPdf )
delete m_solutionPdf;
407 if (m_solutionDomain )
delete m_solutionDomain;
409 P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
410 numEvaluationPointsVec.cwSet(250.);
417 m_likelihoodFunction,
421 m_postRv.setPdf(*m_solutionPdf);
428 m_likelihoodFunction);
432 m_mlSampler->generateSequence(*m_chain,
439 m_postRv.setRealizer(*m_solutionRealizer);
441 if (m_env.subDisplayFile()) {
442 *m_env.subDisplayFile() << std::endl;
445 m_env.fullComm().syncPrintDebugMsg(
"Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
446 m_env.fullComm().Barrier();
451 template <
class P_V,
class P_M>
455 return *m_mhSeqGenerator;
459 template <
class P_V,
class P_M>
466 template <
class P_V,
class P_M>
473 template <
class P_V,
class P_M>
481 template <
class P_V,
class P_M>
486 return *m_logLikelihoodValues;
489 template <
class P_V,
class P_M>
494 return *m_logTargetValues;
497 template <
class P_V,
class P_M>
501 return m_mlSampler->logEvidence();
504 template <
class P_V,
class P_M>
508 return m_mlSampler->meanLogLikelihood();
511 template <
class P_V,
class P_M>
515 return m_mlSampler->eig();
518 template <
class P_V,
class P_M>
A class for handling sequential draws (sampling) from probability density distributions.
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
Class to accommodate arrays of one-dimensional grid.
This class provides options for a Statistical Inverse Problem if no input file is available...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
A class for handling sampled vector MDFs.
~StatisticalInverseProblem()
Destructor.
void setInitialPoint(const GslVector &intialPoint)
Set the point at which the optimization starts.
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
const GenericVectorRV< P_V, P_M > & postRv() const
Returns the Posterior RV; access to private attribute m_postrRv.
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
double logEvidence() const
Returns the logarithm value of the evidence. Related to ML.
void solveWithBayesMLSampling()
Solves with Bayes Multi-Level (ML) sampling.
#define queso_require_equal_to_msg(expr1, expr2, msg)
bool m_userDidNotProvideOptions
A templated class that represents a Metropolis-Hastings generator of samples.
double meanLogLikelihood() const
Returns the mean of the logarithm value of the likelihood. Related to ML.
virtual void minimize(OptimizerMonitor *monitor=NULL)
Minimize the objective function, starting at m_initialPoint.
#define UQ_SIP_FILENAME_FOR_NO_FILE
This class provides options for a Optimizer.
void set_display_output(bool enable_output, bool print_xmin)
Enabling output to std::cout everytime append is called.
const SipOptionsValues * m_optionsObj
void seedWithMAPEstimator()
Seeds the chain with the result of a deterministic optimisation.
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
const GslVector & minimizer() const
Return the point that minimizes the objective function.
This templated class represents a Statistical Inverse Problem.
A base class for handling optimisation of scalar functions.
void solveWithBayesMetropolisHastings(const MhOptionsValues *alternativeOptionsValues, const P_V &initialValues, const P_M *initialProposalCovMatrix)
Solves the problem via Bayes formula and a Metropolis-Hastings algorithm.
const BaseVectorRV< P_V, P_M > & priorRv() const
Returns the Prior RV; access to private attribute m_priorRv.
const BaseVectorRV< P_V, P_M > & m_priorRv
#define queso_require_msg(asserted, msg)
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 BaseVectorSequence< P_V, P_M > & chain() const
Returns the MCMC chain; access to private attribute m_chain.
Object to monitor convergence of optimizers.
This class provides options for the Metropolis-Hastings generator of samples if no input file is avai...
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
const ScalarSequence< double > & logTargetValues() const
Returns log target values; access to private attribute m_logTargetValues.
Class to accommodate arrays of one-dimensional tables.
A templated class that represents a Multilevel generator of samples.
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
std::string m_help
If this string is non-empty, options are print to the output file.
const BaseEnvironment & m_env
void print(std::ostream &os) const
TODO: Prints the sequence.
unsigned int dimLocal() const
const MetropolisHastingsSG< P_V, P_M > & sequenceGenerator() const
Return the underlying MetropolisHastingSG object.
const ScalarSequence< double > & logLikelihoodValues() const
Returns log likelihood values; access to private attribute m_logLikelihoodValues. ...