queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
QUESO::MetropolisHastingsSG< P_V, P_M > Class Template Reference

A templated class that represents a Metropolis-Hastings generator of samples. More...

#include <MetropolisHastingsSG.h>

Public Member Functions

const BaseTKGroup< P_V, P_M > & transitionKernel () const
 Returns the underlying transition kernel for this sequence generator. More...
 
Constructor/Destructor methods
 MetropolisHastingsSG (const char *prefix, const MhOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, const P_M *inputProposalCovMatrix)
 Constructor. More...
 
 MetropolisHastingsSG (const char *prefix, const MhOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, double initialLogPrior, double initialLogLikelihood, const P_M *inputProposalCovMatrix)
 Constructor. More...
 
 MetropolisHastingsSG (const MLSamplingLevelOptions &mlOptions, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, const P_M *inputProposalCovMatrix)
 Constructor. More...
 
 MetropolisHastingsSG (const MLSamplingLevelOptions &mlOptions, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, double initialLogPrior, double initialLogLikelihood, const P_M *inputProposalCovMatrix)
 Constructor. More...
 
 ~MetropolisHastingsSG ()
 Destructor. More...
 
Statistical methods
void generateSequence (BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
 Method to generate the chain. More...
 
void getRawChainInfo (MHRawChainInfoStruct &info) const
 Gets information from the raw chain. More...
 

Private Member Functions

void commonConstructor ()
 Reads the options values from the options input file. More...
 
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. More...
 
void adapt (unsigned int positionId, BaseVectorSequence< P_V, P_M > &workingChain)
 Adaptive Metropolis method that deals with adapting the proposal covariance matrix. More...
 
bool delayedRejection (unsigned int positionId, const MarkovChainPositionData< P_V > &currentPositionData, MarkovChainPositionData< P_V > &currentCandidateData)
 Does delayed rejection. More...
 
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. More...
 
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. More...
 
double alpha (const std::vector< MarkovChainPositionData< P_V > * > &inputPositions, const std::vector< unsigned int > &inputTKStageIds)
 Calculates acceptance ratio. More...
 
bool acceptAlpha (double alpha)
 Decides whether or not to accept alpha. More...
 
int writeInfo (const BaseVectorSequence< P_V, P_M > &workingChain, std::ofstream &ofsvar) const
 Writes information about the Markov chain in a file. More...
 
void transformInitialCovMatrixToGaussianSpace (const BoxSubset< P_V, P_M > &boxSubset)
 

Private Attributes

const BaseEnvironmentm_env
 
const VectorSpace< P_V, P_M > & m_vectorSpace
 
const BaseJointPdf< P_V, P_M > & m_targetPdf
 
P_V m_initialPosition
 
P_M m_initialProposalCovMatrix
 
bool m_nullInputProposalCovMatrix
 
unsigned int m_numDisabledParameters
 
std::vector< bool > m_parameterEnabledStatus
 
ScopedPtr< const
ScalarFunctionSynchronizer
< P_V, P_M > >::Type 
m_targetPdfSynchronizer
 
SharedPtr< BaseTKGroup< P_V,
P_M > >::Type 
m_tk
 
SharedPtr< Algorithm< P_V, P_M >
>::Type 
m_algorithm
 
unsigned int m_positionIdForDebugging
 
unsigned int m_stageIdForDebugging
 
std::vector< unsigned int > m_idsOfUniquePositions
 
std::vector< double > m_logTargets
 
std::vector< double > m_alphaQuotients
 
double m_lastChainSize
 
ScopedPtr< P_V >::Type m_lastMean
 
ScopedPtr< P_M >::Type m_lastAdaptedCovMatrix
 
unsigned int m_numPositionsNotSubWritten
 
MHRawChainInfoStruct m_rawChainInfo
 
ScopedPtr< const
MhOptionsValues >::Type 
m_optionsObj
 
ScopedPtr
< MetropolisHastingsSGOptions >
::Type 
m_oldOptions
 
bool m_computeInitialPriorAndLikelihoodValues
 
double m_initialLogPriorValue
 
double m_initialLogLikelihoodValue
 
bool m_userDidNotProvideOptions
 
unsigned int m_latestDirtyCovMatrixIteration
 

I/O methods

void print (std::ostream &os) const
 TODO: Prints the sequence. More...
 
std::ostream & operator<< (std::ostream &os, const MetropolisHastingsSG< P_V, P_M > &obj)
 

Detailed Description

template<class P_V = GslVector, class P_M = GslMatrix>
class QUESO::MetropolisHastingsSG< P_V, P_M >

A templated class that represents a Metropolis-Hastings generator of samples.

This class implements a Metropolis-Hastings generator of samples. 'SG' stands for 'Sequence Generator'. Options reading is handled by class 'MetropolisHastingsOptions'. If options request data to be written in the output file (MATLAB .m format only, for now), the user can check which MATLAB variables are defined and set by running 'grep zeros <OUTPUT file="" name>="">' after the solution procedures ends. The names of the variables are self explanatory.

Definition at line 126 of file MetropolisHastingsSG.h.

Constructor & Destructor Documentation

template<class P_V , class P_M >
QUESO::MetropolisHastingsSG< P_V, P_M >::MetropolisHastingsSG ( const char *  prefix,
const MhOptionsValues alternativeOptionsValues,
const BaseVectorRV< P_V, P_M > &  sourceRv,
const P_V &  initialPosition,
const P_M *  inputProposalCovMatrix 
)

Constructor.

This method reads the options from the options input file. It calls commonConstructor().

Requirements:

  1. the image set of the vector random variable 'sourceRv' should belong to a vector space of dimension equal to the size of the vector 'initialPosition' and
  2. if 'inputProposalCovMatrix' is not NULL, is should be square and its size should be equal to the size of 'initialPosition'.

If alternativeOptionsValues is NULL and an input file is specified, the constructor reads input options that begin with the string '<prefix>mh_'. For instance, if 'prefix' is 'pROblem_775_ip_', then the constructor will read all options that begin with 'pROblem_775_ip_mh_'.

If alternativeOptionsValues is not NULL, the input file is ignored and construction copies the object pointed to by alternativeOptionsValues to and stores the copy internally. Users may delete the object poined to by alternativeOptionsValues. Users cannot change the options object after MetropolisHastingsSG has been constructed.

Options reading is handled by class 'MetropolisHastingsOptions'.

Parameters
prefixPrefix
alternativeOptionsValuesOptions (if no input file)
sourceRvThe source RV
initialPositionInitial chain position
inputProposalCovMatrixProposal cov. matrix

Definition at line 141 of file MetropolisHastingsSG.C.

References QUESO::queso_require_equal_to_msg.

147  :
148  m_env (sourceRv.env()),
149  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
150  m_targetPdf (sourceRv.pdf()),
151  m_initialPosition (initialPosition),
153  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
154  m_numDisabledParameters (0), // gpmsa2
156  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
157  m_tk (),
158  m_algorithm (),
161  m_idsOfUniquePositions (0),//0.),
162  m_logTargets (0),//0.),
163  m_alphaQuotients (0),//0.),
164  m_lastChainSize (0),
165  m_lastMean (),
168  m_optionsObj (),
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 
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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:90
const BaseJointPdf< P_V, P_M > & m_targetPdf
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
std::vector< bool > m_parameterEnabledStatus
ScopedPtr< P_M >::Type m_lastAdaptedCovMatrix
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:83
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:354
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
std::vector< double > m_alphaQuotients
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
std::vector< unsigned int > m_idsOfUniquePositions
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
std::vector< double > m_logTargets
void commonConstructor()
Reads the options values from the options input file.
ScopedPtr< P_V >::Type m_lastMean
unsigned int dimLocal() const
Definition: VectorSpace.C:155
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"))
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:76
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
QUESO::MetropolisHastingsSG< P_V, P_M >::MetropolisHastingsSG ( const char *  prefix,
const MhOptionsValues alternativeOptionsValues,
const BaseVectorRV< P_V, P_M > &  sourceRv,
const P_V &  initialPosition,
double  initialLogPrior,
double  initialLogLikelihood,
const P_M *  inputProposalCovMatrix 
)

Constructor.

Parameters
prefixPrefix
alternativeOptionsValuesOptions (if no input file)
sourceRvThe source RV
initialPositionInitial chain position
initialLogPriorInitial chain prior
initialLogLikelihoodInitial chain likelihood
inputProposalCovMatrixProposal cov. matrix

Definition at line 221 of file MetropolisHastingsSG.C.

References QUESO::queso_require_equal_to_msg.

229  :
230  m_env (sourceRv.env()),
231  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
232  m_targetPdf (sourceRv.pdf()),
233  m_initialPosition (initialPosition),
235  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
236  m_numDisabledParameters (0), // gpmsa2
238  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
239  m_tk (),
240  m_algorithm (),
243  m_idsOfUniquePositions (0),//0.),
244  m_logTargets (0),//0.),
245  m_alphaQuotients (0),//0.),
246  m_lastChainSize (0),
247  m_lastMean (),
250  m_optionsObj (),
252  m_initialLogPriorValue (initialLogPrior),
253  m_initialLogLikelihoodValue (initialLogLikelihood),
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 
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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:90
const BaseJointPdf< P_V, P_M > & m_targetPdf
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
std::vector< bool > m_parameterEnabledStatus
ScopedPtr< P_M >::Type m_lastAdaptedCovMatrix
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:83
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:354
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
std::vector< double > m_alphaQuotients
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
std::vector< unsigned int > m_idsOfUniquePositions
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
std::vector< double > m_logTargets
void commonConstructor()
Reads the options values from the options input file.
ScopedPtr< P_V >::Type m_lastMean
unsigned int dimLocal() const
Definition: VectorSpace.C:155
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"))
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:76
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
QUESO::MetropolisHastingsSG< P_V, P_M >::MetropolisHastingsSG ( const MLSamplingLevelOptions mlOptions,
const BaseVectorRV< P_V, P_M > &  sourceRv,
const P_V &  initialPosition,
const P_M *  inputProposalCovMatrix 
)

Constructor.

Definition at line 303 of file MetropolisHastingsSG.C.

References QUESO::MetropolisHastingsSG< P_V, P_M >::commonConstructor(), QUESO::MetropolisHastingsSG< P_V, P_M >::m_env, QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialProposalCovMatrix, QUESO::MetropolisHastingsSG< P_V, P_M >::m_oldOptions, QUESO::MetropolisHastingsSG< P_V, P_M >::m_optionsObj, and QUESO::BaseEnvironment::subDisplayFile().

308  :
309  m_env (sourceRv.env()),
310  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
311  m_targetPdf (sourceRv.pdf()),
312  m_initialPosition (initialPosition),
314  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
315  m_numDisabledParameters (0), // gpmsa2
317  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
318  m_tk (),
319  m_algorithm (),
322  m_idsOfUniquePositions (0),//0.),
323  m_logTargets (0),//0.),
324  m_alphaQuotients (0),//0.),
325  m_lastChainSize (0),
326  m_lastMean (),
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 
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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:90
ScopedPtr< MetropolisHastingsSGOptions >::Type m_oldOptions
const BaseJointPdf< P_V, P_M > & m_targetPdf
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
std::vector< bool > m_parameterEnabledStatus
ScopedPtr< P_M >::Type m_lastAdaptedCovMatrix
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:83
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
std::vector< double > m_alphaQuotients
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
std::vector< unsigned int > m_idsOfUniquePositions
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
std::vector< double > m_logTargets
void commonConstructor()
Reads the options values from the options input file.
ScopedPtr< P_V >::Type m_lastMean
unsigned int dimLocal() const
Definition: VectorSpace.C:155
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:76
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
QUESO::MetropolisHastingsSG< P_V, P_M >::MetropolisHastingsSG ( const MLSamplingLevelOptions mlOptions,
const BaseVectorRV< P_V, P_M > &  sourceRv,
const P_V &  initialPosition,
double  initialLogPrior,
double  initialLogLikelihood,
const P_M *  inputProposalCovMatrix 
)

Constructor.

Definition at line 369 of file MetropolisHastingsSG.C.

References QUESO::MetropolisHastingsSG< P_V, P_M >::commonConstructor(), QUESO::MetropolisHastingsSG< P_V, P_M >::m_env, QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialProposalCovMatrix, QUESO::MetropolisHastingsSG< P_V, P_M >::m_oldOptions, QUESO::MetropolisHastingsSG< P_V, P_M >::m_optionsObj, and QUESO::BaseEnvironment::subDisplayFile().

376  :
377  m_env (sourceRv.env()),
378  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
379  m_targetPdf (sourceRv.pdf()),
380  m_initialPosition (initialPosition),
382  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
383  m_numDisabledParameters (0), // gpmsa2
385  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
386  m_tk (),
387  m_algorithm (),
390  m_idsOfUniquePositions (0),//0.),
391  m_logTargets (0),//0.),
392  m_alphaQuotients (0),//0.),
393  m_lastChainSize (0),
394  m_lastMean (),
397  m_initialLogPriorValue (initialLogPrior),
398  m_initialLogLikelihoodValue (initialLogLikelihood),
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 
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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:90
ScopedPtr< MetropolisHastingsSGOptions >::Type m_oldOptions
const BaseJointPdf< P_V, P_M > & m_targetPdf
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
std::vector< bool > m_parameterEnabledStatus
ScopedPtr< P_M >::Type m_lastAdaptedCovMatrix
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:83
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
std::vector< double > m_alphaQuotients
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
std::vector< unsigned int > m_idsOfUniquePositions
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
std::vector< double > m_logTargets
void commonConstructor()
Reads the options values from the options input file.
ScopedPtr< P_V >::Type m_lastMean
unsigned int dimLocal() const
Definition: VectorSpace.C:155
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:76
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
QUESO::MetropolisHastingsSG< P_V, P_M >::~MetropolisHastingsSG ( )

Destructor.

Definition at line 437 of file MetropolisHastingsSG.C.

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;
446  m_alphaQuotients.clear();
447  m_logTargets.clear();
448  m_numDisabledParameters = 0; // gpmsa2
449  m_parameterEnabledStatus.clear(); // gpmsa2
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 }
void reset()
Resets Metropolis-Hastings chain info.
std::vector< bool > m_parameterEnabledStatus
std::vector< double > m_alphaQuotients
std::vector< unsigned int > m_idsOfUniquePositions
std::vector< double > m_logTargets
MHRawChainInfoStruct m_rawChainInfo

