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>
34 template <
class P_V,
class P_M>
42 m_env (priorRv.env()),
44 m_likelihoodFunction (likelihoodFunction),
46 m_solutionDomain (NULL),
48 m_subSolutionMdf (NULL),
49 m_subSolutionCdf (NULL),
50 m_solutionRealizer (NULL),
51 m_mhSeqGenerator (NULL),
54 m_logLikelihoodValues (NULL),
55 m_logTargetValues (NULL),
56 m_alternativeOptionsValues(),
58 m_seedWithMAPEstimator (
false)
60 #ifdef QUESO_MEMORY_DEBUGGING
61 std::cout <<
"Entering Sip" << std::endl;
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()
71 if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
72 if (m_env.optionsInputFileName() ==
"") {
77 m_optionsObj->scanOptionsValues();
79 #ifdef QUESO_MEMORY_DEBUGGING
80 std::cout <<
"In Sip, finished scanning options" << std::endl;
83 UQ_FATAL_TEST_MACRO(priorRv.imageSet().vectorSpace().dimLocal() != likelihoodFunction.domainSet().vectorSpace().dimLocal(),
85 "StatisticalInverseProblem<P_V,P_M>::constructor()",
86 "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");
88 UQ_FATAL_TEST_MACRO(priorRv.imageSet().vectorSpace().dimLocal() != postRv.imageSet().vectorSpace().dimLocal(),
90 "StatisticalInverseProblem<P_V,P_M>::constructor()",
91 "'priorRv' and 'postRv' are related to vector spaces of different dimensions");
93 if (m_env.subDisplayFile()) {
94 *m_env.subDisplayFile() <<
"Leaving StatisticalInverseProblem<P_V,P_M>::constructor()"
95 <<
": prefix = " << m_optionsObj->m_prefix
102 template <
class P_V,
class P_M>
106 const GPMSAFactory<P_V, P_M> & gpmsaFactory,
109 m_env (gpmsaFactory.m_totalPrior->env()),
110 m_priorRv (*(gpmsaFactory.m_totalPrior)),
111 m_likelihoodFunction (gpmsaFactory.getGPMSAEmulator()),
113 m_solutionDomain (NULL),
114 m_solutionPdf (NULL),
115 m_subSolutionMdf (NULL),
116 m_subSolutionCdf (NULL),
117 m_solutionRealizer (NULL),
118 m_mhSeqGenerator (NULL),
121 m_logLikelihoodValues (NULL),
122 m_logTargetValues (NULL),
123 m_alternativeOptionsValues(),
125 m_seedWithMAPEstimator (false)
129 <<
": prefix = " << prefix
130 <<
", alternativeOptionsValues = " << alternativeOptionsValues
146 "StatisticalInverseProblem<P_V,P_M>::constructor()",
147 "'priorRv' and 'likelihoodFunction' are related to vector spaces of different dimensions");
151 "StatisticalInverseProblem<P_V,P_M>::constructor()",
152 "'priorRv' and 'postRv' are related to vector spaces of different dimensions");
164 template <
class P_V,
class P_M>
171 if (m_logLikelihoodValues) {
172 m_logLikelihoodValues->clear();
173 delete m_logLikelihoodValues;
175 if (m_logTargetValues) {
176 m_logTargetValues->clear();
177 delete m_logTargetValues;
179 if (m_mlSampler )
delete m_mlSampler;
180 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
181 if (m_solutionRealizer)
delete m_solutionRealizer;
182 if (m_subSolutionCdf )
delete m_subSolutionCdf;
183 if (m_subSolutionMdf )
delete m_subSolutionMdf;
184 if (m_solutionPdf )
delete m_solutionPdf;
185 if (m_solutionDomain )
delete m_solutionDomain;
186 if (m_optionsObj )
delete m_optionsObj;
189 template <
class P_V,
class P_M>
193 const P_V& initialValues,
194 const P_M* initialProposalCovMatrix)
198 m_env.fullComm().Barrier();
200 m_env.fullComm().syncPrintDebugMsg(
"Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
202 if (m_optionsObj->m_ov.m_computeSolution ==
false) {
203 if ((m_env.subDisplayFile())) {
204 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
205 <<
": avoiding solution, as requested by user"
210 if ((m_env.subDisplayFile())) {
211 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()"
212 <<
": computing solution, as requested by user"
216 UQ_FATAL_TEST_MACRO(m_priorRv.imageSet().vectorSpace().dimLocal() != initialValues.sizeLocal(),
218 "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
219 "'m_priorRv' and 'initialValues' should have equal dimensions");
221 if (initialProposalCovMatrix) {
222 UQ_FATAL_TEST_MACRO(m_priorRv.imageSet().vectorSpace().dimLocal() != initialProposalCovMatrix->numRowsLocal(),
224 "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
225 "'m_priorRv' and 'initialProposalCovMatrix' should have equal dimensions");
226 UQ_FATAL_TEST_MACRO(initialProposalCovMatrix->numCols() != initialProposalCovMatrix->numRowsGlobal(),
228 "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
229 "'initialProposalCovMatrix' should be a square matrix");
232 if (m_mlSampler )
delete m_mlSampler;
233 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
234 if (m_solutionRealizer)
delete m_solutionRealizer;
235 if (m_subSolutionCdf )
delete m_subSolutionCdf;
236 if (m_subSolutionMdf )
delete m_subSolutionMdf;
237 if (m_solutionPdf )
delete m_solutionPdf;
238 if (m_solutionDomain )
delete m_solutionDomain;
240 P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
241 numEvaluationPointsVec.cwSet(250.);
248 m_likelihoodFunction,
252 m_postRv.setPdf(*m_solutionPdf);
259 if (this->m_seedWithMAPEstimator) {
262 optimizer.
setInitialPoint(dynamic_cast<const GslVector &>(initialValues));
267 m_optionsObj->m_prefix.c_str(), alternativeOptionsValues,
268 m_postRv, optimizer.
minimizer(), initialProposalCovMatrix);
273 m_optionsObj->m_prefix.c_str(), alternativeOptionsValues, m_postRv,
274 initialValues, initialProposalCovMatrix);
279 m_optionsObj->m_prefix +
283 m_optionsObj->m_prefix +
287 m_mhSeqGenerator->generateSequence(*m_chain, m_logLikelihoodValues,
293 m_postRv.setRealizer(*m_solutionRealizer);
295 m_env.fullComm().syncPrintDebugMsg(
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings(), code place 1",3,3000000);
299 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
302 m_chain->subUniformlySampledMdf(numEvaluationPointsVec,
308 m_postRv.setMdf(*m_subSolutionMdf);
311 (m_optionsObj->m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_ov.m_dataOutputAllowedSet.end())) {
312 if (m_env.subRank() == 0) {
314 if (m_env.subDisplayFile()) {
315 *m_env.subDisplayFile() <<
"Opening data output file '" << m_optionsObj->m_ov.m_dataOutputFileName
316 <<
"' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
322 std::ofstream* ofsvar =
new std::ofstream((m_optionsObj->m_ov.m_dataOutputFileName+
"_sub"+m_env.subIdString()+
".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
323 if ((ofsvar == NULL ) ||
324 (ofsvar->is_open() ==
false)) {
326 ofsvar =
new std::ofstream((m_optionsObj->m_ov.m_dataOutputFileName+
"_sub"+m_env.subIdString()+
".m").c_str(), std::ofstream::out | std::ofstream::trunc);
330 "StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",
331 "failed to open file");
333 m_postRv.mdf().print(*ofsvar);
338 if (m_env.subDisplayFile()) {
339 *m_env.subDisplayFile() <<
"Closed data output file '" << m_optionsObj->m_ov.m_dataOutputFileName
340 <<
"' for calibration problem with problem with prefix = " << m_optionsObj->m_prefix
346 if (m_env.subDisplayFile()) {
347 *m_env.subDisplayFile() << std::endl;
350 m_env.fullComm().syncPrintDebugMsg(
"Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMetropolisHastings()",1,3000000);
351 m_env.fullComm().Barrier();
356 template <
class P_V,
class P_M>
360 this->m_seedWithMAPEstimator =
true;
363 template <
class P_V,
class P_M>
367 m_env.fullComm().Barrier();
368 m_env.fullComm().syncPrintDebugMsg(
"Entering StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
370 if (m_optionsObj->m_ov.m_computeSolution ==
false) {
371 if ((m_env.subDisplayFile())) {
372 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
373 <<
": avoiding solution, as requested by user"
378 if ((m_env.subDisplayFile())) {
379 *m_env.subDisplayFile() <<
"In StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()"
380 <<
": computing solution, as requested by user"
384 if (m_mlSampler )
delete m_mlSampler;
385 if (m_mhSeqGenerator )
delete m_mhSeqGenerator;
386 if (m_solutionRealizer)
delete m_solutionRealizer;
387 if (m_subSolutionCdf )
delete m_subSolutionCdf;
388 if (m_subSolutionMdf )
delete m_subSolutionMdf;
389 if (m_solutionPdf )
delete m_solutionPdf;
390 if (m_solutionDomain )
delete m_solutionDomain;
392 P_V numEvaluationPointsVec(m_priorRv.imageSet().vectorSpace().zeroVector());
393 numEvaluationPointsVec.cwSet(250.);
400 m_likelihoodFunction,
404 m_postRv.setPdf(*m_solutionPdf);
411 m_likelihoodFunction);
415 m_mlSampler->generateSequence(*m_chain,
422 m_postRv.setRealizer(*m_solutionRealizer);
424 if (m_env.subDisplayFile()) {
425 *m_env.subDisplayFile() << std::endl;
428 m_env.fullComm().syncPrintDebugMsg(
"Leaving StatisticalInverseProblem<P_V,P_M>::solveWithBayesMLSampling()",1,3000000);
429 m_env.fullComm().Barrier();
434 template <
class P_V,
class P_M>
438 return *m_mhSeqGenerator;
442 template <
class P_V,
class P_M>
449 template <
class P_V,
class P_M>
456 template <
class P_V,
class P_M>
462 "StatisticalInverseProblem<P_V,P_M>::chain()",
467 template <
class P_V,
class P_M>
473 "StatisticalInverseProblem<P_V,P_M>::logLikelihoodValues()",
474 "m_logLikelihoodValues is NULL");
475 return *m_logLikelihoodValues;
478 template <
class P_V,
class P_M>
484 "StatisticalInverseProblem<P_V,P_M>::logTargetValues()",
485 "m_logTargetValues is NULL");
486 return *m_logTargetValues;
489 template <
class P_V,
class P_M>
494 "StatisticalInverseProblem<P_V,P_M>::logEvidence()",
495 "m_mlSampler is NULL");
496 return m_mlSampler->logEvidence();
499 template <
class P_V,
class P_M>
504 "StatisticalInverseProblem<P_V,P_M>::meanLogLikelihood()",
505 "m_mlSampler is NULL");
506 return m_mlSampler->meanLogLikelihood();
509 template <
class P_V,
class P_M>
514 "StatisticalInverseProblem<P_V,P_M>::eig()",
515 "m_mlSampler is NULL");
516 return m_mlSampler->eig();
519 template <
class P_V,
class P_M>
This class provides options for the Metropolis-Hastings generator of samples if no input file is avai...
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 grid.
StatisticalInverseProblemOptions * m_optionsObj
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
This class provides options for a Statistical Inverse Problem if no input file is available...
SipOptionsValues m_alternativeOptionsValues
double meanLogLikelihood() const
Returns the mean of the logarithm value of the likelihood. Related to ML.
~StatisticalInverseProblem()
Destructor.
int worldRank() const
Returns the process world rank.
This templated class represents a Statistical Inverse Problem.
unsigned int dimLocal() const
#define UQ_SIP_FILENAME_FOR_NO_FILE
A class for handling sampled vector MDFs.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
const MetropolisHastingsSG< P_V, P_M > & sequenceGenerator() const
Return the underlying MetropolisHastingSG object.
const ScalarSequence< double > & logTargetValues() const
Returns log target values; access to private attribute m_logTargetValues.
const GenericVectorRV< P_V, P_M > & postRv() const
Returns the Posterior RV; access to private attribute m_postrRv.
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
double logEvidence() const
Returns the logarithm value of the evidence. Related to ML.
A templated class that represents a Multilevel generator of samples.
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
virtual void minimize(OptimizerMonitor *monitor=NULL)
Minimize the objective function, starting at m_initialPoint.
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
void seedWithMAPEstimator()
Seeds the chain with the result of a deterministic optimisation.
const BaseEnvironment & m_env
const ScalarSequence< double > & logLikelihoodValues() const
Returns log likelihood values; access to private attribute m_logLikelihoodValues. ...
void setInitialPoint(const GslVector &intialPoint)
Set the point at which the optimization starts.
A class for handling sequential draws (sampling) from probability density distributions.
Class to accommodate arrays of one-dimensional tables.
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
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.
This class reads option values for a Statistical Inverse Problem from an input file.
A base class for handling optimisation of scalar functions.
Class for handling vector samples (sequence of vectors).
const GslVector & minimizer() const
Return the point that minimizes the objective function.
const BaseVectorRV< P_V, P_M > & m_priorRv
void solveWithBayesMLSampling()
Solves with Bayes Multi-Level (ML) sampling.
const BaseVectorSequence< P_V, P_M > & chain() const
Returns the MCMC chain; access to private attribute m_chain.
void scanOptionsValues()
It scans the option values from the options input file.
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.