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>
461 "StatisticalInverseProblem<P_V,P_M>::logEvidence()",
462 "m_mlSampler is NULL");
463 return m_mlSampler->logEvidence();
466 template <
class P_V,
class P_M>
471 "StatisticalInverseProblem<P_V,P_M>::meanLogLikelihood()",
472 "m_mlSampler is NULL");
473 return m_mlSampler->meanLogLikelihood();
476 template <
class P_V,
class P_M>
481 "StatisticalInverseProblem<P_V,P_M>::eig()",
482 "m_mlSampler is NULL");
483 return m_mlSampler->eig();
486 template <
class P_V,
class P_M>
unsigned int dimLocal() const
void solveWithBayesMetropolisHastings(const MhOptionsValues *alternativeOptionsValues, const P_V &initialValues, const P_M *initialProposalCovMatrix)
Solves the problem via Bayes formula and a Metropolis-Hastings algorithm.
Class to accommodate arrays of one-dimensional grid.
A class for handling sampled vector MDFs.
void scanOptionsValues()
It scans the option values from the options input file.
double logEvidence() const
Returns the logarithm value of the evidence. Related to ML.
int worldRank() const
Returns the process world rank.
const BaseVectorRV< P_V, P_M > & m_priorRv
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
void seedWithMAPEstimator()
Seeds the chain with the result of a deterministic optimisation.
const BaseEnvironment & m_env
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 VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
A templated class that represents a Multilevel generator of samples.
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
double meanLogLikelihood() const
Returns the mean of the logarithm value of the likelihood. Related to ML.
A templated class that represents a Metropolis-Hastings generator of samples.
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
#define UQ_SIP_FILENAME_FOR_NO_FILE
const BaseVectorRV< P_V, P_M > & priorRv() const
Returns the Prior RV; access to private attribute m_priorRv.
virtual void minimize(OptimizerMonitor *monitor=NULL)
Minimize the objective function, starting at m_initialPoint.
A class for handling sequential draws (sampling) from probability density distributions.
const BaseScalarFunction< P_V, P_M > & m_likelihoodFunction
void solveWithBayesMLSampling()
Solves with Bayes Multi-Level (ML) sampling.
void setInitialPoint(const GslVector &intialPoint)
Set the point at which the optimization starts.
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
This class reads option values for a Statistical Inverse Problem from an input file.
Class to accommodate arrays of one-dimensional tables.
const GslVector & minimizer() const
Return the point that minimizes the objective function.
A base class for handling optimisation of scalar functions.
void print(std::ostream &os) const
TODO: Prints the sequence.
Class for handling vector samples (sequence of vectors).
StatisticalInverseProblemOptions * m_optionsObj
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
SipOptionsValues m_alternativeOptionsValues
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
~StatisticalInverseProblem()
Destructor.
This templated class represents a Statistical Inverse Problem.
This class provides options for the Metropolis-Hastings generator of samples if no input file is avai...
const MetropolisHastingsSG< P_V, P_M > & sequenceGenerator() const
Return the underlying MetropolisHastingSG object.
This class provides options for a Statistical Inverse Problem if no input file is available...
const GenericVectorRV< P_V, P_M > & postRv() const
Returns the Posterior RV; access to private attribute m_postrRv.