Member Function Documentation

template<class P_V , class P_M >
bool QUESO::MetropolisHastingsSG< P_V, P_M >::acceptAlpha ( double  alpha)
private

Decides whether or not to accept alpha.

If either alpha is negative or greater than one, its value will not be accepted.

Definition at line 753 of file MetropolisHastingsSG.C.

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 }
const BaseEnvironment & m_env
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:471
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
double alpha(const std::vector< MarkovChainPositionData< P_V > * > &inputPositions, const std::vector< unsigned int > &inputTKStageIds)
Calculates acceptance ratio.
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::adapt ( unsigned int  positionId,
BaseVectorSequence< P_V, P_M > &  workingChain 
)
private

Adaptive Metropolis method that deals with adapting the proposal covariance matrix.

Definition at line 1977 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::getPositionValues(), QUESO::MiscGetEllapsedSeconds(), QUESO::queso_require_equal_to_msg, QUESO::SequenceOfVectors< V, M >::resizeSequence(), QUESO::SequenceOfVectors< V, M >::setPositionValues(), QUESO::SequenceOfVectors< V, M >::subSequenceSize(), QUESO::UQ_INCOMPLETE_IMPLEMENTATION_RC, QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC, and QUESO::UQ_OK_RC.

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);
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.
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) {
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,
2081  *m_lastMean,
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) {
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) {
2274  }
2275 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:342
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseEnvironment & m_env
const int UQ_OK_RC
Definition: Defines.h:92
M * newMatrix() const
Creates an empty matrix of size given by Map&amp; map. See template specialization.
std::vector< bool > m_parameterEnabledStatus
ScopedPtr< P_M >::Type m_lastAdaptedCovMatrix
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Definition: Defines.h:100
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
V * newVector() const
Creates an empty vector of size given by Map&amp; map. See template specialization.
const int UQ_INCOMPLETE_IMPLEMENTATION_RC
Definition: Defines.h:93
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.
ScopedPtr< P_V >::Type m_lastMean
unsigned int dimLocal() const
Definition: VectorSpace.C:155
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"))
unsigned int displayVerbosity() const
Definition: Environment.C:450
double MiscGetEllapsedSeconds(struct timeval *timeval0)
M * newDiagMatrix(const V &v) const
Creates a diagonal matrix with the elements and size of vector v.
Definition: VectorSpace.C:190
MHRawChainInfoStruct m_rawChainInfo
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
double QUESO::MetropolisHastingsSG< P_V, P_M >::alpha ( const std::vector< MarkovChainPositionData< P_V > * > &  inputPositions,
const std::vector< unsigned int > &  inputTKStageIds 
)
private

