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