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