Calculates acceptance ratio.

The acceptance ratio is used to decide whether to accept or reject a candidate.

Definition at line 561 of file MetropolisHastingsSG.C.

References QUESO::queso_isnan(), QUESO::queso_require_equal_to_msg, and QUESO::BaseEnvironment::worldRank().

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 }
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:342
const BaseEnvironment & m_env
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
bool queso_isnan(T arg)
Definition: math_macros.h:39
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
int fullRank() const
Returns the rank of the MPI process in QUESO&#39;s full communicator.
Definition: Environment.C:268
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
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"))
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
unsigned int displayVerbosity() const
Definition: Environment.C:450
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
double alpha(const std::vector< MarkovChainPositionData< P_V > * > &inputPositions, const std::vector< unsigned int > &inputTKStageIds)
Calculates acceptance ratio.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::commonConstructor ( )
private

Reads the options values from the options input file.

This method actually reads the options input file, such as the value for the Delayed Rejection scales, the presence of Hessian covariance matrices and reads the user-provided initial proposal covariance matrix (alternative to local Hessians).

Definition at line 463 of file MetropolisHastingsSG.C.

References QUESO::Factory< BaseTKGroup< GslVector, GslMatrix > >::build(), QUESO::Factory< Algorithm< GslVector, GslMatrix > >::build(), QUESO::TransitionKernelFactory::set_dr_scales(), QUESO::AlgorithmFactory::set_environment(), QUESO::TransitionKernelFactory::set_initial_cov_matrix(), QUESO::TransitionKernelFactory::set_options(), QUESO::TransitionKernelFactory::set_pdf_synchronizer(), QUESO::TransitionKernelFactory::set_target_pdf(), QUESO::AlgorithmFactory::set_tk(), and QUESO::TransitionKernelFactory::set_vectorspace().

