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