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