Referenced by QUESO::MetropolisHastingsSG< P_V, P_M >::MetropolisHastingsSG().

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()) {
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
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 
550 
551  // This instantiates all the algorithms with their associated factories
552  AlgorithmFactoryInitializer algorithm_factory_initializer;
553 
557 }
void transformInitialCovMatrixToGaussianSpace(const BoxSubset< P_V, P_M > &boxSubset)
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:342
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseJointPdf< P_V, P_M > & m_targetPdf
static void set_tk(const BaseTKGroup< GslVector, GslMatrix > &tk)
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
std::vector< bool > m_parameterEnabledStatus
static void set_pdf_synchronizer(const ScalarFunctionSynchronizer< GslVector, GslMatrix > &synchronizer)
static void set_vectorspace(const VectorSpace< GslVector, GslMatrix > &v)
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
static void set_target_pdf(const BaseJointPdf< GslVector, GslMatrix > &target_pdf)
static void set_options(const MhOptionsValues &options)
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
static SharedPtr< BaseTKGroup< GslVector, GslMatrix > >::Type build(const std::string &name)
unsigned int dimLocal() const
Definition: VectorSpace.C:155
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
static void set_environment(const BaseEnvironment &env)
static void set_dr_scales(const std::vector< double > &scales)
const std::string & subIdString() const
Access to the attribute m_subIdString; which stores the string for the sub-environment, and it will be used, for instance, to create the output files for each sub-environment.
Definition: Environment.C:348
static void set_initial_cov_matrix(GslMatrix &cov_matrix)
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
bool QUESO::MetropolisHastingsSG< P_V, P_M >::delayedRejection ( unsigned int  positionId,
const MarkovChainPositionData< P_V > &  currentPositionData,
MarkovChainPositionData< P_V > &  currentCandidateData 
)
private

