queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
StatisticalForwardProblem.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/StatisticalForwardProblem.h>
26 #include <queso/SequentialVectorRealizer.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 
30 namespace QUESO {
31 
32 // Default constructor -----------------------------
33 template <class P_V,class P_M,class Q_V,class Q_M>
34 StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::StatisticalForwardProblem( const char* prefix, const SfpOptionsValues* alternativeOptionsValues, // dakota const BaseVectorRV <P_V,P_M>& paramRv, const BaseVectorFunction<P_V,P_M,Q_V,Q_M>& qoiFunction, GenericVectorRV <Q_V,Q_M>& qoiRv)
40  :
41  m_env (paramRv.env()),
42  m_paramRv (paramRv),
43  m_qoiFunction (qoiFunction),
44  m_qoiRv (qoiRv),
45  m_paramChain (NULL),
46  m_qoiChain (NULL),
47  m_mcSeqGenerator (NULL),
48  m_solutionRealizer (NULL),
49 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
50  m_subMdfGrids (NULL),
51  m_subMdfValues (NULL),
52 #endif
53 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
54  m_subSolutionMdf (NULL),
55  m_subCdfGrids (NULL),
56  m_subCdfValues (NULL),
57  m_subSolutionCdf (NULL),
58  m_unifiedCdfGrids (NULL),
59  m_unifiedCdfValues (NULL),
60  m_unifiedSolutionCdf (NULL),
61 #endif
62  m_solutionPdf (NULL),
63  m_optionsObj (alternativeOptionsValues),
64  m_userDidNotProvideOptions(false)
65 {
66  if (m_env.subDisplayFile()) {
67  *m_env.subDisplayFile() << "Entering StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::constructor()"
68  << ": prefix = " << prefix
69  << ", alternativeOptionsValues = " << alternativeOptionsValues
70  << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
71  << std::endl;
72  }
73 
74  // If NULL, we create one
75  if (m_optionsObj == NULL) {
76  SfpOptionsValues * tempOptions = new SfpOptionsValues(&m_env, prefix);
77 
78  // We did this dance because scanOptionsValues is not a const method, but
79  // m_optionsObj is a pointer to const
80  m_optionsObj = tempOptions;
81 
82  // We do this so we don't delete the user's object in the dtor
83  m_userDidNotProvideOptions = true;
84  }
85 
86  if (m_optionsObj->m_help != "") {
87  if (m_env.subDisplayFile()) {
88  *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
89  }
90  }
91 
92  queso_require_equal_to_msg(paramRv.imageSet().vectorSpace().dimLocal(), qoiFunction.domainSet().vectorSpace().dimLocal(), "'paramRv' and 'qoiFunction' are related to vector spaces of different dimensions");
93 
94  queso_require_equal_to_msg(qoiFunction.imageSet().vectorSpace().dimLocal(), qoiRv.imageSet().vectorSpace().dimLocal(), "'qoiFunction' and 'qoiRv' are related to vector spaces of different dimensions");
95 
96  if (m_env.subDisplayFile()) {
97  *m_env.subDisplayFile() << "Leaving StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::constructor()"
98  << ": prefix = " << m_optionsObj->m_prefix
99  << std::endl;
100  }
101 }
102 
103 // Destructor --------------------------------------
104 template <class P_V,class P_M,class Q_V,class Q_M>
106 {
107  if (m_solutionPdf ) delete m_solutionPdf;
108 
109 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
110  if (m_unifiedSolutionCdf) delete m_unifiedSolutionCdf;
111  if (m_unifiedCdfValues ) delete m_unifiedCdfValues;
112  if (m_unifiedCdfGrids ) delete m_unifiedCdfGrids;
113 
114  if (m_subSolutionCdf ) delete m_subSolutionCdf;
115  if (m_subCdfValues ) delete m_subCdfValues;
116  if (m_subCdfGrids ) delete m_subCdfGrids;
117 
118  if (m_subSolutionMdf ) delete m_subSolutionMdf;
119 #endif
120 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
121  if (m_subMdfValues ) delete m_subMdfValues;
122  if (m_subMdfGrids ) delete m_subMdfGrids;
123 #endif
124  if (m_solutionRealizer ) delete m_solutionRealizer;
125 
126  if (m_mcSeqGenerator ) delete m_mcSeqGenerator;
127 
128  if (m_qoiChain) {
129  m_qoiChain->clear();
130  delete m_qoiChain;
131  }
132 
133  if (m_paramChain) {
134  m_paramChain->clear();
135  delete m_paramChain;
136  }
137 
138  if (m_optionsObj ) delete m_optionsObj;
139 }
140 
141 // Statistical methods -----------------------------
142 template <class P_V,class P_M, class Q_V, class Q_M>
143 bool
145 {
146  return m_optionsObj->m_computeSolution;
147 }
148 //--------------------------------------------------
149 template <class P_V,class P_M,class Q_V,class Q_M>
150 void
152  const McOptionsValues* alternativeOptionsValues)
153 {
154  m_env.fullComm().Barrier();
155  m_env.fullComm().syncPrintDebugMsg("Entering StatisticalForwardProblem<P_V,P_M>::solveWithMonteCarlo()",1,3000000);
156 
157  if (m_optionsObj->m_computeSolution == false) {
158  if ((m_env.subDisplayFile())) {
159  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
160  << ": avoiding solution, as requested by user"
161  << std::endl;
162  }
163  return;
164  }
165  if ((m_env.subDisplayFile())) {
166  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
167  << ": computing solution, as requested by user"
168  << std::endl;
169  }
170 
171  if (m_solutionPdf ) delete m_solutionPdf;
172 
173 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
174  if (m_unifiedSolutionCdf) delete m_unifiedSolutionCdf;
175  if (m_unifiedCdfValues ) delete m_unifiedCdfValues;
176  if (m_unifiedCdfGrids ) delete m_unifiedCdfGrids;
177 
178  if (m_subSolutionCdf ) delete m_subSolutionCdf;
179  if (m_subCdfValues ) delete m_subCdfValues;
180  if (m_subCdfGrids ) delete m_subCdfGrids;
181 
182  if (m_subSolutionMdf ) delete m_subSolutionMdf;
183 #endif
184 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
185  if (m_subMdfValues ) delete m_subMdfValues;
186  if (m_subMdfGrids ) delete m_subMdfGrids;
187 #endif
188  if (m_solutionRealizer ) delete m_solutionRealizer;
189 
190  if (m_mcSeqGenerator ) delete m_mcSeqGenerator;
191 
192  if (m_qoiChain) {
193  m_qoiChain->clear();
194  delete m_qoiChain;
195  }
196 
197  if (m_paramChain) {
198  m_paramChain->clear();
199  delete m_paramChain;
200  }
201 
202  Q_V numEvaluationPointsVec(m_qoiRv.imageSet().vectorSpace().zeroVector());
203  numEvaluationPointsVec.cwSet(250.);
204 
205  // Compute output realizer: Monte Carlo approach
206  m_paramChain = new SequenceOfVectors<P_V,P_M>(m_paramRv.imageSet().vectorSpace(),0,m_optionsObj->m_prefix+"paramChain");
207  m_qoiChain = new SequenceOfVectors<Q_V,Q_M>(m_qoiRv.imageSet().vectorSpace(), 0,m_optionsObj->m_prefix+"qoiChain" );
208  m_mcSeqGenerator = new MonteCarloSG<P_V,P_M,Q_V,Q_M>(m_optionsObj->m_prefix.c_str(),
209  alternativeOptionsValues,
210  m_paramRv,
211  m_qoiFunction);
212  //m_qoiRv);
213  m_mcSeqGenerator->generateSequence(*m_paramChain,
214  *m_qoiChain);
215  m_solutionRealizer = new SequentialVectorRealizer<Q_V,Q_M>((m_optionsObj->m_prefix+"Qoi").c_str(),
216  *m_qoiChain);
217  m_qoiRv.setRealizer(*m_solutionRealizer);
218 
219  // Compute output mdf: uniform sampling approach
220 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
221  m_subMdfGrids = new ArrayOfOneDGrids <Q_V,Q_M>((m_optionsObj->m_prefix+"QoiMdf_").c_str(),m_qoiRv.imageSet().vectorSpace());
222  m_subMdfValues = new ArrayOfOneDTables<Q_V,Q_M>((m_optionsObj->m_prefix+"QoiMdf_").c_str(),m_qoiRv.imageSet().vectorSpace());
223  m_qoiChain->subUniformlySampledMdf(numEvaluationPointsVec, // input
224  *m_subMdfGrids, // output
225  *m_subMdfValues); // output
226 
227  m_subSolutionMdf = new SampledVectorMdf<Q_V,Q_M>((m_optionsObj->m_prefix+"Qoi").c_str(),
228  *m_subMdfGrids,
229  *m_subMdfValues);
230  m_qoiRv.setMdf(*m_subSolutionMdf);
231 #endif
232 
233  // Compute output cdf: uniform sampling approach
234 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
235  std::string subCoreName_qoiCdf(m_optionsObj->m_prefix+ "QoiCdf_");
236  std::string uniCoreName_qoiCdf(m_optionsObj->m_prefix+"unifQoiCdf_");
237  if (m_env.numSubEnvironments() == 1) subCoreName_qoiCdf = uniCoreName_qoiCdf;
238 
239  std::string subCoreName_solutionCdf(m_optionsObj->m_prefix+ "Qoi");
240  std::string uniCoreName_solutionCdf(m_optionsObj->m_prefix+"unifQoi");
241  if (m_env.numSubEnvironments() == 1) subCoreName_solutionCdf = uniCoreName_solutionCdf;
242 
243  m_subCdfGrids = new ArrayOfOneDGrids <Q_V,Q_M>(subCoreName_qoiCdf.c_str(),m_qoiRv.imageSet().vectorSpace());
244  m_subCdfValues = new ArrayOfOneDTables<Q_V,Q_M>(subCoreName_qoiCdf.c_str(),m_qoiRv.imageSet().vectorSpace());
245  m_qoiChain->subUniformlySampledCdf(numEvaluationPointsVec, // input
246  *m_subCdfGrids, // output
247  *m_subCdfValues); // output
248 
249  m_subSolutionCdf = new SampledVectorCdf<Q_V,Q_M>(subCoreName_solutionCdf.c_str(),
250  *m_subCdfGrids,
251  *m_subCdfValues);
252  m_qoiRv.setSubCdf(*m_subSolutionCdf);
253 
254  // Compute unified cdf if necessary
255  if (m_env.numSubEnvironments() == 1) {
256  m_qoiRv.setUnifiedCdf(*m_subSolutionCdf);
257  }
258  else {
259  m_unifiedCdfGrids = new ArrayOfOneDGrids <Q_V,Q_M>(uniCoreName_qoiCdf.c_str(),m_qoiRv.imageSet().vectorSpace());
260  m_unifiedCdfValues = new ArrayOfOneDTables<Q_V,Q_M>(uniCoreName_qoiCdf.c_str(),m_qoiRv.imageSet().vectorSpace());
261  m_qoiChain->unifiedUniformlySampledCdf(numEvaluationPointsVec, // input
262  *m_unifiedCdfGrids, // output
263  *m_unifiedCdfValues); // output
264 
265  m_unifiedSolutionCdf = new SampledVectorCdf<Q_V,Q_M>(uniCoreName_solutionCdf.c_str(),
266  *m_unifiedCdfGrids,
267  *m_unifiedCdfValues);
268  m_qoiRv.setUnifiedCdf(*m_unifiedSolutionCdf);
269  }
270 #endif
271  // Compute (just unified one) covariance matrix, if requested
272  // Compute (just unified one) correlation matrix, if requested
273  P_M* pqCovarianceMatrix = NULL;
274  P_M* pqCorrelationMatrix = NULL;
275  if (m_optionsObj->m_computeCovariances || m_optionsObj->m_computeCorrelations) {
276  if (m_env.subDisplayFile()) {
277  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
278  << ", prefix = " << m_optionsObj->m_prefix
279  << ": instantiating cov and corr matrices"
280  << std::endl;
281  }
282 
283  // Only compute correlations on the inter0Comm communicator
284  if (m_env.subRank() == 0) {
285  pqCovarianceMatrix = new P_M(m_env,
286  m_paramRv.imageSet().vectorSpace().map(), // number of rows
287  m_qoiRv.imageSet().vectorSpace().dimGlobal()); // number of cols
288  pqCorrelationMatrix = new P_M(m_env,
289  m_paramRv.imageSet().vectorSpace().map(), // number of rows
290  m_qoiRv.imageSet().vectorSpace().dimGlobal()); // number of cols
292  *m_qoiChain,
293  std::min(m_paramRv.realizer().subPeriod(),m_qoiRv.realizer().subPeriod()), // FIX ME: might be INFINITY
294  *pqCovarianceMatrix,
295  *pqCorrelationMatrix);
296  }
297  }
298 
299  // Write data out
300  if (m_env.subDisplayFile()) {
301  if (pqCovarianceMatrix ) {
302  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
303  << ", prefix = " << m_optionsObj->m_prefix
304  << ": contents of covariance matrix are\n" << *pqCovarianceMatrix
305  << std::endl;
306  }
307  if (pqCorrelationMatrix) {
308  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
309  << ", prefix = " << m_optionsObj->m_prefix
310  << ": contents of correlation matrix are\n" << *pqCorrelationMatrix
311  << std::endl;
312  }
313  }
314 
315  // Open data output file
316  if (m_env.subDisplayFile()) {
317  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
318  << ", prefix = " << m_optionsObj->m_prefix
319  << ": checking necessity of opening data output file '" << m_optionsObj->m_dataOutputFileName
320  << "'"
321  << std::endl;
322  }
323  FilePtrSetStruct filePtrSet;
324  if (m_env.openOutputFile(m_optionsObj->m_dataOutputFileName,
325  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
326  m_optionsObj->m_dataOutputAllowedSet,
327  false,
328  filePtrSet)) {
329 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
330  m_qoiRv.mdf().print(*filePtrSet.ofsVar);
331 #endif
332 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
333  *filePtrSet.ofsVar << m_qoiRv.subCdf();
334 #endif
335 
336  //if (pqCovarianceMatrix ) *filePtrSet.ofsVar << *pqCovarianceMatrix; // FIX ME: output matrix in matlab format
337  //if (pqCorrelationMatrix) *filePtrSet.ofsVar << *pqCorrelationMatrix; // FIX ME: output matrix in matlab format
338 
339  // Write unified cdf if necessary
340 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
341  if (m_env.numSubEnvironments() > 1) {
342  if (m_qoiRv.imageSet().vectorSpace().numOfProcsForStorage() == 1) {
343  if (m_env.inter0Rank() == 0) {
344  *filePtrSet.ofsVar << m_qoiRv.unifiedCdf(); //*m_unifiedSolutionCdf;
345  }
346  }
347  else {
348  queso_error_msg("unified cdf writing, parallel vectors not supported yet");
349  }
350  }
351 #endif
352  // Close data output file
353  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
354  if (m_env.subDisplayFile()) {
355  *m_env.subDisplayFile() << "In StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::solveWithMonteCarlo()"
356  << ", prefix = " << m_optionsObj->m_prefix
357  << ": closed data output file '" << m_optionsObj->m_dataOutputFileName
358  << "'"
359  << std::endl;
360  }
361  }
362  if (m_env.subDisplayFile()) {
363  *m_env.subDisplayFile() << std::endl;
364  }
365 
366  if (pqCovarianceMatrix ) delete pqCovarianceMatrix;
367  if (pqCorrelationMatrix) delete pqCorrelationMatrix;
368 
369  m_env.fullComm().syncPrintDebugMsg("Leaving StatisticalForwardProblem<P_V,P_M>::solveWithMonteCarlo()",1,3000000);
370  m_env.fullComm().Barrier();
371 
372  return;
373 }
374 //--------------------------------------------------
375 template <class P_V,class P_M,class Q_V,class Q_M>
378 {
379  return m_qoiRv;
380 }
381 //--------------------------------------------------
382 
383 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
384 template <class P_V,class P_M,class Q_V,class Q_M>
387 {
388  if (m_env.numSubEnvironments() == 1) {
389  return m_qoiRv.subCdf();
390  }
391 
392  if (m_env.inter0Rank() < 0) {
393  return m_qoiRv.subCdf();
394  }
395 
396 
397  // m_env.worldRank(),
398  // "StatisticalForwardProblem<P_V,P_M,Q_V,Q_M>::qoiRv_unifiedCdf()",
399  // "variable is NULL");
400  return m_qoiRv.unifiedCdf(); //*m_unifiedSolutionCdf;
401 }
402 #endif
403 //--------------------------------------------------
404 template <class P_V,class P_M,class Q_V,class Q_M>
407 {
408 
409  // Make sure this runs after the forward propagation
410  // only then we obtain the actual realizations of the parameters
411  queso_require_msg(m_paramChain, "m_paramChain is NULL");
412 
413  return *m_paramChain;
414 
415 }
416 // I/O methods--------------------------------------
417 template <class P_V,class P_M,class Q_V,class Q_M>
418 void
420 {
421  return;
422 }
423 
424 } // End namespace QUESO
425 
const GenericVectorRV< Q_V, Q_M > & qoiRv() const
Returns the QoI RV; access to private attribute m_qoiRv.
A class for handling sampled vector CDFs.
This class provides options for a Statistical Forward Problem if no input file is available...
std::string m_prefix
Definition: VectorCdf.h:99
void unifiedUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &unifiedCdfGrids, ArrayOfOneDTables< V, M > &unifiedCdfValues) const
Uniformly samples from the CDF from the sub-sequence.
This class provides options for the Monte Carlo sequence generator if no input file is available...
void subUniformlySampledMdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &mdfGrids, ArrayOfOneDTables< V, M > &mdfValues) const
A templated class that implements a Monte Carlo generator of samples.
Definition: MonteCarloSG.h:52
This templated class represents a Statistical Forward Problem.
Class for handling vector samples (sequence of vectors).
bool computeSolutionFlag() const
Whether or not compute the solution.
void solveWithMonteCarlo(const McOptionsValues *alternativeOptionsValues)
Solves the problem through Monte Carlo algorithm.
const BaseVectorCdf< Q_V, Q_M > & qoiRv_unifiedCdf() const
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
void subUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &cdfGrids, ArrayOfOneDTables< V, M > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:86
const BaseVectorSequence< Q_V, Q_M > & getParamChain() const
Returns the parameter chain; access to private attribute m_paramChain.
A class for handling sampled vector MDFs.
Struct for handling data input and output from files.
Definition: Environment.h:77
StatisticalForwardProblem(const char *prefix, const SfpOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > &paramRv, const BaseVectorFunction< P_V, P_M, Q_V, Q_M > &qoiFunction, GenericVectorRV< Q_V, Q_M > &qoiRv)
Constructor.
void ComputeCovCorrMatricesBetweenVectorSequences(const BaseVectorSequence< P_V, P_M > &subPSeq, const BaseVectorSequence< Q_V, Q_M > &subQSeq, unsigned int subNumSamples, P_M &pqCovMatrix, P_M &pqCorrMatrix)
A class for handling sequential draws (sampling) from probability density distributions.
void print(std::ostream &os) const
TODO: Prints the sequence.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

Generated on Tue Jun 5 2018 19:48:56 for queso-0.57.1 by  doxygen 1.8.5