queso-0.52.0
MetropolisHastingsSG.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-2015 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/asserts.h>
26 #include <queso/MetropolisHastingsSG.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 
30 #include <queso/HessianCovMatricesTKGroup.h>
31 #include <queso/ScaledCovMatrixTKGroup.h>
32 #include <queso/TransformedScaledCovMatrixTKGroup.h>
33 
34 #include <queso/InvLogitGaussianJointPdf.h>
35 
36 #include <queso/GaussianJointPdf.h>
37 
38 namespace QUESO {
39 
40 // Default constructor -----------------------------
41 MHRawChainInfoStruct::MHRawChainInfoStruct()
42 {
43  reset();
44 }
45 // Copy constructor----------------------------------
46 MHRawChainInfoStruct::MHRawChainInfoStruct(const MHRawChainInfoStruct& rhs)
47 {
48  this->copy(rhs);
49 }
50 // Destructor ---------------------------------------
51 MHRawChainInfoStruct::~MHRawChainInfoStruct()
52 {
53 }
54 // Set methods---------------------------------------
55 MHRawChainInfoStruct&
56 MHRawChainInfoStruct::operator=(const MHRawChainInfoStruct& rhs)
57 {
58  this->copy(rhs);
59  return *this;
60 }
61 //---------------------------------------------------
62 MHRawChainInfoStruct&
63 MHRawChainInfoStruct::operator+=(const MHRawChainInfoStruct& rhs)
64 {
65  runTime += rhs.runTime;
66  candidateRunTime += rhs.candidateRunTime;
67  targetRunTime += rhs.targetRunTime;
68  mhAlphaRunTime += rhs.mhAlphaRunTime;
69  drAlphaRunTime += rhs.drAlphaRunTime;
70  drRunTime += rhs.drRunTime;
71  amRunTime += rhs.amRunTime;
72 
73  numTargetCalls += rhs.numTargetCalls;
74  numDRs += rhs.numDRs;
75  numOutOfTargetSupport += rhs.numOutOfTargetSupport;
76  numOutOfTargetSupportInDR += rhs.numOutOfTargetSupportInDR;
77  numRejections += rhs.numRejections;
78 
79  return *this;
80 }
81 // Misc methods--------------------------------------------------
82 void
83 MHRawChainInfoStruct::reset()
84 {
85  runTime = 0.;
86  candidateRunTime = 0.;
87  targetRunTime = 0.;
88  mhAlphaRunTime = 0.;
89  drAlphaRunTime = 0.;
90  drRunTime = 0.;
91  amRunTime = 0.;
92 
93  numTargetCalls = 0;
94  numDRs = 0;
95  numOutOfTargetSupport = 0;
96  numOutOfTargetSupportInDR = 0;
97  numRejections = 0;
98 }
99 //---------------------------------------------------
100 void
101 MHRawChainInfoStruct::copy(const MHRawChainInfoStruct& rhs)
102 {
103  runTime = rhs.runTime;
104  candidateRunTime = rhs.candidateRunTime;
105  targetRunTime = rhs.targetRunTime;
106  mhAlphaRunTime = rhs.mhAlphaRunTime;
107  drAlphaRunTime = rhs.drAlphaRunTime;
108  drRunTime = rhs.drRunTime;
109  amRunTime = rhs.amRunTime;
110 
111  numTargetCalls = rhs.numTargetCalls;
112  numDRs = rhs.numDRs;
113  numOutOfTargetSupport = rhs.numOutOfTargetSupport;
114  numOutOfTargetSupportInDR = rhs.numOutOfTargetSupportInDR;
115  numRejections = rhs.numRejections;
116 
117  return;
118 }
119 //---------------------------------------------------
120 void
121 MHRawChainInfoStruct::mpiSum(const MpiComm& comm, MHRawChainInfoStruct& sumInfo) const
122 {
123  comm.Allreduce((void *) &runTime, (void *) &sumInfo.runTime, (int) 7, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
124  "MHRawChainInfoStruct::mpiSum()",
125  "failed MPI.Allreduce() for sum of doubles");
126 
127  comm.Allreduce((void *) &numTargetCalls, (void *) &sumInfo.numTargetCalls, (int) 5, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
128  "MHRawChainInfoStruct::mpiSum()",
129  "failed MPI.Allreduce() for sum of unsigned ints");
130 
131  return;
132 }
133 
134 // Default constructor -----------------------------
135 template<class P_V,class P_M>
136 MetropolisHastingsSG<P_V,P_M>::MetropolisHastingsSG( const char* prefix, const MhOptionsValues* alternativeOptionsValues, // dakota const BaseVectorRV<P_V,P_M>& sourceRv, const P_V& initialPosition, const P_M* inputProposalCovMatrix)
142  :
143  m_env (sourceRv.env()),
144  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
145  m_targetPdf (sourceRv.pdf()),
146  m_initialPosition (initialPosition),
147  m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
148  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
149  m_numDisabledParameters (0), // gpmsa2
150  m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true), // gpmsa2
151  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
152  m_tk (NULL),
153  m_positionIdForDebugging (0),
154  m_stageIdForDebugging (0),
155  m_idsOfUniquePositions (0),//0.),
156  m_logTargets (0),//0.),
157  m_alphaQuotients (0),//0.),
158  m_lastChainSize (0),
159  m_lastMean (NULL),
160  m_lastAdaptedCovMatrix (NULL),
161  m_numPositionsNotSubWritten (0),
162 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
163  m_alternativeOptionsValues (NULL,NULL),
164 #else
165  m_alternativeOptionsValues (),
166 #endif
167  m_optionsObj (NULL),
168  m_computeInitialPriorAndLikelihoodValues(true),
169  m_initialLogPriorValue (0.),
170  m_initialLogLikelihoodValue (0.)
171 {
172  if (inputProposalCovMatrix != NULL) {
173  m_initialProposalCovMatrix = *inputProposalCovMatrix;
174  }
175  if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
176  if (m_env.optionsInputFileName() == "") {
177  m_optionsObj = new MetropolisHastingsSGOptions(m_env,prefix,m_alternativeOptionsValues);
178  }
179  else {
180  m_optionsObj = new MetropolisHastingsSGOptions(m_env,prefix);
181  m_optionsObj->scanOptionsValues();
182  }
183 
184  if ((m_env.subDisplayFile() ) &&
185  (m_optionsObj->m_ov.m_totallyMute == false)) {
186  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::constructor(1)"
187  << ": prefix = " << prefix
188  << ", alternativeOptionsValues = " << alternativeOptionsValues
189  << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
190  << ", m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
191  << std::endl;
192  }
193 
194  UQ_FATAL_TEST_MACRO(sourceRv.imageSet().vectorSpace().dimLocal() != initialPosition.sizeLocal(),
195  m_env.worldRank(),
196  "MetropolisHastingsSG<P_V,P_M>::constructor(1)",
197  "'sourceRv' and 'initialPosition' should have equal dimensions");
198 
199  if (inputProposalCovMatrix) {
200  UQ_FATAL_TEST_MACRO(sourceRv.imageSet().vectorSpace().dimLocal() != inputProposalCovMatrix->numRowsLocal(),
201  m_env.worldRank(),
202  "MetropolisHastingsSG<P_V,P_M>::constructor(1)",
203  "'sourceRv' and 'inputProposalCovMatrix' should have equal dimensions");
204  UQ_FATAL_TEST_MACRO(inputProposalCovMatrix->numCols() != inputProposalCovMatrix->numRowsGlobal(),
205  m_env.worldRank(),
206  "MetropolisHastingsSG<P_V,P_M>::constructor(1)",
207  "'inputProposalCovMatrix' should be a square matrix");
208  }
209 
210  commonConstructor();
211 
212  if ((m_env.subDisplayFile() ) &&
213  (m_optionsObj->m_ov.m_totallyMute == false)) {
214  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::constructor(1)"
215  << std::endl;
216  }
217 }
218 template<class P_V,class P_M>
219 MetropolisHastingsSG<P_V,P_M>::MetropolisHastingsSG( const char* prefix, const MhOptionsValues* alternativeOptionsValues, // dakota const BaseVectorRV<P_V,P_M>& sourceRv, const P_V& initialPosition, double initialLogPrior, double initialLogLikelihood, const P_M* inputProposalCovMatrix)
227  :
228  m_env (sourceRv.env()),
229  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
230  m_targetPdf (sourceRv.pdf()),
231  m_initialPosition (initialPosition),
232  m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
233  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
234  m_numDisabledParameters (0), // gpmsa2
235  m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true), // gpmsa2
236  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
237  m_tk (NULL),
238  m_positionIdForDebugging (0),
239  m_stageIdForDebugging (0),
240  m_idsOfUniquePositions (0),//0.),
241  m_logTargets (0),//0.),
242  m_alphaQuotients (0),//0.),
243  m_lastChainSize (0),
244  m_lastMean (NULL),
245  m_lastAdaptedCovMatrix (NULL),
246  m_numPositionsNotSubWritten (0),
247 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
248  m_alternativeOptionsValues (NULL,NULL),
249 #else
250  m_alternativeOptionsValues (),
251 #endif
252  m_optionsObj (NULL),
253  m_computeInitialPriorAndLikelihoodValues(false),
254  m_initialLogPriorValue (initialLogPrior),
255  m_initialLogLikelihoodValue (initialLogLikelihood)
256 {
257  if (inputProposalCovMatrix != NULL) {
258  m_initialProposalCovMatrix = *inputProposalCovMatrix;
259  }
260  if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
261  if (m_env.optionsInputFileName() == "") {
262  m_optionsObj = new MetropolisHastingsSGOptions(m_env,prefix,m_alternativeOptionsValues);
263  }
264  else {
265  m_optionsObj = new MetropolisHastingsSGOptions(m_env,prefix);
266  m_optionsObj->scanOptionsValues();
267  }
268 
269  if ((m_env.subDisplayFile() ) &&
270  (m_optionsObj->m_ov.m_totallyMute == false)) {
271  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::constructor(2)"
272  << ": prefix = " << prefix
273  << ", alternativeOptionsValues = " << alternativeOptionsValues
274  << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
275  << ", m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
276  << std::endl;
277  }
278 
279  UQ_FATAL_TEST_MACRO(sourceRv.imageSet().vectorSpace().dimLocal() != initialPosition.sizeLocal(),
280  m_env.worldRank(),
281  "MetropolisHastingsSG<P_V,P_M>::constructor(2)",
282  "'sourceRv' and 'initialPosition' should have equal dimensions");
283 
284  if (inputProposalCovMatrix) {
285  UQ_FATAL_TEST_MACRO(sourceRv.imageSet().vectorSpace().dimLocal() != inputProposalCovMatrix->numRowsLocal(),
286  m_env.worldRank(),
287  "MetropolisHastingsSG<P_V,P_M>::constructor(2)",
288  "'sourceRv' and 'inputProposalCovMatrix' should have equal dimensions");
289  UQ_FATAL_TEST_MACRO(inputProposalCovMatrix->numCols() != inputProposalCovMatrix->numRowsGlobal(),
290  m_env.worldRank(),
291  "MetropolisHastingsSG<P_V,P_M>::constructor(2)",
292  "'inputProposalCovMatrix' should be a square matrix");
293  }
294 
295  commonConstructor();
296 
297  if ((m_env.subDisplayFile() ) &&
298  (m_optionsObj->m_ov.m_totallyMute == false)) {
299  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::constructor(2)"
300  << std::endl;
301  }
302 }
303 // Constructor -------------------------------------
304 template<class P_V,class P_M>
305 MetropolisHastingsSG<P_V,P_M>::MetropolisHastingsSG(
306  const MLSamplingLevelOptions& mlOptions,
307  const BaseVectorRV<P_V,P_M>& sourceRv,
308  const P_V& initialPosition, // KEY
309  const P_M* inputProposalCovMatrix)
310  :
311  m_env (sourceRv.env()),
312  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
313  m_targetPdf (sourceRv.pdf()),
314  m_initialPosition (initialPosition),
315  m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
316  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
317  m_numDisabledParameters (0), // gpmsa2
318  m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true), // gpmsa2
319  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
320  m_tk (NULL),
321  m_positionIdForDebugging (0),
322  m_stageIdForDebugging (0),
323  m_idsOfUniquePositions (0),//0.),
324  m_logTargets (0),//0.),
325  m_alphaQuotients (0),//0.),
326  m_lastChainSize (0),
327  m_lastMean (NULL),
328  m_lastAdaptedCovMatrix (NULL),
329 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
330  m_alternativeOptionsValues (NULL,NULL),
331 #else
332  m_alternativeOptionsValues (),
333 #endif
334  m_optionsObj (new MetropolisHastingsSGOptions(mlOptions)),
335  m_computeInitialPriorAndLikelihoodValues(true),
336  m_initialLogPriorValue (0.),
337  m_initialLogLikelihoodValue (0.)
338 {
339  if (inputProposalCovMatrix != NULL) {
340  m_initialProposalCovMatrix = *inputProposalCovMatrix;
341  if ((m_env.subDisplayFile() ) &&
342  (m_optionsObj->m_ov.m_totallyMute == false)) {
343  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::constructor(3)"
344  << ": just set m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
345  << std::endl;
346  }
347  }
348  if ((m_env.subDisplayFile() ) &&
349  (m_optionsObj->m_ov.m_totallyMute == false)) {
350  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::constructor(3)"
351  << std::endl;
352  }
353 
354  commonConstructor();
355 
356  if ((m_env.subDisplayFile() ) &&
357  (m_optionsObj->m_ov.m_totallyMute == false)) {
358  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::constructor(3)"
359  << std::endl;
360  }
361 }
362 // Constructor -------------------------------------
363 template<class P_V,class P_M>
364 MetropolisHastingsSG<P_V,P_M>::MetropolisHastingsSG(
365  const MLSamplingLevelOptions& mlOptions,
366  const BaseVectorRV<P_V,P_M>& sourceRv,
367  const P_V& initialPosition, // KEY
368  double initialLogPrior,
369  double initialLogLikelihood,
370  const P_M* inputProposalCovMatrix)
371  :
372  m_env (sourceRv.env()),
373  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
374  m_targetPdf (sourceRv.pdf()),
375  m_initialPosition (initialPosition),
376  m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
377  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
378  m_numDisabledParameters (0), // gpmsa2
379  m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true), // gpmsa2
380  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
381  m_tk (NULL),
382  m_positionIdForDebugging (0),
383  m_stageIdForDebugging (0),
384  m_idsOfUniquePositions (0),//0.),
385  m_logTargets (0),//0.),
386  m_alphaQuotients (0),//0.),
387  m_lastChainSize (0),
388  m_lastMean (NULL),
389  m_lastAdaptedCovMatrix (NULL),
390 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
391  m_alternativeOptionsValues (NULL,NULL),
392 #else
393  m_alternativeOptionsValues (),
394 #endif
395  m_optionsObj (new MetropolisHastingsSGOptions(mlOptions)),
396  m_computeInitialPriorAndLikelihoodValues(false),
397  m_initialLogPriorValue (initialLogPrior),
398  m_initialLogLikelihoodValue (initialLogLikelihood)
399 {
400  if (inputProposalCovMatrix != NULL) {
401  m_initialProposalCovMatrix = *inputProposalCovMatrix;
402  if ((m_env.subDisplayFile() ) &&
403  (m_optionsObj->m_ov.m_totallyMute == false)) {
404  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::constructor(4)"
405  << ": just set m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
406  << std::endl;
407  }
408  }
409  if ((m_env.subDisplayFile() ) &&
410  (m_optionsObj->m_ov.m_totallyMute == false)) {
411  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::constructor(4)"
412  << std::endl;
413  }
414 
415  commonConstructor();
416 
417  if ((m_env.subDisplayFile() ) &&
418  (m_optionsObj->m_ov.m_totallyMute == false)) {
419  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::constructor(4)"
420  << std::endl;
421  }
422 }
423 // Destructor ---------------------------------------
424 template<class P_V,class P_M>
425 MetropolisHastingsSG<P_V,P_M>::~MetropolisHastingsSG()
426 {
427  //if (m_env.subDisplayFile()) {
428  // *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::destructor()"
429  // << std::endl;
430  //}
431 
432  if (m_lastAdaptedCovMatrix) delete m_lastAdaptedCovMatrix;
433  if (m_lastMean) delete m_lastMean;
434  m_lastChainSize = 0;
435  m_rawChainInfo.reset();
436  m_alphaQuotients.clear();
437  m_logTargets.clear();
438  m_numDisabledParameters = 0; // gpmsa2
439  m_parameterEnabledStatus.clear(); // gpmsa2
440  m_positionIdForDebugging = 0;
441  m_stageIdForDebugging = 0;
442  m_idsOfUniquePositions.clear();
443 
444  if (m_tk ) delete m_tk;
445  if (m_targetPdfSynchronizer) delete m_targetPdfSynchronizer;
446  if (m_optionsObj ) delete m_optionsObj;
447 
448  //if (m_env.subDisplayFile()) {
449  // *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::destructor()"
450  // << std::endl;
451  //}
452 }
453 
454 // Private methods----------------------------------
455 template<class P_V,class P_M>
456 void
457 MetropolisHastingsSG<P_V,P_M>::commonConstructor()
458 {
460  // Instantiate the appropriate TK (transition kernel)
462  if ((m_env.subDisplayFile() ) &&
463  (m_optionsObj->m_ov.m_totallyMute == false)) {
464  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
465  << std::endl;
466  }
467 
468  if (m_optionsObj->m_ov.m_initialPositionDataInputFileName != ".") { // palms
469  std::set<unsigned int> tmpSet;
470  tmpSet.insert(m_env.subId());
471  m_initialPosition.subReadContents((m_optionsObj->m_ov.m_initialPositionDataInputFileName+"_sub"+m_env.subIdString()),
472  m_optionsObj->m_ov.m_initialPositionDataInputFileType,
473  tmpSet);
474  if ((m_env.subDisplayFile() ) &&
475  (m_optionsObj->m_ov.m_totallyMute == false)) {
476  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
477  << ": just read initial position contents = " << m_initialPosition
478  << std::endl;
479  }
480  }
481 
482  if (m_optionsObj->m_ov.m_parameterDisabledSet.size() > 0) { // gpmsa2
483  for (std::set<unsigned int>::iterator setIt = m_optionsObj->m_ov.m_parameterDisabledSet.begin(); setIt != m_optionsObj->m_ov.m_parameterDisabledSet.end(); ++setIt) {
484  unsigned int paramId = *setIt;
485  if (paramId < m_vectorSpace.dimLocal()) {
486  m_numDisabledParameters++;
487  m_parameterEnabledStatus[paramId] = false;
488  }
489  }
490  }
491 
492  std::vector<double> drScalesAll(m_optionsObj->m_ov.m_drScalesForExtraStages.size()+1,1.);
493  for (unsigned int i = 1; i < (m_optionsObj->m_ov.m_drScalesForExtraStages.size()+1); ++i) {
494  drScalesAll[i] = m_optionsObj->m_ov.m_drScalesForExtraStages[i-1];
495  }
496  if (m_optionsObj->m_ov.m_tkUseLocalHessian) { // sep2011
497  m_tk = new HessianCovMatricesTKGroup<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
498  m_vectorSpace,
499  drScalesAll,
500  *m_targetPdfSynchronizer);
501  if ((m_env.subDisplayFile() ) &&
502  (m_optionsObj->m_ov.m_totallyMute == false)) {
503  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
504  << ": just instantiated a 'HessianCovMatrices' TK class"
505  << std::endl;
506  }
507  }
508  else {
509  if (m_optionsObj->m_ov.m_initialProposalCovMatrixDataInputFileName != ".") { // palms
510  std::set<unsigned int> tmpSet;
511  tmpSet.insert(m_env.subId());
512  m_initialProposalCovMatrix.subReadContents((m_optionsObj->m_ov.m_initialProposalCovMatrixDataInputFileName+"_sub"+m_env.subIdString()),
513  m_optionsObj->m_ov.m_initialProposalCovMatrixDataInputFileType,
514  tmpSet);
515  if ((m_env.subDisplayFile() ) &&
516  (m_optionsObj->m_ov.m_totallyMute == false)) {
517  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
518  << ": just read initial proposal cov matrix contents = " << m_initialProposalCovMatrix
519  << std::endl;
520  }
521  }
522  else {
523  UQ_FATAL_TEST_MACRO(m_nullInputProposalCovMatrix,
524  m_env.worldRank(),
525  "MetropolisHastingsSG<P_V,P_M>::commonConstructor()",
526  "proposal cov matrix should have been passed by user, since, according to the input algorithm options, local Hessians will not be used in the proposal");
527  }
528 
529  // Decide whether or not to do logit transform
530  if (m_optionsObj->m_ov.m_doLogitTransform) {
531  // Variable transform initial proposal cov matrix
532  transformInitialCovMatrixToGaussianSpace(
533  dynamic_cast<const BoxSubset<P_V, P_M> & >(m_targetPdf.domainSet()));
534 
535  // We need this dynamic_cast to BoxSubset so that m_tk can inspect the
536  // domain bounds and do the necessary transform
537  m_tk = new TransformedScaledCovMatrixTKGroup<P_V, P_M>(
538  m_optionsObj->m_prefix.c_str(),
539  dynamic_cast<const BoxSubset<P_V, P_M> & >(m_targetPdf.domainSet()),
540  drScalesAll, m_initialProposalCovMatrix);
541  }
542  else {
543  m_tk = new ScaledCovMatrixTKGroup<P_V, P_M>(
544  m_optionsObj->m_prefix.c_str(), m_vectorSpace, drScalesAll,
545  m_initialProposalCovMatrix);
546  }
547 
548  if ((m_env.subDisplayFile() ) &&
549  (m_optionsObj->m_ov.m_totallyMute == false)) {
550  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
551  << ": just instantiated a 'ScaledCovMatrix' TK class"
552  << std::endl;
553  }
554  }
555 
556  if ((m_env.subDisplayFile() ) &&
557  (m_optionsObj->m_ov.m_totallyMute == false)) {
558  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
559  << std::endl;
560  }
561  return;
562 }
563 //--------------------------------------------------
564 template<class P_V,class P_M>
565 double
566 MetropolisHastingsSG<P_V,P_M>::alpha(
567  const MarkovChainPositionData<P_V>& x,
568  const MarkovChainPositionData<P_V>& y,
569  unsigned int xStageId,
570  unsigned int yStageId,
571  double* alphaQuotientPtr)
572 {
573  double alphaQuotient = 0.;
574  if ((x.outOfTargetSupport() == false) &&
575  (y.outOfTargetSupport() == false)) {
576  if ((x.logTarget() == -INFINITY) ||
577  (x.logTarget() == INFINITY) ||
578  ( (boost::math::isnan)(x.logTarget()) )) {
579  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
580  << ", worldRank " << m_env.worldRank()
581  << ", fullRank " << m_env.fullRank()
582  << ", subEnvironment " << m_env.subId()
583  << ", subRank " << m_env.subRank()
584  << ", inter0Rank " << m_env.inter0Rank()
585  << ", positionId = " << m_positionIdForDebugging
586  << ", stageId = " << m_stageIdForDebugging
587  << ": x.logTarget() = " << x.logTarget()
588  << ", x.values() = " << x.vecValues()
589  << ", y.values() = " << y.vecValues()
590  << std::endl;
591  }
592  else if ((y.logTarget() == -INFINITY ) ||
593  (y.logTarget() == INFINITY ) ||
594  ( (boost::math::isnan)(y.logTarget()) )) {
595  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
596  << ", worldRank " << m_env.worldRank()
597  << ", fullRank " << m_env.fullRank()
598  << ", subEnvironment " << m_env.subId()
599  << ", subRank " << m_env.subRank()
600  << ", inter0Rank " << m_env.inter0Rank()
601  << ", positionId = " << m_positionIdForDebugging
602  << ", stageId = " << m_stageIdForDebugging
603  << ": y.logTarget() = " << y.logTarget()
604  << ", x.values() = " << x.vecValues()
605  << ", y.values() = " << y.vecValues()
606  << std::endl;
607  }
608  else {
609  double yLogTargetToUse = y.logTarget();
610 
611  if (m_tk->symmetric()) {
612  alphaQuotient = std::exp(yLogTargetToUse - x.logTarget());
613 
614  if ((m_env.subDisplayFile() ) &&
615  (m_env.displayVerbosity() >= 3 ) &&
616  (m_optionsObj->m_ov.m_totallyMute == false)) {
617  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
618  << ": symmetric proposal case"
619  << ", x = " << x.vecValues()
620  << ", y = " << y.vecValues()
621  << ", yLogTargetToUse = " << yLogTargetToUse
622  << ", x.logTarget() = " << x.logTarget()
623  << ", alpha = " << alphaQuotient
624  << std::endl;
625  }
626  }
627  else {
628  double qyx = m_tk->rv(yStageId).pdf().lnValue(x.vecValues(),NULL,NULL,NULL,NULL);
629  if ((m_env.subDisplayFile() ) &&
630  (m_env.displayVerbosity() >= 10 ) &&
631  (m_optionsObj->m_ov.m_totallyMute == false)) {
632  const InvLogitGaussianJointPdf<P_V,P_M>* pdfYX = dynamic_cast< const InvLogitGaussianJointPdf<P_V,P_M>* >(&(m_tk->rv(yStageId).pdf()));
633  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
634  << ", rvYX.lawExpVector = " << pdfYX->lawExpVector()
635  << ", rvYX.lawVarVector = " << pdfYX->lawVarVector()
636  << ", rvYX.lawCovMatrix = " << pdfYX->lawCovMatrix()
637  << std::endl;
638  }
639  double qxy = m_tk->rv(xStageId).pdf().lnValue(y.vecValues(),NULL,NULL,NULL,NULL);
640  if ((m_env.subDisplayFile() ) &&
641  (m_env.displayVerbosity() >= 10 ) &&
642  (m_optionsObj->m_ov.m_totallyMute == false)) {
643  const InvLogitGaussianJointPdf<P_V,P_M>* pdfXY = dynamic_cast< const InvLogitGaussianJointPdf<P_V,P_M>* >(&(m_tk->rv(xStageId).pdf()));
644  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
645  << ", rvXY.lawExpVector = " << pdfXY->lawExpVector()
646  << ", rvXY.lawVarVector = " << pdfXY->lawVarVector()
647  << ", rvXY.lawCovMatrix = " << pdfXY->lawCovMatrix()
648  << std::endl;
649  }
650  alphaQuotient = std::exp(yLogTargetToUse +
651  qyx -
652  x.logTarget() -
653  qxy);
654  if ((m_env.subDisplayFile() ) &&
655  (m_env.displayVerbosity() >= 3 ) &&
656  (m_optionsObj->m_ov.m_totallyMute == false)) {
657  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
658  << ": asymmetric proposal case"
659  << ", xStageId = " << xStageId
660  << ", yStageId = " << yStageId
661  << ", x = " << x.vecValues()
662  << ", y = " << y.vecValues()
663  << ", yLogTargetToUse = " << yLogTargetToUse
664  << ", q(y,x) = " << qyx
665  << ", x.logTarget() = " << x.logTarget()
666  << ", q(x,y) = " << qxy
667  << ", alpha = " << alphaQuotient
668  << std::endl;
669  }
670  }
671  } // protection logic against logTarget values
672  }
673  else {
674  if ((m_env.subDisplayFile() ) &&
675  (m_env.displayVerbosity() >= 10 ) &&
676  (m_optionsObj->m_ov.m_totallyMute == false)) {
677  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
678  << ": x.outOfTargetSupport = " << x.outOfTargetSupport()
679  << ", y.outOfTargetSupport = " << y.outOfTargetSupport()
680  << std::endl;
681  }
682  }
683  if (alphaQuotientPtr != NULL) *alphaQuotientPtr = alphaQuotient;
684 
685  return std::min(1.,alphaQuotient);
686 }
687 //--------------------------------------------------
688 template<class P_V,class P_M>
689 double
690 MetropolisHastingsSG<P_V,P_M>::alpha(
691  const std::vector<MarkovChainPositionData<P_V>*>& inputPositionsData,
692  const std::vector<unsigned int >& inputTKStageIds)
693 {
694  unsigned int inputSize = inputPositionsData.size();
695  if ((m_env.subDisplayFile() ) &&
696  (m_env.displayVerbosity() >= 10 ) &&
697  (m_optionsObj->m_ov.m_totallyMute == false)) {
698  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
699  << ", inputSize = " << inputSize
700  << std::endl;
701  }
702  UQ_FATAL_TEST_MACRO((inputSize < 2),
703  m_env.worldRank(),
704  "MetropolisHastingsSG<P_V,P_M>::alpha(vec)",
705  "inputPositionsData has size < 2");
706 
707  // If necessary, return 0. right away
708  if (inputPositionsData[0 ]->outOfTargetSupport()) return 0.;
709  if (inputPositionsData[inputSize-1]->outOfTargetSupport()) return 0.;
710 
711  if ((inputPositionsData[0]->logTarget() == -INFINITY ) ||
712  (inputPositionsData[0]->logTarget() == INFINITY ) ||
713  ( (boost::math::isnan)(inputPositionsData[0]->logTarget()) )) {
714  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
715  << ", worldRank " << m_env.worldRank()
716  << ", fullRank " << m_env.fullRank()
717  << ", subEnvironment " << m_env.subId()
718  << ", subRank " << m_env.subRank()
719  << ", inter0Rank " << m_env.inter0Rank()
720  << ", positionId = " << m_positionIdForDebugging
721  << ", stageId = " << m_stageIdForDebugging
722  << ": inputSize = " << inputSize
723  << ", inputPositionsData[0]->logTarget() = " << inputPositionsData[0]->logTarget()
724  << ", [0]->values() = " << inputPositionsData[0]->vecValues()
725  << ", [inputSize - 1]->values() = " << inputPositionsData[inputSize-1]->vecValues()
726  << std::endl;
727  return 0.;
728  }
729  else if ((inputPositionsData[inputSize - 1]->logTarget() == -INFINITY ) ||
730  (inputPositionsData[inputSize - 1]->logTarget() == INFINITY ) ||
731  ( (boost::math::isnan)(inputPositionsData[inputSize - 1]->logTarget()) )) {
732  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
733  << ", worldRank " << m_env.worldRank()
734  << ", fullRank " << m_env.fullRank()
735  << ", subEnvironment " << m_env.subId()
736  << ", subRank " << m_env.subRank()
737  << ", inter0Rank " << m_env.inter0Rank()
738  << ", positionId = " << m_positionIdForDebugging
739  << ", stageId = " << m_stageIdForDebugging
740  << ": inputSize = " << inputSize
741  << ", inputPositionsData[inputSize - 1]->logTarget() = " << inputPositionsData[inputSize-1]->logTarget()
742  << ", [0]->values() = " << inputPositionsData[0]->vecValues()
743  << ", [inputSize - 1]->values() = " << inputPositionsData[inputSize-1]->vecValues()
744  << std::endl;
745  return 0.;
746  }
747 
748  // If inputSize is 2, recursion is not needed
749  if (inputSize == 2) return this->alpha(*(inputPositionsData[0 ]),
750  *(inputPositionsData[inputSize - 1]),
751  inputTKStageIds[0],
752  inputTKStageIds[inputSize-1]);
753 
754  // Prepare two vectors of positions
755  std::vector<MarkovChainPositionData<P_V>*> positionsData (inputSize,NULL);
756  std::vector<MarkovChainPositionData<P_V>*> backwardPositionsData (inputSize,NULL);
757 
758  std::vector<unsigned int > tkStageIds (inputSize,0);
759  std::vector<unsigned int > backwardTKStageIds (inputSize,0);
760 
761  std::vector<unsigned int > tkStageIdsLess1(inputSize,0);
762  std::vector<unsigned int > backwardTKStageIdsLess1(inputSize,0);
763 
764  for (unsigned int i = 0; i < inputSize; ++i) {
765  positionsData [i] = inputPositionsData[i];
766  backwardPositionsData [i] = inputPositionsData[inputSize-i-1];
767 
768  tkStageIds [i] = inputTKStageIds [i];
769  backwardTKStageIds [i] = inputTKStageIds [inputSize-i-1];
770 
771  tkStageIdsLess1[i] = inputTKStageIds [i];
772  backwardTKStageIdsLess1[i] = inputTKStageIds [inputSize-i-1];
773  }
774 
775  tkStageIdsLess1.pop_back();
776  backwardTKStageIdsLess1.pop_back();
777 
778  // Initialize cumulative variables
779  double logNumerator = 0.;
780  double logDenominator = 0.;
781  double alphasNumerator = 1.;
782  double alphasDenominator = 1.;
783 
784  // Compute cumulative variables
785  const P_V& _lastTKPosition = m_tk->preComputingPosition( tkStageIds[inputSize-1]);
786  const P_V& _lastBackwardTKPosition = m_tk->preComputingPosition(backwardTKStageIds[inputSize-1]);
787 
788  double numContrib = m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(_lastBackwardTKPosition,NULL,NULL,NULL,NULL);
789  double denContrib = m_tk->rv(tkStageIdsLess1).pdf().lnValue(_lastTKPosition,NULL,NULL,NULL,NULL);
790  if ((m_env.subDisplayFile() ) &&
791  (m_env.displayVerbosity() >= 10 ) &&
792  (m_optionsObj->m_ov.m_totallyMute == false)) {
793  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
794  << ", inputSize = " << inputSize
795  << ", before loop"
796  << ": numContrib = " << numContrib
797  << ", denContrib = " << denContrib
798  << std::endl;
799  }
800  logNumerator += numContrib;
801  logDenominator += denContrib;
802 
803  for (unsigned int i = 0; i < (inputSize-2); ++i) { // That is why size must be >= 2
804  positionsData.pop_back();
805  backwardPositionsData.pop_back();
806 
807  const P_V& lastTKPosition = m_tk->preComputingPosition( tkStageIds[inputSize-2-i]);
808  const P_V& lastBackwardTKPosition = m_tk->preComputingPosition(backwardTKStageIds[inputSize-2-i]);
809 
810  tkStageIds.pop_back();
811  backwardTKStageIds.pop_back();
812 
813  tkStageIdsLess1.pop_back();
814  backwardTKStageIdsLess1.pop_back();
815 
816  numContrib = m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(lastBackwardTKPosition,NULL,NULL,NULL,NULL);
817  denContrib = m_tk->rv(tkStageIdsLess1).pdf().lnValue(lastTKPosition,NULL,NULL,NULL,NULL);
818  if ((m_env.subDisplayFile() ) &&
819  (m_env.displayVerbosity() >= 10 ) &&
820  (m_optionsObj->m_ov.m_totallyMute == false)) {
821  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
822  << ", inputSize = " << inputSize
823  << ", in loop, i = " << i
824  << ": numContrib = " << numContrib
825  << ", denContrib = " << denContrib
826  << std::endl;
827  }
828  logNumerator += numContrib;
829  logDenominator += denContrib;
830 
831  alphasNumerator *= (1. - this->alpha(backwardPositionsData,backwardTKStageIds));
832  alphasDenominator *= (1. - this->alpha( positionsData, tkStageIds));
833  }
834 
835  double numeratorLogTargetToUse = backwardPositionsData[0]->logTarget();
836  numContrib = numeratorLogTargetToUse;
837  denContrib = positionsData[0]->logTarget();
838  if ((m_env.subDisplayFile() ) &&
839  (m_env.displayVerbosity() >= 10 ) &&
840  (m_optionsObj->m_ov.m_totallyMute == false)) {
841  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
842  << ", inputSize = " << inputSize
843  << ", after loop"
844  << ": numContrib = " << numContrib
845  << ", denContrib = " << denContrib
846  << std::endl;
847  }
848  logNumerator += numContrib;
849  logDenominator += denContrib;
850 
851  if ((m_env.subDisplayFile() ) &&
852  (m_env.displayVerbosity() >= 10 ) &&
853  (m_optionsObj->m_ov.m_totallyMute == false)) {
854  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
855  << ", inputSize = " << inputSize
856  << ": alphasNumerator = " << alphasNumerator
857  << ", alphasDenominator = " << alphasDenominator
858  << ", logNumerator = " << logNumerator
859  << ", logDenominator = " << logDenominator
860  << std::endl;
861  }
862 
863  // Return result
864  return std::min(1.,(alphasNumerator/alphasDenominator)*std::exp(logNumerator-logDenominator));
865 }
866 //--------------------------------------------------
867 template<class P_V,class P_M>
868 bool
869 MetropolisHastingsSG<P_V,P_M>::acceptAlpha(double alpha)
870 {
871  bool result = false;
872 
873  if (alpha <= 0. ) result = false;
874  else if (alpha >= 1. ) result = true;
875  else if (alpha >= m_env.rngObject()->uniformSample()) result = true;
876  else result = false;
877 
878  return result;
879 }
880 //--------------------------------------------------
881 template<class P_V,class P_M>
882 int
883 MetropolisHastingsSG<P_V,P_M>::writeInfo(
884  const BaseVectorSequence<P_V,P_M>& workingChain,
885  std::ofstream& ofsvar) const
886 {
887  if ((m_env.subDisplayFile() ) &&
888  (m_optionsObj->m_ov.m_totallyMute == false)) {
889  *m_env.subDisplayFile() << "\n"
890  << "\n-----------------------------------------------------"
891  << "\n Writing more information about the Markov chain " << workingChain.name() << " to output file ..."
892  << "\n-----------------------------------------------------"
893  << "\n"
894  << std::endl;
895  }
896 
897  int iRC = UQ_OK_RC;
898 
899  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
900  // Write m_logTargets
901  ofsvar << m_optionsObj->m_prefix << "logTargets_sub" << m_env.subIdString() << " = zeros(" << m_logTargets.size()
902  << "," << 1
903  << ");"
904  << std::endl;
905  ofsvar << m_optionsObj->m_prefix << "logTargets_sub" << m_env.subIdString() << " = [";
906  for (unsigned int i = 0; i < m_logTargets.size(); ++i) {
907  ofsvar << m_logTargets[i]
908  << std::endl;
909  }
910  ofsvar << "];\n";
911 
912  // Write m_alphaQuotients
913  ofsvar << m_optionsObj->m_prefix << "alphaQuotients_sub" << m_env.subIdString() << " = zeros(" << m_alphaQuotients.size()
914  << "," << 1
915  << ");"
916  << std::endl;
917  ofsvar << m_optionsObj->m_prefix << "alphaQuotients_sub" << m_env.subIdString() << " = [";
918  for (unsigned int i = 0; i < m_alphaQuotients.size(); ++i) {
919  ofsvar << m_alphaQuotients[i]
920  << std::endl;
921  }
922  ofsvar << "];\n";
923  }
924 
925  // Write number of rejections
926  ofsvar << m_optionsObj->m_prefix << "rejected = " << (double) m_rawChainInfo.numRejections/(double) (workingChain.subSequenceSize()-1)
927  << ";\n"
928  << std::endl;
929 
930  if (false) { // Don't see need for code below. Let it there though, compiling, in case it is needed in the future.
931  // Write names of components
932  ofsvar << m_optionsObj->m_prefix << "componentNames = {";
933  m_vectorSpace.printComponentsNames(ofsvar,false);
934  ofsvar << "};\n";
935 
936  // Write number of out of target support
937  ofsvar << m_optionsObj->m_prefix << "outTargetSupport = " << (double) m_rawChainInfo.numOutOfTargetSupport/(double) (workingChain.subSequenceSize()-1)
938  << ";\n"
939  << std::endl;
940 
941  // Write chain run time
942  ofsvar << m_optionsObj->m_prefix << "runTime = " << m_rawChainInfo.runTime
943  << ";\n"
944  << std::endl;
945  }
946 
947  if ((m_env.subDisplayFile() ) &&
948  (m_optionsObj->m_ov.m_totallyMute == false)) {
949  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
950  << "\n Finished writing more information about the Markov chain " << workingChain.name()
951  << "\n-----------------------------------------------------"
952  << "\n"
953  << std::endl;
954  }
955 
956  return iRC;
957 }
958 
959 template <class P_V, class P_M>
960 const BaseTKGroup<P_V, P_M> &
961 MetropolisHastingsSG<P_V, P_M>::transitionKernel() const
962 {
963  return *m_tk;
964 }
965 
966 // Statistical methods -----------------------------
967 /* This operation currently implements the DRAM algorithm (Heikki Haario, Marko
968  * Laine, Antonietta Mira and Eero Saksman, "DRAM: Efficient Adaptive MCMC",
969  * Statistics and Computing (2006), 16:339-354). It also provides support for
970  * Stochastic Newton algorithm through the TK (transition kernel) class. Stochastic
971  * Newton is not totally implemented yet though, since it is being researched by
972  * James Martin and Omar Ghattas at ICES at the University of Texas at Austin.*/
973 template <class P_V,class P_M>
974 void
975 MetropolisHastingsSG<P_V,P_M>::generateSequence(
976  BaseVectorSequence<P_V,P_M>& workingChain,
977  ScalarSequence<double>* workingLogLikelihoodValues, // KEY: add LogPriorValues
978  ScalarSequence<double>* workingLogTargetValues)
979 {
980  if ((m_env.subDisplayFile() ) &&
981  (m_env.displayVerbosity() >= 5 ) &&
982  (m_optionsObj->m_ov.m_totallyMute == false)) {
983  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::generateSequence()..."
984  << std::endl;
985  }
986 
987  if (m_vectorSpace.dimLocal() != workingChain.vectorSizeLocal()) {
988  std::cerr << "'m_vectorSpace' and 'workingChain' are related to vector"
989  << "spaces of different dimensions"
990  << std::endl;
991  queso_error();
992  }
993 
994  // Set a flag to write out log likelihood or not
995  bool writeLogLikelihood;
996  if ((workingLogLikelihoodValues != NULL) &&
997  (m_optionsObj->m_ov.m_outputLogLikelihood)) {
998  writeLogLikelihood = true;
999  }
1000  else {
1001  writeLogLikelihood = false;
1002  }
1003 
1004  // Set a flag to write out log target or not
1005  bool writeLogTarget;
1006  if ((workingLogTargetValues != NULL) &&
1007  (m_optionsObj->m_ov.m_outputLogTarget)) {
1008  writeLogTarget = true;
1009  }
1010  else {
1011  writeLogTarget = false;
1012  }
1013 
1014  MiscCheckTheParallelEnvironment<P_V,P_V>(m_initialPosition,
1015  m_initialPosition);
1016 
1017  P_V valuesOf1stPosition(m_initialPosition);
1018  int iRC = UQ_OK_RC;
1019 
1020  workingChain.setName(m_optionsObj->m_prefix + "rawChain");
1021 
1022  //****************************************************
1023  // Generate chain
1024  //****************************************************
1025  if (m_optionsObj->m_ov.m_rawChainDataInputFileName == UQ_MH_SG_FILENAME_FOR_NO_FILE) {
1026  generateFullChain(valuesOf1stPosition,
1027  m_optionsObj->m_ov.m_rawChainSize,
1028  workingChain,
1029  workingLogLikelihoodValues,
1030  workingLogTargetValues);
1031  }
1032  else {
1033  readFullChain(m_optionsObj->m_ov.m_rawChainDataInputFileName,
1034  m_optionsObj->m_ov.m_rawChainDataInputFileType,
1035  m_optionsObj->m_ov.m_rawChainSize,
1036  workingChain);
1037  }
1038 
1039  //****************************************************
1040  // Open generic output file
1041  //****************************************************
1042  if ((m_env.subDisplayFile() ) &&
1043  (m_optionsObj->m_ov.m_totallyMute == false)) {
1044  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1045  << ", prefix = " << m_optionsObj->m_prefix
1046  << ", chain name = " << workingChain.name()
1047  << ": about to try to open generic output file '" << m_optionsObj->m_ov.m_dataOutputFileName
1048  << "." << UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT // Yes, always ".m"
1049  << "', subId = " << m_env.subId()
1050  << ", subenv is allowed to write (1/true or 0/false) = " << (m_optionsObj->m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_ov.m_dataOutputAllowedSet.end())
1051  << "..."
1052  << std::endl;
1053  }
1054 
1055  FilePtrSetStruct genericFilePtrSet;
1056  m_env.openOutputFile(m_optionsObj->m_ov.m_dataOutputFileName,
1057  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
1058  m_optionsObj->m_ov.m_dataOutputAllowedSet,
1059  false,
1060  genericFilePtrSet);
1061 
1062  if ((m_env.subDisplayFile() ) &&
1063  (m_optionsObj->m_ov.m_totallyMute == false)) {
1064  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1065  << ", prefix = " << m_optionsObj->m_prefix
1066  << ", raw chain name = " << workingChain.name()
1067  << ": returned from opening generic output file '" << m_optionsObj->m_ov.m_dataOutputFileName
1068  << "." << UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT // Yes, always ".m"
1069  << "', subId = " << m_env.subId()
1070  << std::endl;
1071  }
1072 
1073  //****************************************************************************************
1074  // Eventually:
1075  // --> write raw chain
1076  // --> compute statistics on it
1077  //****************************************************************************************
1078  if ((m_optionsObj->m_ov.m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) &&
1079  (m_optionsObj->m_ov.m_totallyMute == false )) {
1080 
1081  // Take "sub" care of raw chain
1082  if ((m_env.subDisplayFile() ) &&
1083  (m_optionsObj->m_ov.m_totallyMute == false)) {
1084  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1085  << ", prefix = " << m_optionsObj->m_prefix
1086  << ", raw chain name = " << workingChain.name()
1087  << ": about to try to write raw sub chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
1088  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
1089  << "', subId = " << m_env.subId()
1090  << ", subenv is allowed to write 1/true or 0/false) = " << (m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet.end())
1091  << "..."
1092  << std::endl;
1093  }
1094 
1095  if ((m_numPositionsNotSubWritten > 0 ) &&
1096  (m_optionsObj->m_ov.m_rawChainDataOutputFileName != ".")) {
1097  workingChain.subWriteContents(m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten,
1098  m_numPositionsNotSubWritten,
1099  m_optionsObj->m_ov.m_rawChainDataOutputFileName,
1100  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1101  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1102  if ((m_env.subDisplayFile() ) &&
1103  (m_optionsObj->m_ov.m_totallyMute == false)) {
1104  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1105  << ": just wrote (per period request) remaining " << m_numPositionsNotSubWritten << " chain positions "
1106  << ", " << m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten << " <= pos <= " << m_optionsObj->m_ov.m_rawChainSize - 1
1107  << std::endl;
1108  }
1109 
1110  if (writeLogLikelihood) {
1111  workingLogLikelihoodValues->subWriteContents(m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten,
1112  m_numPositionsNotSubWritten,
1113  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_loglikelihood",
1114  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1115  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1116  }
1117 
1118  if (writeLogTarget) {
1119  workingLogTargetValues->subWriteContents(m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten,
1120  m_numPositionsNotSubWritten,
1121  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_logtarget",
1122  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1123  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1124  }
1125 
1126  m_numPositionsNotSubWritten = 0;
1127  }
1128 
1129  // Compute raw sub MLE
1130  if (workingLogLikelihoodValues) {
1131  SequenceOfVectors<P_V,P_M> rawSubMLEpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawSubMLEseq");
1132  double rawSubMLEvalue = workingChain.subPositionsOfMaximum(*workingLogLikelihoodValues,
1133  rawSubMLEpositions);
1134  UQ_FATAL_TEST_MACRO(rawSubMLEpositions.subSequenceSize() == 0,
1135  m_env.worldRank(),
1136  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
1137  "rawSubMLEpositions.subSequenceSize() = 0");
1138 
1139  if ((m_env.subDisplayFile() ) &&
1140  (m_optionsObj->m_ov.m_totallyMute == false)) {
1141  P_V tmpVec(m_vectorSpace.zeroVector());
1142  rawSubMLEpositions.getPositionValues(0,tmpVec);
1143  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1144  << ": just computed MLE"
1145  << ", rawSubMLEvalue = " << rawSubMLEvalue
1146  << ", rawSubMLEpositions.subSequenceSize() = " << rawSubMLEpositions.subSequenceSize()
1147  << ", rawSubMLEpositions[0] = " << tmpVec
1148  << std::endl;
1149  }
1150  }
1151 
1152  // Compute raw sub MAP
1153  if (workingLogTargetValues) {
1154  SequenceOfVectors<P_V,P_M> rawSubMAPpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawSubMAPseq");
1155  double rawSubMAPvalue = workingChain.subPositionsOfMaximum(*workingLogTargetValues,
1156  rawSubMAPpositions);
1157  UQ_FATAL_TEST_MACRO(rawSubMAPpositions.subSequenceSize() == 0,
1158  m_env.worldRank(),
1159  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
1160  "rawSubMAPpositions.subSequenceSize() = 0");
1161 
1162  if ((m_env.subDisplayFile() ) &&
1163  (m_optionsObj->m_ov.m_totallyMute == false)) {
1164  P_V tmpVec(m_vectorSpace.zeroVector());
1165  rawSubMAPpositions.getPositionValues(0,tmpVec);
1166  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1167  << ": just computed MAP"
1168  << ", rawSubMAPvalue = " << rawSubMAPvalue
1169  << ", rawSubMAPpositions.subSequenceSize() = " << rawSubMAPpositions.subSequenceSize()
1170  << ", rawSubMAPpositions[0] = " << tmpVec
1171  << std::endl;
1172  }
1173  }
1174 
1175  if ((m_env.subDisplayFile() ) &&
1176  (m_optionsObj->m_ov.m_totallyMute == false)) {
1177  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1178  << ", prefix = " << m_optionsObj->m_prefix
1179  << ", raw chain name = " << workingChain.name()
1180  << ": returned from writing raw sub chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
1181  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
1182  << "', subId = " << m_env.subId()
1183  << std::endl;
1184  }
1185 
1186  // Take "unified" care of raw chain
1187  if ((m_env.subDisplayFile() ) &&
1188  (m_optionsObj->m_ov.m_totallyMute == false)) {
1189  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1190  << ", prefix = " << m_optionsObj->m_prefix
1191  << ", raw chain name = " << workingChain.name()
1192  << ": about to try to write raw unified chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
1193  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
1194  << "', subId = " << m_env.subId()
1195  << "..."
1196  << std::endl;
1197  }
1198 
1199  workingChain.unifiedWriteContents(m_optionsObj->m_ov.m_rawChainDataOutputFileName,
1200  m_optionsObj->m_ov.m_rawChainDataOutputFileType);
1201  if ((m_env.subDisplayFile() ) &&
1202  (m_optionsObj->m_ov.m_totallyMute == false)) {
1203  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1204  << ", prefix = " << m_optionsObj->m_prefix
1205  << ", raw chain name = " << workingChain.name()
1206  << ": returned from writing raw unified chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
1207  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
1208  << "', subId = " << m_env.subId()
1209  << std::endl;
1210  }
1211 
1212  if (writeLogLikelihood) {
1213  workingLogLikelihoodValues->unifiedWriteContents(m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_loglikelihood",
1214  m_optionsObj->m_ov.m_rawChainDataOutputFileType);
1215  }
1216 
1217  if (writeLogTarget) {
1218  workingLogTargetValues->unifiedWriteContents(m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_logtarget",
1219  m_optionsObj->m_ov.m_rawChainDataOutputFileType);
1220  }
1221 
1222  // Compute raw unified MLE
1223  if (workingLogLikelihoodValues) {
1224  SequenceOfVectors<P_V,P_M> rawUnifiedMLEpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawUnifiedMLEseq");
1225 
1226  double rawUnifiedMLEvalue = workingChain.unifiedPositionsOfMaximum(*workingLogLikelihoodValues,
1227  rawUnifiedMLEpositions);
1228 
1229  if ((m_env.subDisplayFile() ) &&
1230  (m_optionsObj->m_ov.m_totallyMute == false)) {
1231  P_V tmpVec(m_vectorSpace.zeroVector());
1232 
1233  // Make sure the positions vector (which only contains stuff on process
1234  // zero) actually contains positions
1235  if (rawUnifiedMLEpositions.subSequenceSize() > 0) {
1236  rawUnifiedMLEpositions.getPositionValues(0,tmpVec);
1237  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1238  << ": just computed MLE"
1239  << ", rawUnifiedMLEvalue = " << rawUnifiedMLEvalue
1240  << ", rawUnifiedMLEpositions.subSequenceSize() = " << rawUnifiedMLEpositions.subSequenceSize()
1241  << ", rawUnifiedMLEpositions[0] = " << tmpVec
1242  << std::endl;
1243  }
1244  }
1245  }
1246 
1247  // Compute raw unified MAP
1248  if (workingLogTargetValues) {
1249  SequenceOfVectors<P_V,P_M> rawUnifiedMAPpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawUnifiedMAPseq");
1250  double rawUnifiedMAPvalue = workingChain.unifiedPositionsOfMaximum(*workingLogTargetValues,
1251  rawUnifiedMAPpositions);
1252 
1253  if ((m_env.subDisplayFile() ) &&
1254  (m_optionsObj->m_ov.m_totallyMute == false)) {
1255  P_V tmpVec(m_vectorSpace.zeroVector());
1256 
1257  // Make sure the positions vector (which only contains stuff on process
1258  // zero) actually contains positions
1259  if (rawUnifiedMAPpositions.subSequenceSize() > 0) {
1260  rawUnifiedMAPpositions.getPositionValues(0,tmpVec);
1261  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1262  << ": just computed MAP"
1263  << ", rawUnifiedMAPvalue = " << rawUnifiedMAPvalue
1264  << ", rawUnifiedMAPpositions.subSequenceSize() = " << rawUnifiedMAPpositions.subSequenceSize()
1265  << ", rawUnifiedMAPpositions[0] = " << tmpVec
1266  << std::endl;
1267  }
1268  }
1269  }
1270  }
1271 
1272  // Take care of other aspects of raw chain
1273  if ((genericFilePtrSet.ofsVar ) &&
1274  (m_optionsObj->m_ov.m_totallyMute == false)) {
1275  // Write likelihoodValues and alphaValues, if they were requested by user
1276  iRC = writeInfo(workingChain,
1277  *genericFilePtrSet.ofsVar);
1278  UQ_FATAL_RC_MACRO(iRC,
1279  m_env.worldRank(),
1280  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
1281  "improper writeInfo() return");
1282  }
1283 
1284 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
1285  if (m_optionsObj->m_ov.m_rawChainComputeStats) {
1286  workingChain.computeStatistics(*m_optionsObj->m_rawChainStatisticalOptionsObj,
1287  genericFilePtrSet.ofsVar);
1288  }
1289 #endif
1290 
1291  //****************************************************************************************
1292  // Eventually:
1293  // --> filter the raw chain
1294  // --> write it
1295  // --> compute statistics on it
1296  //****************************************************************************************
1297  if (m_optionsObj->m_ov.m_filteredChainGenerate) {
1298  // Compute filter parameters
1299  unsigned int filterInitialPos = (unsigned int) (m_optionsObj->m_ov.m_filteredChainDiscardedPortion * (double) workingChain.subSequenceSize());
1300  unsigned int filterSpacing = m_optionsObj->m_ov.m_filteredChainLag;
1301  if (filterSpacing == 0) {
1302  workingChain.computeFilterParams(genericFilePtrSet.ofsVar,
1303  filterInitialPos,
1304  filterSpacing);
1305  }
1306 
1307  // Filter positions from the converged portion of the chain
1308  workingChain.filter(filterInitialPos,
1309  filterSpacing);
1310  workingChain.setName(m_optionsObj->m_prefix + "filtChain");
1311 
1312  if (workingLogLikelihoodValues) workingLogLikelihoodValues->filter(filterInitialPos,
1313  filterSpacing);
1314 
1315  if (workingLogTargetValues) workingLogTargetValues->filter(filterInitialPos,
1316  filterSpacing);
1317 
1318  // Write filtered chain
1319  if ((m_env.subDisplayFile() ) &&
1320  (m_optionsObj->m_ov.m_totallyMute == false)) {
1321  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1322  << ", prefix = " << m_optionsObj->m_prefix
1323  << ": checking necessity of opening output files for filtered chain " << workingChain.name()
1324  << "..."
1325  << std::endl;
1326  }
1327 
1328  // Take "sub" care of filtered chain
1329  if ((m_optionsObj->m_ov.m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) &&
1330  (m_optionsObj->m_ov.m_totallyMute == false )) {
1331  workingChain.subWriteContents(0,
1332  workingChain.subSequenceSize(),
1333  m_optionsObj->m_ov.m_filteredChainDataOutputFileName,
1334  m_optionsObj->m_ov.m_filteredChainDataOutputFileType,
1335  m_optionsObj->m_ov.m_filteredChainDataOutputAllowedSet);
1336  if ((m_env.subDisplayFile() ) &&
1337  (m_optionsObj->m_ov.m_totallyMute == false)) {
1338  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1339  << ", prefix = " << m_optionsObj->m_prefix
1340  << ": closed sub output file '" << m_optionsObj->m_ov.m_filteredChainDataOutputFileName
1341  << "' for filtered chain " << workingChain.name()
1342  << std::endl;
1343  }
1344 
1345  if (writeLogLikelihood) {
1346  workingLogLikelihoodValues->subWriteContents(0,
1347  workingChain.subSequenceSize(),
1348  m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_loglikelihood",
1349  m_optionsObj->m_ov.m_filteredChainDataOutputFileType,
1350  m_optionsObj->m_ov.m_filteredChainDataOutputAllowedSet);
1351  }
1352 
1353  if (writeLogTarget) {
1354  workingLogTargetValues->subWriteContents(0,
1355  workingChain.subSequenceSize(),
1356  m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_logtarget",
1357  m_optionsObj->m_ov.m_filteredChainDataOutputFileType,
1358  m_optionsObj->m_ov.m_filteredChainDataOutputAllowedSet);
1359  }
1360  }
1361 
1362  // Compute sub filtered MLE and sub filtered MAP
1363 
1364  // Take "unified" care of filtered chain
1365  if ((m_optionsObj->m_ov.m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) &&
1366  (m_optionsObj->m_ov.m_totallyMute == false )) {
1367  workingChain.unifiedWriteContents(m_optionsObj->m_ov.m_filteredChainDataOutputFileName,
1368  m_optionsObj->m_ov.m_filteredChainDataOutputFileType);
1369  if ((m_env.subDisplayFile() ) &&
1370  (m_optionsObj->m_ov.m_totallyMute == false)) {
1371  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1372  << ", prefix = " << m_optionsObj->m_prefix
1373  << ": closed unified output file '" << m_optionsObj->m_ov.m_filteredChainDataOutputFileName
1374  << "' for filtered chain " << workingChain.name()
1375  << std::endl;
1376  }
1377 
1378  if (writeLogLikelihood) {
1379  workingLogLikelihoodValues->unifiedWriteContents(m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_loglikelihood",
1380  m_optionsObj->m_ov.m_filteredChainDataOutputFileType);
1381  }
1382 
1383  if (writeLogTarget) {
1384  workingLogTargetValues->unifiedWriteContents(m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_logtarget",
1385  m_optionsObj->m_ov.m_filteredChainDataOutputFileType);
1386  }
1387  }
1388 
1389  // Compute unified filtered MLE and unified filtered MAP
1390 
1391  // Compute statistics
1392 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
1393  if (m_optionsObj->m_ov.m_filteredChainComputeStats) {
1394  workingChain.computeStatistics(*m_optionsObj->m_filteredChainStatisticalOptionsObj,
1395  genericFilePtrSet.ofsVar);
1396  }
1397 #endif
1398  }
1399 
1400  //****************************************************
1401  // Close generic output file
1402  //****************************************************
1403  if (genericFilePtrSet.ofsVar) {
1404  //genericFilePtrSet.ofsVar->close();
1405  delete genericFilePtrSet.ofsVar;
1406  if ((m_env.subDisplayFile() ) &&
1407  (m_optionsObj->m_ov.m_totallyMute == false)) {
1408  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1409  << ", prefix = " << m_optionsObj->m_prefix
1410  << ": closed generic output file '" << m_optionsObj->m_ov.m_dataOutputFileName
1411  << "' (chain name is " << workingChain.name()
1412  << ")"
1413  << std::endl;
1414  }
1415  }
1416 
1417  if ((m_env.subDisplayFile() ) &&
1418  (m_optionsObj->m_ov.m_totallyMute == false)) {
1419  *m_env.subDisplayFile() << std::endl;
1420  }
1421 
1422  //m_env.syncPrintDebugMsg("Leaving MetropolisHastingsSG<P_V,P_M>::generateSequence()",2,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
1423 
1424  if ((m_env.subDisplayFile() ) &&
1425  (m_env.displayVerbosity() >= 5 ) &&
1426  (m_optionsObj->m_ov.m_totallyMute == false)) {
1427  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1428  << std::endl;
1429  }
1430 
1431  return;
1432 }
1433 
1434 // -------------------------------------------------
1435 template<class P_V,class P_M>
1436 void
1437 MetropolisHastingsSG<P_V,P_M>::getRawChainInfo(MHRawChainInfoStruct& info) const
1438 {
1439  info = m_rawChainInfo;
1440  return;
1441 }
1442 //--------------------------------------------------
1443 template <class P_V,class P_M>
1444 void
1445 MetropolisHastingsSG<P_V,P_M>::readFullChain(
1446  const std::string& inputFileName,
1447  const std::string& inputFileType,
1448  unsigned int chainSize,
1449  BaseVectorSequence<P_V,P_M>& workingChain)
1450 {
1451  workingChain.unifiedReadContents(inputFileName,inputFileType,chainSize);
1452  return;
1453 }
1454 // Private methods ---------------------------------
1455 template <class P_V,class P_M>
1456 void
1457 MetropolisHastingsSG<P_V,P_M>::generateFullChain(
1458  const P_V& valuesOf1stPosition,
1459  unsigned int chainSize,
1460  BaseVectorSequence<P_V,P_M>& workingChain,
1461  ScalarSequence<double>* workingLogLikelihoodValues,
1462  ScalarSequence<double>* workingLogTargetValues)
1463 {
1464  //m_env.syncPrintDebugMsg("Entering MetropolisHastingsSG<P_V,P_M>::generateFullChain()",3,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
1465 
1466  if ((m_env.subDisplayFile() ) &&
1467  (m_optionsObj->m_ov.m_totallyMute == false)) {
1468  *m_env.subDisplayFile() << "Starting the generation of Markov chain " << workingChain.name()
1469  << ", with " << chainSize
1470  << " positions..."
1471  << std::endl;
1472  }
1473 
1474  int iRC = UQ_OK_RC;
1475  struct timeval timevalChain;
1476  struct timeval timevalCandidate;
1477  struct timeval timevalTarget;
1478  struct timeval timevalMhAlpha;
1479  struct timeval timevalDrAlpha;
1480  struct timeval timevalDR;
1481  struct timeval timevalAM;
1482 
1483  m_positionIdForDebugging = 0;
1484  m_stageIdForDebugging = 0;
1485 
1486  m_rawChainInfo.reset();
1487 
1488  iRC = gettimeofday(&timevalChain, NULL);
1489 
1490  if ((m_env.subDisplayFile() ) &&
1491  (m_optionsObj->m_ov.m_totallyMute == false)) {
1492  *m_env.subDisplayFile() << "\nIn MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1493  << ": contents of initial position are:";
1494  *m_env.subDisplayFile() << valuesOf1stPosition; // FIX ME: might need parallelism
1495  *m_env.subDisplayFile() << "\nIn MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1496  << ": targetPdf.domaintSet() info is:"
1497  << m_targetPdf.domainSet();
1498  *m_env.subDisplayFile() << std::endl;
1499  }
1500 
1501  // Set a flag to write out log likelihood or not
1502  bool writeLogLikelihood;
1503  if ((workingLogLikelihoodValues != NULL) &&
1504  (m_optionsObj->m_ov.m_outputLogLikelihood)) {
1505  writeLogLikelihood = true;
1506  }
1507  else {
1508  writeLogLikelihood = false;
1509  }
1510 
1511  // Set a flag to write out log target or not
1512  bool writeLogTarget;
1513  if ((workingLogTargetValues != NULL) &&
1514  (m_optionsObj->m_ov.m_outputLogTarget)) {
1515  writeLogTarget = true;
1516  }
1517  else {
1518  writeLogTarget = false;
1519  }
1520 
1521  bool outOfTargetSupport = !m_targetPdf.domainSet().contains(valuesOf1stPosition);
1522  if ((m_env.subDisplayFile()) &&
1523  (outOfTargetSupport )) {
1524  *m_env.subDisplayFile() << "ERROR: In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1525  << ": contents of initial position are:\n";
1526  *m_env.subDisplayFile() << valuesOf1stPosition; // FIX ME: might need parallelism
1527  *m_env.subDisplayFile() << "\nERROR: In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1528  << ": targetPdf.domaintSet() info is:\n"
1529  << m_targetPdf.domainSet();
1530  *m_env.subDisplayFile() << std::endl;
1531  }
1532  UQ_FATAL_TEST_MACRO(outOfTargetSupport,
1533  m_env.worldRank(),
1534  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
1535  "initial position should not be out of target pdf support");
1536 
1537  double logPrior = 0.;
1538  double logLikelihood = 0.;
1539  double logTarget = 0.;
1540  if (m_computeInitialPriorAndLikelihoodValues) {
1541  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalTarget, NULL);
1542  logTarget = m_targetPdfSynchronizer->callFunction(&valuesOf1stPosition,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment // KEY
1543  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) {
1544  m_rawChainInfo.targetRunTime += MiscGetEllapsedSeconds(&timevalTarget);
1545  }
1546  m_rawChainInfo.numTargetCalls++;
1547  if ((m_env.subDisplayFile() ) &&
1548  (m_env.displayVerbosity() >= 3 ) &&
1549  (m_optionsObj->m_ov.m_totallyMute == false)) {
1550  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1551  << ": just returned from likelihood() for initial chain position"
1552  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1553  << ", logPrior = " << logPrior
1554  << ", logLikelihood = " << logLikelihood
1555  << ", logTarget = " << logTarget
1556  << std::endl;
1557  }
1558  }
1559  else {
1560  logPrior = m_initialLogPriorValue;
1561  logLikelihood = m_initialLogLikelihoodValue;
1562  logTarget = logPrior + logLikelihood;
1563  if ((m_env.subDisplayFile() ) &&
1564  (m_env.displayVerbosity() >= 3 ) &&
1565  (m_optionsObj->m_ov.m_totallyMute == false)) {
1566  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1567  << ": used input prior and likelihood for initial chain position"
1568  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1569  << ", logPrior = " << logPrior
1570  << ", logLikelihood = " << logLikelihood
1571  << ", logTarget = " << logTarget
1572  << std::endl;
1573  }
1574  }
1575 
1576  //*m_env.subDisplayFile() << "AQUI 001" << std::endl;
1577  MarkovChainPositionData<P_V> currentPositionData(m_env,
1578  valuesOf1stPosition,
1579  outOfTargetSupport,
1580  logLikelihood,
1581  logTarget);
1582 
1583  P_V gaussianVector(m_vectorSpace.zeroVector());
1584  P_V tmpVecValues(m_vectorSpace.zeroVector());
1585  MarkovChainPositionData<P_V> currentCandidateData(m_env);
1586 
1587  //****************************************************
1588  // Set chain position with positionId = 0
1589  //****************************************************
1590  workingChain.resizeSequence(chainSize);
1591  m_numPositionsNotSubWritten = 0;
1592  if (workingLogLikelihoodValues) workingLogLikelihoodValues->resizeSequence(chainSize);
1593  if (workingLogTargetValues ) workingLogTargetValues->resizeSequence (chainSize);
1594  if (true/*m_uniqueChainGenerate*/) m_idsOfUniquePositions.resize(chainSize,0);
1595  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1596  m_logTargets.resize (chainSize,0.);
1597  m_alphaQuotients.resize(chainSize,0.);
1598  }
1599 
1600  unsigned int uniquePos = 0;
1601  workingChain.setPositionValues(0,currentPositionData.vecValues());
1602  m_numPositionsNotSubWritten++;
1603  if ((m_optionsObj->m_ov.m_rawChainDataOutputPeriod > 0 ) &&
1604  (((0+1) % m_optionsObj->m_ov.m_rawChainDataOutputPeriod) == 0 ) &&
1605  (m_optionsObj->m_ov.m_rawChainDataOutputFileName != ".")) {
1606  workingChain.subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1607  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1608  m_optionsObj->m_ov.m_rawChainDataOutputFileName,
1609  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1610  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1611  if ((m_env.subDisplayFile() ) &&
1612  (m_optionsObj->m_ov.m_totallyMute == false)) {
1613  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1614  << ": just wrote (per period request) " << m_numPositionsNotSubWritten << " chain positions "
1615  << ", " << 0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod << " <= pos <= " << 0
1616  << std::endl;
1617  }
1618 
1619  if (writeLogLikelihood) {
1620  workingLogLikelihoodValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1621  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1622  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_loglikelihood",
1623  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1624  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1625  }
1626 
1627  if (writeLogTarget) {
1628  workingLogTargetValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1629  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1630  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_logtarget",
1631  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1632  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1633  }
1634 
1635  m_numPositionsNotSubWritten = 0;
1636  }
1637 
1638  if (workingLogLikelihoodValues) (*workingLogLikelihoodValues)[0] = currentPositionData.logLikelihood();
1639  if (workingLogTargetValues ) (*workingLogTargetValues )[0] = currentPositionData.logTarget();
1640  if (true/*m_uniqueChainGenerate*/) m_idsOfUniquePositions[uniquePos++] = 0;
1641  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1642  m_logTargets [0] = currentPositionData.logTarget();
1643  m_alphaQuotients[0] = 1.;
1644  }
1645  //*m_env.subDisplayFile() << "AQUI 002" << std::endl;
1646 
1647  if ((m_env.subDisplayFile() ) &&
1648  (m_env.displayVerbosity() >= 10 ) &&
1649  (m_optionsObj->m_ov.m_totallyMute == false)) {
1650  *m_env.subDisplayFile() << "\n"
1651  << "\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
1652  << "\n"
1653  << std::endl;
1654  }
1655 
1656  //m_env.syncPrintDebugMsg("In MetropolisHastingsSG<P_V,P_M>::generateFullChain(), right before main loop",3,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
1657 
1658  //****************************************************
1659  // Begin chain loop from positionId = 1
1660  //****************************************************
1661  if ((m_env.numSubEnvironments() < (unsigned int) m_env.fullComm().NumProc()) &&
1662  (m_initialPosition.numOfProcsForStorage() == 1 ) &&
1663  (m_env.subRank() != 0 )) {
1664  // subRank != 0 --> Enter the barrier and wait for processor 0 to decide to call the targetPdf
1665  double aux = 0.;
1666  aux = m_targetPdfSynchronizer->callFunction(NULL,
1667  NULL,
1668  NULL,
1669  NULL,
1670  NULL,
1671  NULL,
1672  NULL);
1673  if (aux) {}; // just to remove compiler warning
1674  for (unsigned int positionId = 1; positionId < workingChain.subSequenceSize(); ++positionId) {
1675  // Multiply by position values by 'positionId' in order to avoid a constant sequence,
1676  // which would cause zero variance and eventually OVERFLOW flags raised
1677  workingChain.setPositionValues(positionId,((double) positionId) * currentPositionData.vecValues());
1678  m_rawChainInfo.numRejections++;
1679  }
1680  }
1681  else for (unsigned int positionId = 1; positionId < workingChain.subSequenceSize(); ++positionId) {
1682  //****************************************************
1683  // Point 1/6 of logic for new position
1684  // Loop: initialize variables and print some information
1685  //****************************************************
1686  m_positionIdForDebugging = positionId;
1687  if ((m_env.subDisplayFile() ) &&
1688  (m_env.displayVerbosity() >= 3 ) &&
1689  (m_optionsObj->m_ov.m_totallyMute == false)) {
1690  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1691  << ": beginning chain position of id = " << positionId
1692  << ", m_optionsObj->m_ov.m_drMaxNumExtraStages = " << m_optionsObj->m_ov.m_drMaxNumExtraStages
1693  << std::endl;
1694  }
1695  unsigned int stageId = 0;
1696  m_stageIdForDebugging = stageId;
1697 
1698  m_tk->clearPreComputingPositions();
1699 
1700  if ((m_env.subDisplayFile() ) &&
1701  (m_env.displayVerbosity() >= 5 ) &&
1702  (m_optionsObj->m_ov.m_totallyMute == false)) {
1703  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1704  << ": about to set TK pre computing position of local id " << 0
1705  << ", values = " << currentPositionData.vecValues()
1706  << std::endl;
1707  }
1708  bool validPreComputingPosition = m_tk->setPreComputingPosition(currentPositionData.vecValues(),0);
1709  if ((m_env.subDisplayFile() ) &&
1710  (m_env.displayVerbosity() >= 5 ) &&
1711  (m_optionsObj->m_ov.m_totallyMute == false)) {
1712  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1713  << ": returned from setting TK pre computing position of local id " << 0
1714  << ", values = " << currentPositionData.vecValues()
1715  << ", valid = " << validPreComputingPosition
1716  << std::endl;
1717  }
1718  UQ_FATAL_TEST_MACRO(validPreComputingPosition == false,
1719  m_env.worldRank(),
1720  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
1721  "initial position should not be an invalid pre computing position");
1722 
1723  //****************************************************
1724  // Point 2/6 of logic for new position
1725  // Loop: generate new position
1726  //****************************************************
1727  bool keepGeneratingCandidates = true;
1728  while (keepGeneratingCandidates) {
1729  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) {
1730  iRC = gettimeofday(&timevalCandidate, NULL);
1731  }
1732 
1733  m_tk->rv(0).realizer().realization(tmpVecValues);
1734 
1735  if (m_numDisabledParameters > 0) { // gpmsa2
1736  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
1737  if (m_parameterEnabledStatus[paramId] == false) {
1738  tmpVecValues[paramId] = m_initialPosition[paramId];
1739  }
1740  }
1741  }
1742  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.candidateRunTime += MiscGetEllapsedSeconds(&timevalCandidate);
1743 
1744  outOfTargetSupport = !m_targetPdf.domainSet().contains(tmpVecValues);
1745 
1746  bool displayDetail = (m_env.displayVerbosity() >= 10/*99*/) || m_optionsObj->m_ov.m_displayCandidates;
1747  if ((m_env.subDisplayFile() ) &&
1748  (displayDetail ) &&
1749  (m_optionsObj->m_ov.m_totallyMute == false)) {
1750  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1751  << ": for chain position of id = " << positionId
1752  << ", candidate = " << tmpVecValues // FIX ME: might need parallelism
1753  << ", outOfTargetSupport = " << outOfTargetSupport
1754  << std::endl;
1755  }
1756 
1757  if (m_optionsObj->m_ov.m_putOutOfBoundsInChain) keepGeneratingCandidates = false;
1758  else keepGeneratingCandidates = outOfTargetSupport;
1759  }
1760 
1761  if ((m_env.subDisplayFile() ) &&
1762  (m_env.displayVerbosity() >= 5 ) &&
1763  (m_optionsObj->m_ov.m_totallyMute == false)) {
1764  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1765  << ": about to set TK pre computing position of local id " << stageId+1
1766  << ", values = " << tmpVecValues
1767  << std::endl;
1768  }
1769  validPreComputingPosition = m_tk->setPreComputingPosition(tmpVecValues,stageId+1);
1770  if ((m_env.subDisplayFile() ) &&
1771  (m_env.displayVerbosity() >= 5 ) &&
1772  (m_optionsObj->m_ov.m_totallyMute == false)) {
1773  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1774  << ": returned from setting TK pre computing position of local id " << stageId+1
1775  << ", values = " << tmpVecValues
1776  << ", valid = " << validPreComputingPosition
1777  << std::endl;
1778  }
1779 
1780  if (outOfTargetSupport) {
1781  m_rawChainInfo.numOutOfTargetSupport++;
1782  logPrior = -INFINITY;
1783  logLikelihood = -INFINITY;
1784  logTarget = -INFINITY;
1785  }
1786  else {
1787  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalTarget, NULL);
1788  logTarget = m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1789  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime += MiscGetEllapsedSeconds(&timevalTarget);
1790  m_rawChainInfo.numTargetCalls++;
1791  if ((m_env.subDisplayFile() ) &&
1792  (m_env.displayVerbosity() >= 3 ) &&
1793  (m_optionsObj->m_ov.m_totallyMute == false)) {
1794  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1795  << ": just returned from likelihood() for chain position of id " << positionId
1796  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1797  << ", logPrior = " << logPrior
1798  << ", logLikelihood = " << logLikelihood
1799  << ", logTarget = " << logTarget
1800  << std::endl;
1801  }
1802  }
1803  currentCandidateData.set(tmpVecValues,
1804  outOfTargetSupport,
1805  logLikelihood,
1806  logTarget);
1807 
1808  if ((m_env.subDisplayFile() ) &&
1809  (m_env.displayVerbosity() >= 10 ) &&
1810  (m_optionsObj->m_ov.m_totallyMute == false)) {
1811  *m_env.subDisplayFile() << "\n"
1812  << "\n-----------------------------------------------------------\n"
1813  << "\n"
1814  << std::endl;
1815  }
1816  bool accept = false;
1817  double alphaFirstCandidate = 0.;
1818  if (outOfTargetSupport) {
1819  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1820  m_alphaQuotients[positionId] = 0.;
1821  }
1822  }
1823  else {
1824  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalMhAlpha, NULL);
1825  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1826  alphaFirstCandidate = this->alpha(currentPositionData,currentCandidateData,0,1,&m_alphaQuotients[positionId]);
1827  }
1828  else {
1829  alphaFirstCandidate = this->alpha(currentPositionData,currentCandidateData,0,1,NULL);
1830  }
1831  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.mhAlphaRunTime += MiscGetEllapsedSeconds(&timevalMhAlpha);
1832  if ((m_env.subDisplayFile() ) &&
1833  (m_env.displayVerbosity() >= 10 ) &&
1834  (m_optionsObj->m_ov.m_totallyMute == false)) {
1835  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1836  << ": for chain position of id = " << positionId
1837  << std::endl;
1838  }
1839  accept = acceptAlpha(alphaFirstCandidate);
1840  }
1841 
1842  bool displayDetail = (m_env.displayVerbosity() >= 10/*99*/) || m_optionsObj->m_ov.m_displayCandidates;
1843  if ((m_env.subDisplayFile() ) &&
1844  (displayDetail ) &&
1845  (m_optionsObj->m_ov.m_totallyMute == false)) {
1846  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1847  << ": for chain position of id = " << positionId
1848  << ", outOfTargetSupport = " << outOfTargetSupport
1849  << ", alpha = " << alphaFirstCandidate
1850  << ", accept = " << accept
1851  << ", currentCandidateData.vecValues() = ";
1852  *m_env.subDisplayFile() << currentCandidateData.vecValues(); // FIX ME: might need parallelism
1853  *m_env.subDisplayFile() << "\n"
1854  << "\n curLogTarget = " << currentPositionData.logTarget()
1855  << "\n"
1856  << "\n canLogTarget = " << currentCandidateData.logTarget()
1857  << "\n"
1858  << std::endl;
1859  }
1860  if ((m_env.subDisplayFile() ) &&
1861  (m_env.displayVerbosity() >= 10 ) &&
1862  (m_optionsObj->m_ov.m_totallyMute == false)) {
1863  *m_env.subDisplayFile() << "\n"
1864  << "\n-----------------------------------------------------------\n"
1865  << "\n"
1866  << std::endl;
1867  }
1868 
1869  //****************************************************
1870  // Point 3/6 of logic for new position
1871  // Loop: delayed rejection
1872  //****************************************************
1873  // sep2011
1874  std::vector<MarkovChainPositionData<P_V>*> drPositionsData(stageId+2,NULL);
1875  std::vector<unsigned int> tkStageIds (stageId+2,0);
1876  if ((accept == false) &&
1877  (outOfTargetSupport == false) && // IMPORTANT
1878  (m_optionsObj->m_ov.m_drMaxNumExtraStages > 0 )) {
1879  if ((m_optionsObj->m_ov.m_drDuringAmNonAdaptiveInt == false ) &&
1880  (m_optionsObj->m_ov.m_tkUseLocalHessian == false ) &&
1881  (m_optionsObj->m_ov.m_amInitialNonAdaptInterval > 0 ) &&
1882  (m_optionsObj->m_ov.m_amAdaptInterval > 0 ) &&
1883  (positionId <= m_optionsObj->m_ov.m_amInitialNonAdaptInterval)) {
1884  // Avoid DR now
1885  }
1886  else {
1887  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalDR, NULL);
1888 
1889  drPositionsData[0] = new MarkovChainPositionData<P_V>(currentPositionData );
1890  drPositionsData[1] = new MarkovChainPositionData<P_V>(currentCandidateData);
1891 
1892  tkStageIds[0] = 0;
1893  tkStageIds[1] = 1;
1894 
1895  while ((validPreComputingPosition == true ) &&
1896  (accept == false ) &&
1897  (stageId < m_optionsObj->m_ov.m_drMaxNumExtraStages)) {
1898  if ((m_env.subDisplayFile() ) &&
1899  (m_env.displayVerbosity() >= 10 ) &&
1900  (m_optionsObj->m_ov.m_totallyMute == false)) {
1901  *m_env.subDisplayFile() << "\n"
1902  << "\n+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"
1903  << "\n"
1904  << std::endl;
1905  }
1906  m_rawChainInfo.numDRs++;
1907  stageId++;
1908  m_stageIdForDebugging = stageId;
1909  if ((m_env.subDisplayFile() ) &&
1910  (m_env.displayVerbosity() >= 10 ) &&
1911  (m_optionsObj->m_ov.m_totallyMute == false)) {
1912  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1913  << ": for chain position of id = " << positionId
1914  << ", beginning stageId = " << stageId
1915  << std::endl;
1916  }
1917 
1918  keepGeneratingCandidates = true;
1919  while (keepGeneratingCandidates) {
1920  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalCandidate, NULL);
1921  m_tk->rv(tkStageIds).realizer().realization(tmpVecValues);
1922  if (m_numDisabledParameters > 0) { // gpmsa2
1923  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
1924  if (m_parameterEnabledStatus[paramId] == false) {
1925  tmpVecValues[paramId] = m_initialPosition[paramId];
1926  }
1927  }
1928  }
1929  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.candidateRunTime += MiscGetEllapsedSeconds(&timevalCandidate);
1930 
1931  outOfTargetSupport = !m_targetPdf.domainSet().contains(tmpVecValues);
1932 
1933  if (m_optionsObj->m_ov.m_putOutOfBoundsInChain) keepGeneratingCandidates = false;
1934  else keepGeneratingCandidates = outOfTargetSupport;
1935  }
1936 
1937  if ((m_env.subDisplayFile() ) &&
1938  (m_env.displayVerbosity() >= 5 ) &&
1939  (m_optionsObj->m_ov.m_totallyMute == false)) {
1940  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1941  << ": about to set TK pre computing position of local id " << stageId+1
1942  << ", values = " << tmpVecValues
1943  << std::endl;
1944  }
1945  validPreComputingPosition = m_tk->setPreComputingPosition(tmpVecValues,stageId+1);
1946  if ((m_env.subDisplayFile() ) &&
1947  (m_env.displayVerbosity() >= 5 ) &&
1948  (m_optionsObj->m_ov.m_totallyMute == false)) {
1949  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1950  << ": returned from setting TK pre computing position of local id " << stageId+1
1951  << ", values = " << tmpVecValues
1952  << ", valid = " << validPreComputingPosition
1953  << std::endl;
1954  }
1955 
1956  if (outOfTargetSupport) {
1957  m_rawChainInfo.numOutOfTargetSupportInDR++; // new 2010/May/12
1958  logPrior = -INFINITY;
1959  logLikelihood = -INFINITY;
1960  logTarget = -INFINITY;
1961  }
1962  else {
1963  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalTarget, NULL);
1964  logTarget = m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1965  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime += MiscGetEllapsedSeconds(&timevalTarget);
1966  m_rawChainInfo.numTargetCalls++;
1967  if ((m_env.subDisplayFile() ) &&
1968  (m_env.displayVerbosity() >= 3 ) &&
1969  (m_optionsObj->m_ov.m_totallyMute == false)) {
1970  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1971  << ": just returned from likelihood() for chain position of id " << positionId
1972  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1973  << ", stageId = " << stageId
1974  << ", logPrior = " << logPrior
1975  << ", logLikelihood = " << logLikelihood
1976  << ", logTarget = " << logTarget
1977  << std::endl;
1978  }
1979  }
1980  currentCandidateData.set(tmpVecValues,
1981  outOfTargetSupport,
1982  logLikelihood,
1983  logTarget);
1984 
1985  drPositionsData.push_back(new MarkovChainPositionData<P_V>(currentCandidateData));
1986  tkStageIds.push_back (stageId+1);
1987 
1988  double alphaDR = 0.;
1989  if (outOfTargetSupport == false) {
1990  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalDrAlpha, NULL);
1991  alphaDR = this->alpha(drPositionsData,tkStageIds);
1992  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.drAlphaRunTime += MiscGetEllapsedSeconds(&timevalDrAlpha);
1993  accept = acceptAlpha(alphaDR);
1994  }
1995 
1996  displayDetail = (m_env.displayVerbosity() >= 10/*99*/) || m_optionsObj->m_ov.m_displayCandidates;
1997  if ((m_env.subDisplayFile() ) &&
1998  (displayDetail ) &&
1999  (m_optionsObj->m_ov.m_totallyMute == false)) {
2000  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2001  << ": for chain position of id = " << positionId
2002  << " and stageId = " << stageId
2003  << ", outOfTargetSupport = " << outOfTargetSupport
2004  << ", alpha = " << alphaDR
2005  << ", accept = " << accept
2006  << ", currentCandidateData.vecValues() = ";
2007  *m_env.subDisplayFile() << currentCandidateData.vecValues(); // FIX ME: might need parallelism
2008  *m_env.subDisplayFile() << std::endl;
2009  }
2010  } // while
2011 
2012  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.drRunTime += MiscGetEllapsedSeconds(&timevalDR);
2013  } // if-else "Avoid DR now"
2014  } // end of 'delayed rejection' logic
2015 
2016  for (unsigned int i = 0; i < drPositionsData.size(); ++i) {
2017  if (drPositionsData[i]) delete drPositionsData[i];
2018  }
2019 
2020  //****************************************************
2021  // Point 4/6 of logic for new position
2022  // Loop: update chain
2023  //****************************************************
2024  if (accept) {
2025  workingChain.setPositionValues(positionId,currentCandidateData.vecValues());
2026  if (true/*m_uniqueChainGenerate*/) m_idsOfUniquePositions[uniquePos++] = positionId;
2027  currentPositionData = currentCandidateData;
2028  }
2029  else {
2030  workingChain.setPositionValues(positionId,currentPositionData.vecValues());
2031  m_rawChainInfo.numRejections++;
2032  }
2033  m_numPositionsNotSubWritten++;
2034  if ((m_optionsObj->m_ov.m_rawChainDataOutputPeriod > 0 ) &&
2035  (((positionId+1) % m_optionsObj->m_ov.m_rawChainDataOutputPeriod) == 0 ) &&
2036  (m_optionsObj->m_ov.m_rawChainDataOutputFileName != ".")) {
2037  if ((m_env.subDisplayFile() ) &&
2038  (m_env.displayVerbosity() >= 10 ) &&
2039  (m_optionsObj->m_ov.m_totallyMute == false)) {
2040  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2041  << ", for chain position of id = " << positionId
2042  << ": about to write (per period request) " << m_numPositionsNotSubWritten << " chain positions "
2043  << ", " << positionId + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod << " <= pos <= " << positionId
2044  << std::endl;
2045  }
2046  workingChain.subWriteContents(positionId + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
2047  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
2048  m_optionsObj->m_ov.m_rawChainDataOutputFileName,
2049  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
2050  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
2051  if ((m_env.subDisplayFile() ) &&
2052  (m_optionsObj->m_ov.m_totallyMute == false)) {
2053  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2054  << ", for chain position of id = " << positionId
2055  << ": just wrote (per period request) " << m_numPositionsNotSubWritten << " chain positions "
2056  << ", " << positionId + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod << " <= pos <= " << positionId
2057  << std::endl;
2058  }
2059 
2060  if (writeLogLikelihood) {
2061  workingLogLikelihoodValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
2062  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
2063  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_loglikelihood",
2064  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
2065  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
2066  }
2067 
2068  if (writeLogTarget) {
2069  workingLogTargetValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
2070  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
2071  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_logtarget",
2072  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
2073  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
2074  }
2075 
2076  m_numPositionsNotSubWritten = 0;
2077  }
2078 
2079 
2080  if (workingLogLikelihoodValues) (*workingLogLikelihoodValues)[positionId] = currentPositionData.logLikelihood();
2081  if (workingLogTargetValues ) (*workingLogTargetValues )[positionId] = currentPositionData.logTarget();
2082 
2083  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
2084  m_logTargets[positionId] = currentPositionData.logTarget();
2085  }
2086 
2087  if( m_optionsObj->m_ov.m_enableBrooksGelmanConvMonitor > 0 ) {
2088  if( positionId%m_optionsObj->m_ov.m_enableBrooksGelmanConvMonitor == 0 &&
2089  positionId > m_optionsObj->m_ov.m_BrooksGelmanLag+1 ) { //+1 to help ensure there are at least 2 samples to use
2090 
2091  double conv_est = workingChain.estimateConvBrooksGelman( m_optionsObj->m_ov.m_BrooksGelmanLag,
2092  positionId - m_optionsObj->m_ov.m_BrooksGelmanLag );
2093 
2094  if ( m_env.subDisplayFile() ) {
2095  *m_env.subDisplayFile() << "positionId = " << positionId
2096  << ", conv_est = " << conv_est << std::endl;
2097  (*m_env.subDisplayFile()).flush();
2098  }
2099  }
2100  }
2101 
2102  //****************************************************
2103  // Point 5/6 of logic for new position
2104  // Loop: adaptive Metropolis (adaptation of covariance matrix)
2105  //****************************************************
2106  // sep2011
2107  if ((m_optionsObj->m_ov.m_tkUseLocalHessian == false) && // IMPORTANT
2108  (m_optionsObj->m_ov.m_amInitialNonAdaptInterval > 0) &&
2109  (m_optionsObj->m_ov.m_amAdaptInterval > 0)) {
2110  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalAM, NULL);
2111 
2112  // Now might be the moment to adapt
2113  unsigned int idOfFirstPositionInSubChain = 0;
2114  SequenceOfVectors<P_V,P_M> partialChain(m_vectorSpace,0,m_optionsObj->m_prefix+"partialChain");
2115 
2116  // Check if now is indeed the moment to adapt
2117  bool printAdaptedMatrix = false;
2118  if (positionId < m_optionsObj->m_ov.m_amInitialNonAdaptInterval) {
2119  // Do nothing
2120  }
2121  else if (positionId == m_optionsObj->m_ov.m_amInitialNonAdaptInterval) {
2122  idOfFirstPositionInSubChain = 0;
2123  partialChain.resizeSequence(m_optionsObj->m_ov.m_amInitialNonAdaptInterval+1);
2124  m_lastMean = m_vectorSpace.newVector();
2125  m_lastAdaptedCovMatrix = m_vectorSpace.newMatrix();
2126  printAdaptedMatrix = true;
2127  }
2128  else {
2129  unsigned int interval = positionId - m_optionsObj->m_ov.m_amInitialNonAdaptInterval;
2130  if ((interval % m_optionsObj->m_ov.m_amAdaptInterval) == 0) {
2131  idOfFirstPositionInSubChain = positionId - m_optionsObj->m_ov.m_amAdaptInterval;
2132  partialChain.resizeSequence(m_optionsObj->m_ov.m_amAdaptInterval);
2133 
2134  if (m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputPeriod > 0) {
2135  if ((interval % m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputPeriod) == 0) {
2136  printAdaptedMatrix = true;
2137  }
2138  }
2139  }
2140  }
2141 
2142  // If now is indeed the moment to adapt, then do it!
2143  if (partialChain.subSequenceSize() > 0) {
2144  P_V transporterVec(m_vectorSpace.zeroVector());
2145  for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2146  workingChain.getPositionValues(idOfFirstPositionInSubChain+i,transporterVec);
2147 
2148  // Transform to the space without boundaries. This is the space
2149  // where the proposal distribution is Gaussian
2150  if (this->m_optionsObj->m_ov.m_doLogitTransform == true) {
2151  // Only do this when we don't use the Hessian (this may change in
2152  // future, but transformToGaussianSpace() is only implemented in
2153  // TransformedScaledCovMatrixTKGroup
2154  P_V transformedTransporterVec(m_vectorSpace.zeroVector());
2155  dynamic_cast<TransformedScaledCovMatrixTKGroup<P_V, P_M>* >(
2156  m_tk)->transformToGaussianSpace(transporterVec,
2157  transformedTransporterVec);
2158  partialChain.setPositionValues(i, transformedTransporterVec);
2159  }
2160  else {
2161  partialChain.setPositionValues(i, transporterVec);
2162  }
2163  }
2164  updateAdaptedCovMatrix(partialChain,
2165  idOfFirstPositionInSubChain,
2166  m_lastChainSize,
2167  *m_lastMean,
2168  *m_lastAdaptedCovMatrix);
2169 
2170  if ((printAdaptedMatrix == true) &&
2171  (m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputFileName != "." )) { // palms
2172  char varNamePrefix[64];
2173  sprintf(varNamePrefix,"mat_am%d",positionId);
2174 
2175  char tmpChar[64];
2176  sprintf(tmpChar,"_am%d",positionId);
2177 
2178  std::set<unsigned int> tmpSet;
2179  tmpSet.insert(m_env.subId());
2180 
2181  m_lastAdaptedCovMatrix->subWriteContents(varNamePrefix,
2182  (m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputFileName+tmpChar),
2183  m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputFileType,
2184  tmpSet);
2185  if ((m_env.subDisplayFile() ) &&
2186  (m_optionsObj->m_ov.m_totallyMute == false)) {
2187  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2188  << ": just wrote last adapted proposal cov matrix contents = " << *m_lastAdaptedCovMatrix
2189  << std::endl;
2190  }
2191  } // if (printAdaptedMatrix && ...)
2192 
2193  bool tmpCholIsPositiveDefinite = false;
2194  P_M tmpChol(*m_lastAdaptedCovMatrix);
2195  P_M attemptedMatrix(tmpChol);
2196  if ((m_env.subDisplayFile() ) &&
2197  (m_env.displayVerbosity() >= 10)) {
2198  //(m_optionsObj->m_ov.m_totallyMute == false)) {
2199  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2200  << ", positionId = " << positionId
2201  << ": 'am' calling first tmpChol.chol()"
2202  << std::endl;
2203  }
2204  iRC = tmpChol.chol();
2205  if (iRC) {
2206  std::cerr << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain(): first chol failed\n";
2207  }
2208  if ((m_env.subDisplayFile() ) &&
2209  (m_env.displayVerbosity() >= 10)) {
2210  //(m_optionsObj->m_ov.m_totallyMute == false)) {
2211  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2212  << ", positionId = " << positionId
2213  << ": 'am' got first tmpChol.chol() with iRC = " << iRC
2214  << std::endl;
2215  if (iRC == 0) {
2216  double diagMult = 1.;
2217  for (unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2218  diagMult *= tmpChol(j,j);
2219  }
2220  *m_env.subDisplayFile() << "diagMult = " << diagMult
2221  << std::endl;
2222  }
2223  }
2224 #if 0 // tentative logic
2225  if (iRC == 0) {
2226  double diagMult = 1.;
2227  for (unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2228  diagMult *= tmpChol(j,j);
2229  }
2230  if (diagMult < 1.e-40) {
2231  iRC = UQ_MATRIX_IS_NOT_POS_DEFINITE_RC;
2232  }
2233  }
2234 #endif
2235 
2236  if (iRC) {
2237  UQ_FATAL_TEST_MACRO(iRC != UQ_MATRIX_IS_NOT_POS_DEFINITE_RC,
2238  m_env.worldRank(),
2239  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
2240  "invalid iRC returned from first chol()");
2241  // Matrix is not positive definite
2242  P_M* tmpDiag = m_vectorSpace.newDiagMatrix(m_optionsObj->m_ov.m_amEpsilon);
2243  if (m_numDisabledParameters > 0) { // gpmsa2
2244  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2245  if (m_parameterEnabledStatus[paramId] == false) {
2246  (*tmpDiag)(paramId,paramId) = 0.;
2247  }
2248  }
2249  }
2250  tmpChol = *m_lastAdaptedCovMatrix + *tmpDiag;
2251  attemptedMatrix = tmpChol;
2252  delete tmpDiag;
2253  if ((m_env.subDisplayFile() ) &&
2254  (m_env.displayVerbosity() >= 10)) {
2255  //(m_optionsObj->m_ov.m_totallyMute == false)) {
2256  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2257  << ", positionId = " << positionId
2258  << ": 'am' calling second tmpChol.chol()"
2259  << std::endl;
2260  }
2261  iRC = tmpChol.chol();
2262  if (iRC) {
2263  std::cerr << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain(): second chol failed\n";
2264  }
2265  if ((m_env.subDisplayFile() ) &&
2266  (m_env.displayVerbosity() >= 10)) {
2267  //(m_optionsObj->m_ov.m_totallyMute == false)) {
2268  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2269  << ", positionId = " << positionId
2270  << ": 'am' got second tmpChol.chol() with iRC = " << iRC
2271  << std::endl;
2272  if (iRC == 0) {
2273  double diagMult = 1.;
2274  for (unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2275  diagMult *= tmpChol(j,j);
2276  }
2277  *m_env.subDisplayFile() << "diagMult = " << diagMult
2278  << std::endl;
2279  }
2280  else {
2281  *m_env.subDisplayFile() << "attemptedMatrix = " << attemptedMatrix // FIX ME: might demand parallelism
2282  << std::endl;
2283  }
2284  }
2285  if (iRC) {
2286  UQ_FATAL_TEST_MACRO(iRC != UQ_MATRIX_IS_NOT_POS_DEFINITE_RC,
2287  m_env.worldRank(),
2288  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
2289  "invalid iRC returned from second chol()");
2290  // Do nothing
2291  }
2292  else {
2293  tmpCholIsPositiveDefinite = true;
2294  }
2295  }
2296  else {
2297  tmpCholIsPositiveDefinite = true;
2298  }
2299  if (tmpCholIsPositiveDefinite) {
2300  P_M tmpMatrix(m_optionsObj->m_ov.m_amEta*attemptedMatrix);
2301  if (m_numDisabledParameters > 0) { // gpmsa2
2302  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2303  if (m_parameterEnabledStatus[paramId] == false) {
2304  for (unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2305  tmpMatrix(i,paramId) = 0.;
2306  }
2307  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2308  tmpMatrix(paramId,j) = 0.;
2309  }
2310  tmpMatrix(paramId,paramId) = 1.;
2311  }
2312  }
2313  }
2314  if (this->m_optionsObj->m_ov.m_doLogitTransform) {
2315  (dynamic_cast<TransformedScaledCovMatrixTKGroup<P_V,P_M>* >(m_tk))
2316  ->updateLawCovMatrix(tmpMatrix);
2317  }
2318  else{
2319  (dynamic_cast<ScaledCovMatrixTKGroup<P_V,P_M>* >(m_tk))
2320  ->updateLawCovMatrix(tmpMatrix);
2321  }
2322 
2323 #ifdef UQ_DRAM_MCG_REQUIRES_INVERTED_COV_MATRICES
2324  UQ_FATAL_RC_MACRO(UQ_INCOMPLETE_IMPLEMENTATION_RC,
2325  m_env.worldRank(),
2326  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
2327  "need to code the update of m_upperCholProposalPrecMatrices");
2328 #endif
2329  }
2330 
2331  //for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2332  // if (partialChain[i]) delete partialChain[i];
2333  //}
2334  } // if (partialChain.subSequenceSize() > 0)
2335 
2336  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.amRunTime += MiscGetEllapsedSeconds(&timevalAM);
2337  } // End of 'adaptive Metropolis' logic
2338 
2339  //****************************************************
2340  // Point 6/6 of logic for new position
2341  // Loop: print some information before going to the next chain position
2342  //****************************************************
2343  if ((m_env.subDisplayFile() ) &&
2344  (m_env.displayVerbosity() >= 3 ) &&
2345  (m_optionsObj->m_ov.m_totallyMute == false)) {
2346  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2347  << ": finishing chain position of id = " << positionId
2348  << ", accept = " << accept
2349  << ", curLogTarget = " << currentPositionData.logTarget()
2350  << ", canLogTarget = " << currentCandidateData.logTarget()
2351  << std::endl;
2352  }
2353 
2354  if ((m_optionsObj->m_ov.m_rawChainDisplayPeriod > 0) &&
2355  (((positionId+1) % m_optionsObj->m_ov.m_rawChainDisplayPeriod) == 0)) {
2356  if ((m_env.subDisplayFile() ) &&
2357  (m_optionsObj->m_ov.m_totallyMute == false)) {
2358  *m_env.subDisplayFile() << "Finished generating " << positionId+1
2359  << " positions" // root
2360  << ", curret rejection percentage = " << (100. * ((double) m_rawChainInfo.numRejections)/((double) (positionId+1)))
2361  << " %"
2362  << std::endl;
2363  }
2364  }
2365 
2366  if ((m_env.subDisplayFile() ) &&
2367  (m_env.displayVerbosity() >= 10 ) &&
2368  (m_optionsObj->m_ov.m_totallyMute == false)) {
2369  *m_env.subDisplayFile() << "\n"
2370  << "\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
2371  << "\n"
2372  << std::endl;
2373  }
2374  } // end chain loop [for (unsigned int positionId = 1; positionId < workingChain.subSequenceSize(); ++positionId) {]
2375 
2376  if ((m_env.numSubEnvironments() < (unsigned int) m_env.fullComm().NumProc()) &&
2377  (m_initialPosition.numOfProcsForStorage() == 1 ) &&
2378  (m_env.subRank() == 0 )) {
2379  // subRank == 0 --> Tell all other processors to exit barrier now that the chain has been fully generated
2380  double aux = 0.;
2381  aux = m_targetPdfSynchronizer->callFunction(NULL,
2382  NULL,
2383  NULL,
2384  NULL,
2385  NULL,
2386  NULL,
2387  NULL);
2388  if (aux) {}; // just to remove compiler warning
2389  }
2390 
2391  //****************************************************
2392  // Print basic information about the chain
2393  //****************************************************
2394  m_rawChainInfo.runTime += MiscGetEllapsedSeconds(&timevalChain);
2395  if ((m_env.subDisplayFile() ) &&
2396  (m_optionsObj->m_ov.m_totallyMute == false)) {
2397  *m_env.subDisplayFile() << "Finished the generation of Markov chain " << workingChain.name()
2398  << ", with sub " << workingChain.subSequenceSize()
2399  << " positions";
2400  *m_env.subDisplayFile() << "\nSome information about this chain:"
2401  << "\n Chain run time = " << m_rawChainInfo.runTime
2402  << " seconds";
2403  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) {
2404  *m_env.subDisplayFile() << "\n\n Breaking of the chain run time:\n";
2405  *m_env.subDisplayFile() << "\n Candidate run time = " << m_rawChainInfo.candidateRunTime
2406  << " seconds (" << 100.*m_rawChainInfo.candidateRunTime/m_rawChainInfo.runTime
2407  << "%)";
2408  *m_env.subDisplayFile() << "\n Num target calls = " << m_rawChainInfo.numTargetCalls;
2409  *m_env.subDisplayFile() << "\n Target d. run time = " << m_rawChainInfo.targetRunTime
2410  << " seconds (" << 100.*m_rawChainInfo.targetRunTime/m_rawChainInfo.runTime
2411  << "%)";
2412  *m_env.subDisplayFile() << "\n Avg target run time = " << m_rawChainInfo.targetRunTime/((double) m_rawChainInfo.numTargetCalls)
2413  << " seconds";
2414  *m_env.subDisplayFile() << "\n Mh alpha run time = " << m_rawChainInfo.mhAlphaRunTime
2415  << " seconds (" << 100.*m_rawChainInfo.mhAlphaRunTime/m_rawChainInfo.runTime
2416  << "%)";
2417  *m_env.subDisplayFile() << "\n Dr alpha run time = " << m_rawChainInfo.drAlphaRunTime
2418  << " seconds (" << 100.*m_rawChainInfo.drAlphaRunTime/m_rawChainInfo.runTime
2419  << "%)";
2420  *m_env.subDisplayFile() << "\n---------------------- --------------";
2421  double sumRunTime = m_rawChainInfo.candidateRunTime + m_rawChainInfo.targetRunTime + m_rawChainInfo.mhAlphaRunTime + m_rawChainInfo.drAlphaRunTime;
2422  *m_env.subDisplayFile() << "\n Sum = " << sumRunTime
2423  << " seconds (" << 100.*sumRunTime/m_rawChainInfo.runTime
2424  << "%)";
2425  *m_env.subDisplayFile() << "\n\n Other run times:";
2426  *m_env.subDisplayFile() << "\n DR run time = " << m_rawChainInfo.drRunTime
2427  << " seconds (" << 100.*m_rawChainInfo.drRunTime/m_rawChainInfo.runTime
2428  << "%)";
2429  *m_env.subDisplayFile() << "\n AM run time = " << m_rawChainInfo.amRunTime
2430  << " seconds (" << 100.*m_rawChainInfo.amRunTime/m_rawChainInfo.runTime
2431  << "%)";
2432  }
2433  *m_env.subDisplayFile() << "\n Number of DRs = " << m_rawChainInfo.numDRs << "(num_DRs/chain_size = " << (double) m_rawChainInfo.numDRs/(double) workingChain.subSequenceSize()
2434  << ")";
2435  *m_env.subDisplayFile() << "\n Out of target support in DR = " << m_rawChainInfo.numOutOfTargetSupportInDR;
2436  *m_env.subDisplayFile() << "\n Rejection percentage = " << 100. * (double) m_rawChainInfo.numRejections/(double) workingChain.subSequenceSize()
2437  << " %";
2438  *m_env.subDisplayFile() << "\n Out of target support percentage = " << 100. * (double) m_rawChainInfo.numOutOfTargetSupport/(double) workingChain.subSequenceSize()
2439  << " %";
2440  *m_env.subDisplayFile() << std::endl;
2441  }
2442 
2443  //****************************************************
2444  // Release memory before leaving routine
2445  //****************************************************
2446  //m_env.syncPrintDebugMsg("Leaving MetropolisHastingsSG<P_V,P_M>::generateFullChain()",3,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
2447 
2448  return;
2449 }
2450 //--------------------------------------------------
2451 template <class P_V,class P_M>
2452 void
2453 MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix(
2454  const BaseVectorSequence<P_V,P_M>& partialChain,
2455  unsigned int idOfFirstPositionInSubChain,
2456  double& lastChainSize,
2457  P_V& lastMean,
2458  P_M& lastAdaptedCovMatrix)
2459 {
2460  double doubleSubChainSize = (double) partialChain.subSequenceSize();
2461  if (lastChainSize == 0) {
2462  UQ_FATAL_TEST_MACRO(partialChain.subSequenceSize() < 2,
2463  m_env.worldRank(),
2464  "MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix()",
2465  "'partialChain.subSequenceSize()' should be >= 2");
2466 
2467 #if 1 // prudenci-2012-07-06
2468  lastMean = partialChain.subMeanPlain();
2469 #else
2470  partialChain.subMeanExtra(0,partialChain.subSequenceSize(),lastMean);
2471 #endif
2472 
2473  P_V tmpVec(m_vectorSpace.zeroVector());
2474  lastAdaptedCovMatrix = -doubleSubChainSize * matrixProduct(lastMean,lastMean);
2475  for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2476  partialChain.getPositionValues(i,tmpVec);
2477  lastAdaptedCovMatrix += matrixProduct(tmpVec,tmpVec);
2478  }
2479  lastAdaptedCovMatrix /= (doubleSubChainSize - 1.); // That is why partialChain size must be >= 2
2480  }
2481  else {
2482  UQ_FATAL_TEST_MACRO(partialChain.subSequenceSize() < 1,
2483  m_env.worldRank(),
2484  "MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix()",
2485  "'partialChain.subSequenceSize()' should be >= 1");
2486 
2487  UQ_FATAL_TEST_MACRO(idOfFirstPositionInSubChain < 1,
2488  m_env.worldRank(),
2489  "MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix()",
2490  "'idOfFirstPositionInSubChain' should be >= 1");
2491 
2492  P_V tmpVec (m_vectorSpace.zeroVector());
2493  P_V diffVec(m_vectorSpace.zeroVector());
2494  for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2495  double doubleCurrentId = (double) (idOfFirstPositionInSubChain+i);
2496  partialChain.getPositionValues(i,tmpVec);
2497  diffVec = tmpVec - lastMean;
2498 
2499  double ratio1 = (1. - 1./doubleCurrentId); // That is why idOfFirstPositionInSubChain must be >= 1
2500  double ratio2 = (1./(1.+doubleCurrentId));
2501  lastAdaptedCovMatrix = ratio1 * lastAdaptedCovMatrix + ratio2 * matrixProduct(diffVec,diffVec);
2502  lastMean += ratio2 * diffVec;
2503  }
2504  }
2505  lastChainSize += doubleSubChainSize;
2506 
2507  if (m_numDisabledParameters > 0) { // gpmsa2
2508  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2509  if (m_parameterEnabledStatus[paramId] == false) {
2510  for (unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2511  lastAdaptedCovMatrix(i,paramId) = 0.;
2512  }
2513  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2514  lastAdaptedCovMatrix(paramId,j) = 0.;
2515  }
2516  lastAdaptedCovMatrix(paramId,paramId) = 1.;
2517  }
2518  }
2519  }
2520 
2521  return;
2522 }
2523 
2524 template <class P_V, class P_M>
2525 void
2526 MetropolisHastingsSG<P_V, P_M>::transformInitialCovMatrixToGaussianSpace(
2527  const BoxSubset<P_V, P_M> & boxSubset)
2528 {
2529  P_V min_domain_bounds(boxSubset.minValues());
2530  P_V max_domain_bounds(boxSubset.maxValues());
2531 
2532  for (unsigned int i = 0; i < min_domain_bounds.sizeLocal(); i++) {
2533  double min_val = min_domain_bounds[i];
2534  double max_val = max_domain_bounds[i];
2535 
2536  if (boost::math::isfinite(min_val) && boost::math::isfinite(max_val)) {
2537  if (m_initialProposalCovMatrix(i, i) >= max_val - min_val) {
2538  // User is trying to specify a uniform proposal distribution, which
2539  // is unsupported. Throw an error for now.
2540  std::cerr << "Proposal variance element "
2541  << i
2542  << " is "
2543  << m_initialProposalCovMatrix(i, i)
2544  << " but domain is of size "
2545  << max_val - min_val
2546  << std::endl;
2547  std::cerr << "QUESO does not support uniform-like proposal "
2548  << "distributions. Try making the proposal variance smaller"
2549  << std::endl;
2550  }
2551 
2552  // The jacobian at the midpoint of the domain
2553  double transformJacobian = 4.0 / (max_val - min_val);
2554 
2555  // Just do the multiplication by hand for now. There's no method in
2556  // Gsl(Vector|Matrix) to do this for me.
2557  for (unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
2558  // Multiply column j by element j
2559  m_initialProposalCovMatrix(j, i) *= transformJacobian;
2560  }
2561  for (unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
2562  // Multiply row j by element j
2563  m_initialProposalCovMatrix(i, j) *= transformJacobian;
2564  }
2565  }
2566  }
2567 }
2568 
2569 } // End namespace QUESO
2570 
2571 template class QUESO::MetropolisHastingsSG<QUESO::GslVector, QUESO::GslMatrix>;
QUESO::MarkovChainPositionData::outOfTargetSupport
bool outOfTargetSupport() const
Whether or not a position is out of target support; access to private attribute m_outOfTargetSupport...
Definition: MarkovChainPositionData.C:103
QUESO::BaseVectorSequence::setName
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
Definition: VectorSequence.C:128
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:287
QUESO::MhOptionsValues
This class provides options for the Metropolis-Hastings generator of samples if no input file is avai...
Definition: MetropolisHastingsSGOptions.h:110
QUESO::MetropolisHastingsSG
A templated class that represents a Metropolis-Hastings generator of samples.
Definition: MetropolisHastingsSG.h:120
QUESO::BaseVectorSequence::estimateConvBrooksGelman
virtual double estimateConvBrooksGelman(unsigned int initialPos, unsigned int numPos) const =0
Estimates convergence rate using Brooks &amp; Gelman method. See template specialization.
QUESO::BaseVectorSequence::setPositionValues
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
QUESO::MetropolisHastingsSG::transformInitialCovMatrixToGaussianSpace
void transformInitialCovMatrixToGaussianSpace(const BoxSubset< P_V, P_M > &boxSubset)
Definition: MetropolisHastingsSG.C:2526
QUESO::MHRawChainInfoStruct::reset
void reset()
Resets Metropolis-Hastings chain info.
Definition: MetropolisHastingsSG.C:83
QUESO::BoxSubset::maxValues
const V & maxValues() const
Vector of the maximum values of the box subset.
Definition: BoxSubset.C:88
QUESO::BaseVectorSequence::unifiedReadContents
virtual void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)=0
Reads info of the unified sequence from a file. See template specialization.
QUESO::MetropolisHastingsSG::generateSequence
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
Definition: MetropolisHastingsSG.C:975
QUESO::MLSamplingLevelOptions
This class provides options for each level of the Multilevel sequence generator if no input file is a...
Definition: MLSamplingLevelOptions.h:115
QUESO::MetropolisHastingsSG::transitionKernel
const BaseTKGroup< P_V, P_M > & transitionKernel() const
Returns the underlying transition kernel for this sequence generator.
Definition: MetropolisHastingsSG.C:961
QUESO::MetropolisHastingsSG::generateFullChain
void generateFullChain(const P_V &valuesOf1stPosition, unsigned int chainSize, BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
This method generates the chain.
Definition: MetropolisHastingsSG.C:1457
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:298
QUESO::BaseVectorSequence::unifiedWriteContents
virtual void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const =0
Writes info of the unified sequence to a file. See template specialization.
QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Definition: Defines.h:84
QUESO::BaseVectorSequence::subWriteContents
virtual void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const =0
Writes info of the sub-sequence to a file. See template specialization.
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:66
QUESO::MetropolisHastingsSG::commonConstructor
void commonConstructor()
Reads the options values from the options input file.
Definition: MetropolisHastingsSG.C:457
QUESO::MetropolisHastingsSG::updateAdaptedCovMatrix
void updateAdaptedCovMatrix(const BaseVectorSequence< P_V, P_M > &subChain, unsigned int idOfFirstPositionInSubChain, double &lastChainSize, P_V &lastMean, P_M &lastAdaptedCovMatrix)
This method updates the adapted covariance matrix.
Definition: MetropolisHastingsSG.C:2453
QUESO::BaseVectorSequence::getPositionValues
virtual void getPositionValues(unsigned int posId, V &vec) const =0
Gets the values of the sequence at position posId and stores them at vec. See template specialization...
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2411
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
RawValue_MPI_DOUBLE
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:48
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:274
RawValue_MPI_UNSIGNED
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:49
QUESO::MetropolisHastingsSGOptions
This class reads the options for the Metropolis-Hastings generator of samples from an input file...
Definition: MetropolisHastingsSGOptions.h:270
UQ_FATAL_RC_MACRO
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
Definition: Defines.h:207
QUESO::SequenceOfVectors::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
Definition: SequenceOfVectors.C:103
QUESO::InvLogitGaussianJointPdf::lawVarVector
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
Definition: InvLogitGaussianJointPdf.C:84
QUESO::HessianCovMatricesTKGroup
This class allows the representation of a transition kernel with Hessians.
Definition: HessianCovMatricesTKGroup.h:41
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:76
QUESO::BaseVectorSequence::vectorSizeLocal
unsigned int vectorSizeLocal() const
Local dimension (size) of the vector space.
Definition: VectorSequence.C:100
QUESO::ScalarSequence::subWriteContents
void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Writes the sub-sequence to a file.
Definition: ScalarSequence.C:2742
QUESO::BoxSubset::minValues
const V & minValues() const
Vector of the minimum values of the box subset.
Definition: BoxSubset.C:82
QUESO::InvLogitGaussianJointPdf::lawCovMatrix
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
Definition: InvLogitGaussianJointPdf.C:245
QUESO::MetropolisHastingsSG::writeInfo
int writeInfo(const BaseVectorSequence< P_V, P_M > &workingChain, std::ofstream &ofsvar) const
Writes information about the Markov chain in a file.
Definition: MetropolisHastingsSG.C:883
QUESO::MarkovChainPositionData::logTarget
double logTarget() const
Logarithm of the value of the target; access to private attribute m_logTarget.
Definition: MarkovChainPositionData.C:119
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:75
QUESO::MetropolisHastingsSG::alpha
double alpha(const MarkovChainPositionData< P_V > &x, const MarkovChainPositionData< P_V > &y, unsigned int xStageId, unsigned int yStageId, double *alphaQuotientPtr=NULL)
Calculates acceptance ration.
Definition: MetropolisHastingsSG.C:566
UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:89
QUESO::UQ_INCOMPLETE_IMPLEMENTATION_RC
const int UQ_INCOMPLETE_IMPLEMENTATION_RC
Definition: Defines.h:77
QUESO::MetropolisHastingsSG::MetropolisHastingsSG
MetropolisHastingsSG(const char *prefix, const MhOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, const P_M *inputProposalCovMatrix)
Constructor.
Definition: MetropolisHastingsSG.C:136
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:181
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:223
QUESO::MarkovChainPositionData::vecValues
const V & vecValues() const
Values of the chain (vector); access to private attribute m_vecValues.
Definition: MarkovChainPositionData.C:91
QUESO::TransformedScaledCovMatrixTKGroup
This class represents a transition kernel with a scaled covariance matrix on hybrid bounded/unbounded...
Definition: TransformedScaledCovMatrixTKGroup.h:45
QUESO::MpiComm::Allreduce
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
Definition: MpiComm.C:131
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2575
QUESO::BaseVectorSequence::subMeanPlain
const V & subMeanPlain() const
Finds the mean value of the sub-sequence.
Definition: VectorSequence.C:201
QUESO::MetropolisHastingsSG::getRawChainInfo
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
Definition: MetropolisHastingsSG.C:1437
QUESO::InvLogitGaussianJointPdf
A class for handling hybrid (transformed) Gaussians with bounds.
Definition: InvLogitGaussianJointPdf.h:55
queso_error
#define queso_error()
Definition: asserts.h:87
QUESO::SequenceOfVectors::getPositionValues
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
Definition: SequenceOfVectors.C:197
QUESO::MpiComm
The QUESO MPI Communicator Class.
Definition: MpiComm.h:75
QUESO::BaseVectorSequence::subMeanExtra
virtual void subMeanExtra(unsigned int initialPos, unsigned int numPos, V &meanVec) const =0
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
QUESO::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2799
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:299
QUESO::InvLogitGaussianJointPdf::lawExpVector
const V & lawExpVector() const
Access to the vector of mean values of the Gaussian (not transformed) and private attribute: m_lawExp...
Definition: InvLogitGaussianJointPdf.C:77
QUESO::SequenceOfVectors::resizeSequence
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
Definition: SequenceOfVectors.C:110
QUESO::BaseVectorSequence::subPositionsOfMaximum
double subPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
Definition: VectorSequence.C:386
QUESO::MHRawChainInfoStruct
A struct that represents a Metropolis-Hastings sample.
Definition: MetropolisHastingsSG.h:55
QUESO::ScaledCovMatrixTKGroup
This class allows the representation of a transition kernel with a scaled covariance matrix...
Definition: ScaledCovMatrixTKGroup.h:41
QUESO::MHRawChainInfoStruct::operator=
MHRawChainInfoStruct & operator=(const MHRawChainInfoStruct &rhs)
Assignment operator.
Definition: MetropolisHastingsSG.C:56
QUESO::BaseVectorSequence::name
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
Definition: VectorSequence.C:121
QUESO::MetropolisHastingsSG::readFullChain
void readFullChain(const std::string &inputFileName, const std::string &inputFileType, unsigned int chainSize, BaseVectorSequence< P_V, P_M > &workingChain)
This method reads the chain contents.
Definition: MetropolisHastingsSG.C:1445
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:52
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
QUESO::MHRawChainInfoStruct::mpiSum
void mpiSum(const MpiComm &comm, MHRawChainInfoStruct &sumInfo) const
Calculates the MPI sum of this.
Definition: MetropolisHastingsSG.C:121
QUESO::MHRawChainInfoStruct::copy
void copy(const MHRawChainInfoStruct &src)
Copies Metropolis-Hastings chain info from src to this.
Definition: MetropolisHastingsSG.C:101
QUESO::BaseVectorSequence::computeFilterParams
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
Definition: VectorSequence.C:641
QUESO::MarkovChainPositionData::logLikelihood
double logLikelihood() const
Logarithm of the value of the likelihood; access to private attribute m_logLikelihood.
Definition: MarkovChainPositionData.C:111
QUESO::MHRawChainInfoStruct::operator+=
MHRawChainInfoStruct & operator+=(const MHRawChainInfoStruct &rhs)
Addition assignment operator.
Definition: MetropolisHastingsSG.C:63
QUESO::SequenceOfVectors
Class for handling vector samples (sequence of vectors).
Definition: SequenceOfVectors.h:46
QUESO::BaseVectorSequence::unifiedPositionsOfMaximum
double unifiedPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
Definition: VectorSequence.C:433
UQ_MH_SG_FILENAME_FOR_NO_FILE
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
Definition: MetropolisHastingsSGOptions.h:35
QUESO::SequenceOfVectors::setPositionValues
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
Definition: SequenceOfVectors.C:221
QUESO::BaseVectorSequence::resizeSequence
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
QUESO::BaseVectorSequence::filter
virtual void filter(unsigned int initialPos, unsigned int spacing)=0
Filters positions in the sequence of vectors, starting at initialPos, and with spacing given by spaci...
QUESO::MetropolisHastingsSG::acceptAlpha
bool acceptAlpha(double alpha)
Decides whether or not to accept alpha.
Definition: MetropolisHastingsSG.C:869
QUESO::MarkovChainPositionData
A templated class that represents a Markov Chain.
Definition: MarkovChainPositionData.h:47
Generated on Thu Apr 23 2015 19:30:55 for queso-0.52.0 by  doxygen 1.8.5