Does delayed rejection.

When faced with an imminent rejection, this method computes a series of new candidates in the same proposal direction but with smaller proposal step sizes. For each of these new candidates, we check if it will be accepted, if not we repeat the process until all the candidates have been tested.

If there is a candidate that will be accepted, this method will return true, otherwise it returns false.

currentCandidateData is updated whenever a new proposal is generated throughout the delayed rejection procedure.

delayedRejection promises not to change currentPositionData, because changing the current position of the Markov chain would be grossly inappropriate.

Definition at line 2279 of file MetropolisHastingsSG.C.

References QUESO::if, QUESO::MiscGetEllapsedSeconds(), QUESO::queso_require_equal_to_msg, QUESO::MarkovChainPositionData< V >::set(), QUESO::UQ_OK_RC, and QUESO::MarkovChainPositionData< V >::vecValues().

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  }
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);
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 }
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseJointPdf< P_V, P_M > & m_targetPdf
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
const int UQ_OK_RC
Definition: Defines.h:92
std::vector< bool > m_parameterEnabledStatus
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
unsigned int dimLocal() const
Definition: VectorSpace.C:155
McOptionsValues::McOptionsValues(#ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS const SsOptionsValues *alternativePSsOptionsValues, const SsOptionsValues *alternativeQSsOptionsValues#endif if)(m_alternativeQSsOptionsValues=*alternativeQSsOptionsValues)
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"))
unsigned int displayVerbosity() const
Definition: Environment.C:450
double MiscGetEllapsedSeconds(struct timeval *timeval0)
double alpha(const std::vector< MarkovChainPositionData< P_V > * > &inputPositions, const std::vector< unsigned int > &inputTKStageIds)
Calculates acceptance ratio.
bool acceptAlpha(double alpha)
Decides whether or not to accept alpha.
MHRawChainInfoStruct m_rawChainInfo
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::generateFullChain ( const P_V &  valuesOf1stPosition,
unsigned int  chainSize,
BaseVectorSequence< P_V, P_M > &  workingChain,
ScalarSequence< double > *  workingLogLikelihoodValues,
ScalarSequence< double > *  workingLogTargetValues 
)
private

This method generates the chain.

Given the size of the chain, the values of the first position in the chain. It checks if the position is out of target pdf support, otherwise it calculates the value of both its likelihood and prior and adds the position to the chain. For the next positions, once they are generated, some tests are performed (such as unicity and the value of alpha) and the steps for the first position are repeated, including the optional Delayed Rejection and the adaptive Metropolis (adaptation of covariance matrix) steps.

Definition at line 1332 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::estimateConvBrooksGelman(), QUESO::if, QUESO::MarkovChainPositionData< V >::logLikelihood(), QUESO::MarkovChainPositionData< V >::logTarget(), QUESO::MiscGetEllapsedSeconds(), QUESO::BaseVectorSequence< V, M >::name(), QUESO::queso_require_equal_to_msg, QUESO::ScalarSequence< T >::resizeSequence(), QUESO::BaseVectorSequence< V, M >::resizeSequence(), QUESO::BaseVectorSequence< V, M >::setPositionValues(), QUESO::BaseVectorSequence< V, M >::subSequenceSize(), QUESO::BaseVectorSequence< V, M >::subWriteContents(), QUESO::ScalarSequence< T >::subWriteContents(), QUESO::UQ_OK_RC, and QUESO::MarkovChainPositionData< V >::vecValues().

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 
1357 
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.;
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) {
1418  }
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);
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());
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 
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());
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) {
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);
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());
1774  }
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 
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---------------------- --------------";
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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
const VectorSpace< P_V, P_M > & m_vectorSpace
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
const BaseJointPdf< P_V, P_M > & m_targetPdf
void reset()
Resets Metropolis-Hastings chain info.
ScopedPtr< const ScalarFunctionSynchronizer< P_V, P_M > >::Type m_targetPdfSynchronizer
const BaseEnvironment & m_env
void adapt(unsigned int positionId, BaseVectorSequence< P_V, P_M > &workingChain)
Adaptive Metropolis method that deals with adapting the proposal covariance matrix.
const int UQ_OK_RC
Definition: Defines.h:92
std::vector< bool > m_parameterEnabledStatus
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
bool delayedRejection(unsigned int positionId, const MarkovChainPositionData< P_V > &currentPositionData, MarkovChainPositionData< P_V > &currentCandidateData)
Does delayed rejection.
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
std::vector< double > m_alphaQuotients
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
std::vector< unsigned int > m_idsOfUniquePositions
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
std::vector< double > m_logTargets
unsigned int dimLocal() const
Definition: VectorSpace.C:155
McOptionsValues::McOptionsValues(#ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS const SsOptionsValues *alternativePSsOptionsValues, const SsOptionsValues *alternativeQSsOptionsValues#endif if)(m_alternativeQSsOptionsValues=*alternativeQSsOptionsValues)
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"))
SharedPtr< Algorithm< P_V, P_M > >::Type m_algorithm
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.
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:133
unsigned int displayVerbosity() const
Definition: Environment.C:450
const MpiComm & fullComm() const
Access function for the communicator that was passed to QUESO&#39;s environment.
Definition: Environment.C:274
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:335
double MiscGetEllapsedSeconds(struct timeval *timeval0)
bool acceptAlpha(double alpha)
Decides whether or not to accept alpha.
MHRawChainInfoStruct m_rawChainInfo
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::generateSequence ( BaseVectorSequence< P_V, P_M > &  workingChain,
ScalarSequence< double > *  workingLogLikelihoodValues,
ScalarSequence< double > *  workingLogTargetValues 
)

Method to generate the chain.

Requirement: the vector space 'm_vectorSpace' should have dimension equal to the size of a vector in 'workingChain'. If the requirement is satisfied, this operation sets the size and the contents of 'workingChain' using the algorithm options set in the constructor. If not NULL, 'workingLogLikelihoodValues' and 'workingLogTargetValues' are set accordingly. This operation currently implements the DRAM algorithm (Heikki Haario, Marko Laine, Antonietta Mira and Eero Saksman, "DRAM: Efficient Adaptive MCMC", Statistics and Computing (2006), 16:339-354), as a translation of the core routine at the MCMC toolbox for MATLAB, available at helios.fmi.fi/~lainema/mcmc/‎ (accessed in July 3rd, 2013). Indeed, the example available in examples/statisticalInverseProblem is closely related to the 'normal example' in the toolbox. Over time, though:

  1. the whole set of QUESO classes took shape, focusing not only on Markov Chains, but on statistical forward problems and model validation as well;
  2. the interfaces to this Metropolis-Hastings class changed;
  3. QUESO has parallel capabilities;
  4. TK (transition kernel) class has been added in order to have both DRAM with Stochastic Newton capabilities.

Definition at line 859 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::computeFilterParams(), QUESO::BaseVectorSequence< V, M >::computeStatistics(), QUESO::BaseVectorSequence< V, M >::filter(), QUESO::ScalarSequence< T >::filter(), QUESO::SequenceOfVectors< V, M >::getPositionValues(), QUESO::BaseVectorSequence< V, M >::name(), QUESO::FilePtrSetStruct::ofsVar, QUESO::queso_require_not_equal_to_msg, QUESO::BaseVectorSequence< V, M >::setName(), QUESO::BaseVectorSequence< V, M >::subPositionsOfMaximum(), QUESO::BaseVectorSequence< V, M >::subSequenceSize(), QUESO::SequenceOfVectors< V, M >::subSequenceSize(), QUESO::BaseVectorSequence< V, M >::subWriteContents(), QUESO::ScalarSequence< T >::subWriteContents(), QUESO::BaseVectorSequence< V, M >::unifiedPositionsOfMaximum(), QUESO::BaseVectorSequence< V, M >::unifiedWriteContents(), QUESO::ScalarSequence< T >::unifiedWriteContents(), QUESO::UQ_OK_RC, and QUESO::BaseVectorSequence< V, M >::vectorSizeLocal().

Referenced by QUESO::MLSampling< P_V, P_M >::generateBalLinkedChains_all(), and QUESO::MLSampling< P_V, P_M >::generateUnbLinkedChains_all().

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,
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,
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,
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,
1005  m_optionsObj->m_rawChainDataOutputFileName + "_logtarget",
1006  m_optionsObj->m_rawChainDataOutputFileType,
1007  m_optionsObj->m_rawChainDataOutputAllowedSet);
1008  }
1009 
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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:342
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseEnvironment & m_env
const int UQ_OK_RC
Definition: Defines.h:92
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"))
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
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.
int writeInfo(const BaseVectorSequence< P_V, P_M > &workingChain, std::ofstream &ofsvar) const
Writes information about the Markov chain in a file.
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
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.
unsigned int dimLocal() const
Definition: VectorSpace.C:155
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.
unsigned int displayVerbosity() const
Definition: Environment.C:450
bool openOutputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, bool writeOver, FilePtrSetStruct &filePtrSet) const
Opens an output file for each sub-environment that was chosen to send data to the file...
Definition: Environment.C:521
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::getRawChainInfo ( MHRawChainInfoStruct info) const

Gets information from the raw chain.

Definition at line 1312 of file MetropolisHastingsSG.C.

Referenced by QUESO::MLSampling< P_V, P_M >::generateBalLinkedChains_all(), and QUESO::MLSampling< P_V, P_M >::generateUnbLinkedChains_all().

1313 {
1314  info = m_rawChainInfo;
1315  return;
1316 }
MHRawChainInfoStruct m_rawChainInfo
template<class P_V = GslVector, class P_M = GslMatrix>
void QUESO::MetropolisHastingsSG< P_V, P_M >::print ( std::ostream &  os) const

TODO: Prints the sequence.

Todo:
: implement me!
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::readFullChain ( const std::string &  inputFileName,
const std::string &  inputFileType,
unsigned int  chainSize,
BaseVectorSequence< P_V, P_M > &  workingChain 
)
private

This method reads the chain contents.

Definition at line 1320 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::unifiedReadContents().

1325 {
1326  workingChain.unifiedReadContents(inputFileName,inputFileType,chainSize);
1327  return;
1328 }
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::transformInitialCovMatrixToGaussianSpace ( const BoxSubset< P_V, P_M > &  boxSubset)
private

Definition at line 2532 of file MetropolisHastingsSG.C.

References QUESO::BoxSubset< V, M >::maxValues(), QUESO::BoxSubset< V, M >::minValues(), and QUESO::queso_isfinite().

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 "
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 }
bool queso_isfinite(T arg)
Definition: math_macros.h:51
template<class P_V , class P_M >
const BaseTKGroup< P_V, P_M > & QUESO::MetropolisHastingsSG< P_V, P_M >::transitionKernel ( ) const

Returns the underlying transition kernel for this sequence generator.

Definition at line 845 of file MetropolisHastingsSG.C.

846 {
847  return *m_tk;
848 }
SharedPtr< BaseTKGroup< P_V, P_M > >::Type m_tk
template<class P_V , class P_M >
void QUESO::MetropolisHastingsSG< P_V, P_M >::updateAdaptedCovMatrix ( const BaseVectorSequence< P_V, P_M > &  subChain,
unsigned int  idOfFirstPositionInSubChain,
double &  lastChainSize,
P_V &  lastMean,
P_M &  lastAdaptedCovMatrix 
)
private

This method updates the adapted covariance matrix.

This function is called is the option to used adaptive Metropolis was chosen by the user (via options input file). It performs an adaptation of covariance matrix.

Definition at line 2469 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::getPositionValues(), QUESO::matrixProduct(), QUESO::BaseVectorSequence< V, M >::subMeanExtra(), QUESO::BaseVectorSequence< V, M >::subMeanPlain(), and QUESO::BaseVectorSequence< V, M >::subSequenceSize().

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 }
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:174
const VectorSpace< P_V, P_M > & m_vectorSpace
std::vector< bool > m_parameterEnabledStatus
unsigned int dimLocal() const
Definition: VectorSpace.C:155
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2036
template<class P_V , class P_M >
int QUESO::MetropolisHastingsSG< P_V, P_M >::writeInfo ( const BaseVectorSequence< P_V, P_M > &  workingChain,
std::ofstream &  ofsvar 
) const
private

Writes information about the Markov chain in a file.

It writes down the alpha quotients, the number of rejected positions, number of positions out of target support, the name of the components and the chain runtime.

Definition at line 767 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::name(), QUESO::BaseVectorSequence< V, M >::subSequenceSize(), and QUESO::UQ_OK_RC.

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 }
void printComponentsNames(std::ostream &os, bool printHorizontally) const
Prints the local component names.
Definition: VectorSpace.C:289
const VectorSpace< P_V, P_M > & m_vectorSpace
const BaseEnvironment & m_env
const int UQ_OK_RC
Definition: Defines.h:92
std::vector< double > m_alphaQuotients
ScopedPtr< const MhOptionsValues >::Type m_optionsObj
std::vector< double > m_logTargets
const std::string & subIdString() const
Access to the attribute m_subIdString; which stores the string for the sub-environment, and it will be used, for instance, to create the output files for each sub-environment.
Definition: Environment.C:348
MHRawChainInfoStruct m_rawChainInfo
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

Friends And Related Function Documentation

template<class P_V = GslVector, class P_M = GslMatrix>
std::ostream& operator<< ( std::ostream &  os,
const MetropolisHastingsSG< P_V, P_M > &  obj 
)
friend

Definition at line 228 of file MetropolisHastingsSG.h.

230  {
231  obj.print(os);
232 
233  return os;
234  }

Member Data Documentation

template<class P_V = GslVector, class P_M = GslMatrix>
SharedPtr<Algorithm<P_V, P_M> >::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_algorithm
private

Definition at line 324 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
std::vector<double> QUESO::MetropolisHastingsSG< P_V, P_M >::m_alphaQuotients
private

Definition at line 329 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
bool QUESO::MetropolisHastingsSG< P_V, P_M >::m_computeInitialPriorAndLikelihoodValues
private

Definition at line 340 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
const BaseEnvironment& QUESO::MetropolisHastingsSG< P_V, P_M >::m_env
private
template<class P_V = GslVector, class P_M = GslMatrix>
std::vector<unsigned int> QUESO::MetropolisHastingsSG< P_V, P_M >::m_idsOfUniquePositions
private

Definition at line 327 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
double QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialLogLikelihoodValue
private

Definition at line 342 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
double QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialLogPriorValue
private

Definition at line 341 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
P_V QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialPosition
private

Definition at line 316 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
P_M QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialProposalCovMatrix
private
template<class P_V = GslVector, class P_M = GslMatrix>
ScopedPtr<P_M>::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_lastAdaptedCovMatrix
private

Definition at line 332 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
double QUESO::MetropolisHastingsSG< P_V, P_M >::m_lastChainSize
private

Definition at line 330 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
ScopedPtr<P_V>::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_lastMean
private

Definition at line 331 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
unsigned int QUESO::MetropolisHastingsSG< P_V, P_M >::m_latestDirtyCovMatrixIteration
private

Definition at line 349 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
std::vector<double> QUESO::MetropolisHastingsSG< P_V, P_M >::m_logTargets
private

Definition at line 328 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
bool QUESO::MetropolisHastingsSG< P_V, P_M >::m_nullInputProposalCovMatrix
private

Definition at line 318 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
unsigned int QUESO::MetropolisHastingsSG< P_V, P_M >::m_numDisabledParameters
private

Definition at line 319 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
unsigned int QUESO::MetropolisHastingsSG< P_V, P_M >::m_numPositionsNotSubWritten
private

Definition at line 333 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
ScopedPtr<MetropolisHastingsSGOptions>::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_oldOptions
private
template<class P_V = GslVector, class P_M = GslMatrix>
ScopedPtr<const MhOptionsValues>::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_optionsObj
private
template<class P_V = GslVector, class P_M = GslMatrix>
std::vector<bool> QUESO::MetropolisHastingsSG< P_V, P_M >::m_parameterEnabledStatus
private

Definition at line 320 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
unsigned int QUESO::MetropolisHastingsSG< P_V, P_M >::m_positionIdForDebugging
private

Definition at line 325 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
MHRawChainInfoStruct QUESO::MetropolisHastingsSG< P_V, P_M >::m_rawChainInfo
private

Definition at line 335 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
unsigned int QUESO::MetropolisHastingsSG< P_V, P_M >::m_stageIdForDebugging
private

Definition at line 326 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
const BaseJointPdf<P_V,P_M>& QUESO::MetropolisHastingsSG< P_V, P_M >::m_targetPdf
private

Definition at line 315 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
ScopedPtr<const ScalarFunctionSynchronizer<P_V,P_M> >::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_targetPdfSynchronizer
private

Definition at line 321 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
SharedPtr<BaseTKGroup<P_V,P_M> >::Type QUESO::MetropolisHastingsSG< P_V, P_M >::m_tk
private

Definition at line 323 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
bool QUESO::MetropolisHastingsSG< P_V, P_M >::m_userDidNotProvideOptions
private

Definition at line 347 of file MetropolisHastingsSG.h.

template<class P_V = GslVector, class P_M = GslMatrix>
const VectorSpace<P_V,P_M>& QUESO::MetropolisHastingsSG< P_V, P_M >::m_vectorSpace
private

Definition at line 314 of file MetropolisHastingsSG.h.


The documentation for this class was generated from the following files:

Generated on Tue Jun 5 2018 19:49:34 for queso-0.57.1 by  doxygen 1.8.5