queso-0.50.1
Private Member Functions | Private Attributes | List of all members
QUESO::MetropolisHastingsSG< P_V, P_M > Class Template Reference

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

#include <MetropolisHastingsSG.h>

Collaboration diagram for QUESO::MetropolisHastingsSG< P_V, P_M >:
Collaboration graph
[legend]

Public Member Functions

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 MLSamplingLevelOptions &mlOptions, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, 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 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 MarkovChainPositionData< P_V > &x, const MarkovChainPositionData< P_V > &y, unsigned int xStageId, unsigned int yStageId, double *alphaQuotientPtr=NULL)
 Calculates acceptance ration. More...
 
double alpha (const std::vector< MarkovChainPositionData< P_V > * > &inputPositions, const std::vector< unsigned int > &inputTKStageIds)
 Calculates acceptance ration. 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...
 

Private Attributes

const BaseEnvironment & m_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
 
const
ScalarFunctionSynchronizer
< P_V, P_M > * 		m_targetPdfSynchronizer
 
BaseTKGroup< P_V, P_M > * m_tk
 
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
 
P_V * m_lastMean
 
P_M * m_lastAdaptedCovMatrix
 
unsigned int m_numPositionsNotSubWritten
 
MHRawChainInfoStruct m_rawChainInfo
 
MhOptionsValues m_alternativeOptionsValues
 
MetropolisHastingsSGOptions * m_optionsObj
 

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, class P_M>
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 120 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 the requirements are satisfied, the constructor then 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_'. 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 130 of file MetropolisHastingsSG.C.

References UQ_FATAL_TEST_MACRO.

136  :
137  m_env (sourceRv.env()),
138  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
139  m_targetPdf (sourceRv.pdf()),
140  m_initialPosition (initialPosition),
141  m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
142  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
143  m_numDisabledParameters (0), // gpmsa2
144  m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true), // gpmsa2
145  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
146  m_tk (NULL),
147  m_positionIdForDebugging (0),
148  m_stageIdForDebugging (0),
149  m_idsOfUniquePositions (0),//0.),
150  m_logTargets (0),//0.),
151  m_alphaQuotients (0),//0.),
152  m_lastChainSize (0),
153  m_lastMean (NULL),
154  m_lastAdaptedCovMatrix (NULL),
155  m_numPositionsNotSubWritten (0),
156 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
157  m_alternativeOptionsValues (NULL,NULL),
158 #else
159  m_alternativeOptionsValues (),
160 #endif
161  m_optionsObj (NULL)
162 {
163  if (inputProposalCovMatrix != NULL) {
164  m_initialProposalCovMatrix = *inputProposalCovMatrix;
165  }
166  if (alternativeOptionsValues) m_alternativeOptionsValues = *alternativeOptionsValues;
167  if (m_env.optionsInputFileName() == "") {
168  m_optionsObj = new MetropolisHastingsSGOptions(m_env,prefix,m_alternativeOptionsValues);
169  }
170  else {
171  m_optionsObj = new MetropolisHastingsSGOptions(m_env,prefix);
172  m_optionsObj->scanOptionsValues();
173  }
174 
175  if ((m_env.subDisplayFile() ) &&
176  (m_optionsObj->m_ov.m_totallyMute == false)) {
177  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::constructor(1)"
178  << ": prefix = " << prefix
179  << ", alternativeOptionsValues = " << alternativeOptionsValues
180  << ", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
181  << ", m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
182  << std::endl;
183  }
184 
185  UQ_FATAL_TEST_MACRO(sourceRv.imageSet().vectorSpace().dimLocal() != initialPosition.sizeLocal(),
186  m_env.worldRank(),
187  "MetropolisHastingsSG<P_V,P_M>::constructor(1)",
188  "'sourceRv' and 'initialPosition' should have equal dimensions");
189 
190  if (inputProposalCovMatrix) {
191  UQ_FATAL_TEST_MACRO(sourceRv.imageSet().vectorSpace().dimLocal() != inputProposalCovMatrix->numRowsLocal(),
192  m_env.worldRank(),
193  "MetropolisHastingsSG<P_V,P_M>::constructor(1)",
194  "'sourceRv' and 'inputProposalCovMatrix' should have equal dimensions");
195  UQ_FATAL_TEST_MACRO(inputProposalCovMatrix->numCols() != inputProposalCovMatrix->numRowsGlobal(),
196  m_env.worldRank(),
197  "MetropolisHastingsSG<P_V,P_M>::constructor(1)",
198  "'inputProposalCovMatrix' should be a square matrix");
199  }
200 
201  commonConstructor();
202 
203  if ((m_env.subDisplayFile() ) &&
204  (m_optionsObj->m_ov.m_totallyMute == false)) {
205  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::constructor(1)"
206  << std::endl;
207  }
208 }
QUESO::MetropolisHastingsSG::m_logTargets
std::vector< double > m_logTargets
Definition: MetropolisHastingsSG.h:268
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::BaseVectorRV::pdf
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:86
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::MetropolisHastingsSGOptions::scanOptionsValues
void scanOptionsValues()
It scans the option values from the options input file.
Definition: MetropolisHastingsSGOptions.C:519
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:218
QUESO::MetropolisHastingsSG::m_targetPdfSynchronizer
const ScalarFunctionSynchronizer< P_V, P_M > * m_targetPdfSynchronizer
Definition: MetropolisHastingsSG.h:262
QUESO::VectorSet::vectorSpace
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
QUESO::BaseVectorRV::env
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:72
QUESO::BaseEnvironment::optionsInputFileName
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:341
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::MetropolisHastingsSG::m_targetPdf
const BaseJointPdf< P_V, P_M > & m_targetPdf
Definition: MetropolisHastingsSG.h:256
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::BaseVectorRV::imageSet
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:79
QUESO::MetropolisHastingsSG::m_alternativeOptionsValues
MhOptionsValues m_alternativeOptionsValues
Definition: MetropolisHastingsSG.h:277
QUESO::MetropolisHastingsSG::commonConstructor
void commonConstructor()
Reads the options values from the options input file.
Definition: MetropolisHastingsSG.C:299
QUESO::MetropolisHastingsSG::m_idsOfUniquePositions
std::vector< unsigned int > m_idsOfUniquePositions
Definition: MetropolisHastingsSG.h:267
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::MetropolisHastingsSG::m_parameterEnabledStatus
std::vector< bool > m_parameterEnabledStatus
Definition: MetropolisHastingsSG.h:261
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::MetropolisHastingsSG::m_alphaQuotients
std::vector< double > m_alphaQuotients
Definition: MetropolisHastingsSG.h:269
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 211 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_optionsObj, QUESO::MetropolisHastingsSGOptions::m_ov, QUESO::MhOptionsValues::m_totallyMute, and QUESO::BaseEnvironment::subDisplayFile().

216  :
217  m_env (sourceRv.env()),
218  m_vectorSpace (sourceRv.imageSet().vectorSpace()),
219  m_targetPdf (sourceRv.pdf()),
220  m_initialPosition (initialPosition),
221  m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
222  m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
223  m_numDisabledParameters (0), // gpmsa2
224  m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true), // gpmsa2
225  m_targetPdfSynchronizer (new ScalarFunctionSynchronizer<P_V,P_M>(m_targetPdf,m_initialPosition)),
226  m_tk (NULL),
227  m_positionIdForDebugging (0),
228  m_stageIdForDebugging (0),
229  m_idsOfUniquePositions (0),//0.),
230  m_logTargets (0),//0.),
231  m_alphaQuotients (0),//0.),
232  m_lastChainSize (0),
233  m_lastMean (NULL),
234  m_lastAdaptedCovMatrix (NULL),
235 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
236  m_alternativeOptionsValues (NULL,NULL),
237 #else
238  m_alternativeOptionsValues (),
239 #endif
240  m_optionsObj (new MetropolisHastingsSGOptions(mlOptions))
241 {
242  if (inputProposalCovMatrix != NULL) {
243  m_initialProposalCovMatrix = *inputProposalCovMatrix;
244  if ((m_env.subDisplayFile() ) &&
245  (m_optionsObj->m_ov.m_totallyMute == false)) {
246  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::constructor(2)"
247  << ": just set m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
248  << std::endl;
249  }
250  }
251  if ((m_env.subDisplayFile() ) &&
252  (m_optionsObj->m_ov.m_totallyMute == false)) {
253  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::constructor(2)"
254  << std::endl;
255  }
256 
257  commonConstructor();
258 
259  if ((m_env.subDisplayFile() ) &&
260  (m_optionsObj->m_ov.m_totallyMute == false)) {
261  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::constructor(2)"
262  << std::endl;
263  }
264 }
QUESO::MetropolisHastingsSG::m_logTargets
std::vector< double > m_logTargets
Definition: MetropolisHastingsSG.h:268
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::BaseVectorRV::pdf
const BaseJointPdf< V, M > & pdf() const
Posterior Density Function of the vector RV; access to private attribute m_pdf.
Definition: VectorRV.C:86
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:218
QUESO::MetropolisHastingsSG::m_targetPdfSynchronizer
const ScalarFunctionSynchronizer< P_V, P_M > * m_targetPdfSynchronizer
Definition: MetropolisHastingsSG.h:262
QUESO::VectorSet::vectorSpace
virtual const VectorSpace< V, M > & vectorSpace() const =0
Vector space to which this set belongs to. See template specialization.
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
QUESO::BaseVectorRV::env
const BaseEnvironment & env() const
QUESO environment; access to private attribute m_env.
Definition: VectorRV.C:72
QUESO::MetropolisHastingsSG::m_targetPdf
const BaseJointPdf< P_V, P_M > & m_targetPdf
Definition: MetropolisHastingsSG.h:256
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::BaseVectorRV::imageSet
const VectorSet< V, M > & imageSet() const
Image set of the vector RV; access to private attribute m_imageSet.
Definition: VectorRV.C:79
QUESO::MetropolisHastingsSG::m_alternativeOptionsValues
MhOptionsValues m_alternativeOptionsValues
Definition: MetropolisHastingsSG.h:277
QUESO::MetropolisHastingsSG::commonConstructor
void commonConstructor()
Reads the options values from the options input file.
Definition: MetropolisHastingsSG.C:299
QUESO::MetropolisHastingsSG::m_idsOfUniquePositions
std::vector< unsigned int > m_idsOfUniquePositions
Definition: MetropolisHastingsSG.h:267
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::MetropolisHastingsSG::m_parameterEnabledStatus
std::vector< bool > m_parameterEnabledStatus
Definition: MetropolisHastingsSG.h:261
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::MetropolisHastingsSG::m_alphaQuotients
std::vector< double > m_alphaQuotients
Definition: MetropolisHastingsSG.h:269
template<class P_V , class P_M >
QUESO::MetropolisHastingsSG< P_V, P_M >::~MetropolisHastingsSG ( )

Destructor.

Definition at line 267 of file MetropolisHastingsSG.C.

268 {
269  //if (m_env.subDisplayFile()) {
270  // *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::destructor()"
271  // << std::endl;
272  //}
273 
274  if (m_lastAdaptedCovMatrix) delete m_lastAdaptedCovMatrix;
275  if (m_lastMean) delete m_lastMean;
276  m_lastChainSize = 0;
277  m_rawChainInfo.reset();
278  m_alphaQuotients.clear();
279  m_logTargets.clear();
280  m_numDisabledParameters = 0; // gpmsa2
281  m_parameterEnabledStatus.clear(); // gpmsa2
282  m_positionIdForDebugging = 0;
283  m_stageIdForDebugging = 0;
284  m_idsOfUniquePositions.clear();
285 
286  if (m_tk ) delete m_tk;
287  if (m_targetPdfSynchronizer) delete m_targetPdfSynchronizer;
288  if (m_optionsObj ) delete m_optionsObj;
289 
290  //if (m_env.subDisplayFile()) {
291  // *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::destructor()"
292  // << std::endl;
293  //}
294 }
QUESO::MetropolisHastingsSG::m_logTargets
std::vector< double > m_logTargets
Definition: MetropolisHastingsSG.h:268
QUESO::MetropolisHastingsSG::m_targetPdfSynchronizer
const ScalarFunctionSynchronizer< P_V, P_M > * m_targetPdfSynchronizer
Definition: MetropolisHastingsSG.h:262
QUESO::MHRawChainInfoStruct::reset
void reset()
Resets Metropolis-Hastings chain info.
Definition: MetropolisHastingsSG.C:77
QUESO::MetropolisHastingsSG::m_rawChainInfo
MHRawChainInfoStruct m_rawChainInfo
Definition: MetropolisHastingsSG.h:275
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::MetropolisHastingsSG::m_idsOfUniquePositions
std::vector< unsigned int > m_idsOfUniquePositions
Definition: MetropolisHastingsSG.h:267
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::MetropolisHastingsSG::m_parameterEnabledStatus
std::vector< bool > m_parameterEnabledStatus
Definition: MetropolisHastingsSG.h:261
QUESO::MetropolisHastingsSG::m_alphaQuotients
std::vector< double > m_alphaQuotients
Definition: MetropolisHastingsSG.h:269

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 712 of file MetropolisHastingsSG.C.

713 {
714  bool result = false;
715 
716  if (alpha <= 0. ) result = false;
717  else if (alpha >= 1. ) result = true;
718  else if (alpha >= m_env.rngObject()->uniformSample()) result = true;
719  else result = false;
720 
721  return result;
722 }
QUESO::BaseEnvironment::rngObject
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:466
QUESO::RngBase::uniformSample
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
QUESO::MetropolisHastingsSG::alpha
double alpha(const MarkovChainPositionData< P_V > &x, const MarkovChainPositionData< P_V > &y, unsigned int xStageId, unsigned int yStageId, double *alphaQuotientPtr=NULL)
Calculates acceptance ration.
Definition: MetropolisHastingsSG.C:393
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
template<class P_V , class P_M >
double QUESO::MetropolisHastingsSG< P_V, P_M >::alpha ( const MarkovChainPositionData< P_V > &  x,
const MarkovChainPositionData< P_V > &  y,
unsigned int  xStageId,
unsigned int  yStageId,
double *  alphaQuotientPtr = NULL 
)
private

Calculates acceptance ration.

It is called by alpha(const std::vector<MarkovChainPositionData<P_V>*>& inputPositions, const std::vector<unsigned int>& inputTKStageIds);

Definition at line 393 of file MetropolisHastingsSG.C.

References QUESO::GaussianJointPdf< V, M >::lawCovMatrix(), QUESO::GaussianJointPdf< V, M >::lawExpVector(), QUESO::GaussianJointPdf< V, M >::lawVarVector(), QUESO::MarkovChainPositionData< V >::logTarget(), QUESO::MarkovChainPositionData< V >::outOfTargetSupport(), and QUESO::MarkovChainPositionData< V >::vecValues().

399 {
400  double alphaQuotient = 0.;
401  if ((x.outOfTargetSupport() == false) &&
402  (y.outOfTargetSupport() == false)) {
403  if ((x.logTarget() == -INFINITY) ||
404  (x.logTarget() == INFINITY) ||
405  ( (boost::math::isnan)(x.logTarget()) )) {
406  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
407  << ", worldRank " << m_env.worldRank()
408  << ", fullRank " << m_env.fullRank()
409  << ", subEnvironment " << m_env.subId()
410  << ", subRank " << m_env.subRank()
411  << ", inter0Rank " << m_env.inter0Rank()
412  << ", positionId = " << m_positionIdForDebugging
413  << ", stageId = " << m_stageIdForDebugging
414  << ": x.logTarget() = " << x.logTarget()
415  << ", x.values() = " << x.vecValues()
416  << ", y.values() = " << y.vecValues()
417  << std::endl;
418  }
419  else if ((y.logTarget() == -INFINITY ) ||
420  (y.logTarget() == INFINITY ) ||
421  ( (boost::math::isnan)(y.logTarget()) )) {
422  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
423  << ", worldRank " << m_env.worldRank()
424  << ", fullRank " << m_env.fullRank()
425  << ", subEnvironment " << m_env.subId()
426  << ", subRank " << m_env.subRank()
427  << ", inter0Rank " << m_env.inter0Rank()
428  << ", positionId = " << m_positionIdForDebugging
429  << ", stageId = " << m_stageIdForDebugging
430  << ": y.logTarget() = " << y.logTarget()
431  << ", x.values() = " << x.vecValues()
432  << ", y.values() = " << y.vecValues()
433  << std::endl;
434  }
435  else {
436  double yLogTargetToUse = y.logTarget();
437  if (m_tk->symmetric()) {
438  alphaQuotient = std::exp(yLogTargetToUse - x.logTarget());
439  if ((m_env.subDisplayFile() ) &&
440  (m_env.displayVerbosity() >= 3 ) &&
441  (m_optionsObj->m_ov.m_totallyMute == false)) {
442  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
443  << ": symmetric proposal case"
444  << ", x = " << x.vecValues()
445  << ", y = " << y.vecValues()
446  << ", yLogTargetToUse = " << yLogTargetToUse
447  << ", x.logTarget() = " << x.logTarget()
448  << ", alpha = " << alphaQuotient
449  << std::endl;
450  }
451  }
452  else {
453 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
454  double qyx = m_tk->rv(yStageId).pdf().lnValue(x.vecValues(),NULL,NULL,NULL,NULL);
455 #else
456  double qyx = -.5 * m_tk->rv(yStageId).pdf().lnValue(x.vecValues(),NULL,NULL,NULL,NULL);
457 #endif
458  if ((m_env.subDisplayFile() ) &&
459  (m_env.displayVerbosity() >= 10 ) &&
460  (m_optionsObj->m_ov.m_totallyMute == false)) {
461  const GaussianJointPdf<P_V,P_M>* pdfYX = dynamic_cast< const GaussianJointPdf<P_V,P_M>* >(&(m_tk->rv(yStageId).pdf()));
462  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
463  << ", rvYX.lawExpVector = " << pdfYX->lawExpVector()
464  << ", rvYX.lawVarVector = " << pdfYX->lawVarVector()
465  << ", rvYX.lawCovMatrix = " << pdfYX->lawCovMatrix()
466  << std::endl;
467  }
468 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
469  double qxy = m_tk->rv(xStageId).pdf().lnValue(y.vecValues(),NULL,NULL,NULL,NULL);
470 #else
471  double qxy = -.5 * m_tk->rv(xStageId).pdf().lnValue(y.vecValues(),NULL,NULL,NULL,NULL);
472 #endif
473  if ((m_env.subDisplayFile() ) &&
474  (m_env.displayVerbosity() >= 10 ) &&
475  (m_optionsObj->m_ov.m_totallyMute == false)) {
476  const GaussianJointPdf<P_V,P_M>* pdfXY = dynamic_cast< const GaussianJointPdf<P_V,P_M>* >(&(m_tk->rv(xStageId).pdf()));
477  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
478  << ", rvXY.lawExpVector = " << pdfXY->lawExpVector()
479  << ", rvXY.lawVarVector = " << pdfXY->lawVarVector()
480  << ", rvXY.lawCovMatrix = " << pdfXY->lawCovMatrix()
481  << std::endl;
482  }
483  alphaQuotient = std::exp(yLogTargetToUse +
484  qyx -
485  x.logTarget() -
486  qxy);
487  if ((m_env.subDisplayFile() ) &&
488  (m_env.displayVerbosity() >= 3 ) &&
489  (m_optionsObj->m_ov.m_totallyMute == false)) {
490  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
491  << ": asymmetric proposal case"
492  << ", xStageId = " << xStageId
493  << ", yStageId = " << yStageId
494  << ", x = " << x.vecValues()
495  << ", y = " << y.vecValues()
496  << ", yLogTargetToUse = " << yLogTargetToUse
497  << ", q(y,x) = " << qyx
498  << ", x.logTarget() = " << x.logTarget()
499  << ", q(x,y) = " << qxy
500  << ", alpha = " << alphaQuotient
501  << std::endl;
502  }
503  }
504  } // protection logic against logTarget values
505  }
506  else {
507  if ((m_env.subDisplayFile() ) &&
508  (m_env.displayVerbosity() >= 10 ) &&
509  (m_optionsObj->m_ov.m_totallyMute == false)) {
510  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(x,y)"
511  << ": x.outOfTargetSupport = " << x.outOfTargetSupport()
512  << ", y.outOfTargetSupport = " << y.outOfTargetSupport()
513  << std::endl;
514  }
515  }
516  if (alphaQuotientPtr != NULL) *alphaQuotientPtr = alphaQuotient;
517 
518  return std::min(1.,alphaQuotient);
519 }
QUESO::BaseEnvironment::subRank
int subRank() const
Access function for sub-rank.
Definition: Environment.C:263
QUESO::BaseEnvironment::inter0Rank
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:289
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::BaseEnvironment::fullRank
int fullRank() const
Returns the process full rank.
Definition: Environment.C:241
QUESO::BaseTKGroup::rv
virtual const GaussianVectorRV< V, M > & rv(unsigned int stageId)=0
Gaussian increment property to construct a transition kernel. See template specialization.
QUESO::BaseTKGroup::symmetric
virtual bool symmetric() const =0
Whether or not the matrix is symmetric. See template specialization.
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:436
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:329
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 ration.

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

Definition at line 523 of file MetropolisHastingsSG.C.

References UQ_FATAL_TEST_MACRO.

526 {
527  unsigned int inputSize = inputPositionsData.size();
528  if ((m_env.subDisplayFile() ) &&
529  (m_env.displayVerbosity() >= 10 ) &&
530  (m_optionsObj->m_ov.m_totallyMute == false)) {
531  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
532  << ", inputSize = " << inputSize
533  << std::endl;
534  }
535  UQ_FATAL_TEST_MACRO((inputSize < 2),
536  m_env.worldRank(),
537  "MetropolisHastingsSG<P_V,P_M>::alpha(vec)",
538  "inputPositionsData has size < 2");
539 
540  // If necessary, return 0. right away
541  if (inputPositionsData[0 ]->outOfTargetSupport()) return 0.;
542  if (inputPositionsData[inputSize-1]->outOfTargetSupport()) return 0.;
543 
544  if ((inputPositionsData[0]->logTarget() == -INFINITY ) ||
545  (inputPositionsData[0]->logTarget() == INFINITY ) ||
546  ( (boost::math::isnan)(inputPositionsData[0]->logTarget()) )) {
547  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
548  << ", worldRank " << m_env.worldRank()
549  << ", fullRank " << m_env.fullRank()
550  << ", subEnvironment " << m_env.subId()
551  << ", subRank " << m_env.subRank()
552  << ", inter0Rank " << m_env.inter0Rank()
553  << ", positionId = " << m_positionIdForDebugging
554  << ", stageId = " << m_stageIdForDebugging
555  << ": inputSize = " << inputSize
556  << ", inputPositionsData[0]->logTarget() = " << inputPositionsData[0]->logTarget()
557  << ", [0]->values() = " << inputPositionsData[0]->vecValues()
558  << ", [inputSize - 1]->values() = " << inputPositionsData[inputSize-1]->vecValues()
559  << std::endl;
560  return 0.;
561  }
562  else if ((inputPositionsData[inputSize - 1]->logTarget() == -INFINITY ) ||
563  (inputPositionsData[inputSize - 1]->logTarget() == INFINITY ) ||
564  ( (boost::math::isnan)(inputPositionsData[inputSize - 1]->logTarget()) )) {
565  std::cerr << "WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
566  << ", worldRank " << m_env.worldRank()
567  << ", fullRank " << m_env.fullRank()
568  << ", subEnvironment " << m_env.subId()
569  << ", subRank " << m_env.subRank()
570  << ", inter0Rank " << m_env.inter0Rank()
571  << ", positionId = " << m_positionIdForDebugging
572  << ", stageId = " << m_stageIdForDebugging
573  << ": inputSize = " << inputSize
574  << ", inputPositionsData[inputSize - 1]->logTarget() = " << inputPositionsData[inputSize-1]->logTarget()
575  << ", [0]->values() = " << inputPositionsData[0]->vecValues()
576  << ", [inputSize - 1]->values() = " << inputPositionsData[inputSize-1]->vecValues()
577  << std::endl;
578  return 0.;
579  }
580 
581  // If inputSize is 2, recursion is not needed
582  if (inputSize == 2) return this->alpha(*(inputPositionsData[0 ]),
583  *(inputPositionsData[inputSize - 1]),
584  inputTKStageIds[0],
585  inputTKStageIds[inputSize-1]);
586 
587  // Prepare two vectors of positions
588  std::vector<MarkovChainPositionData<P_V>*> positionsData (inputSize,NULL);
589  std::vector<MarkovChainPositionData<P_V>*> backwardPositionsData (inputSize,NULL);
590 
591  std::vector<unsigned int > tkStageIds (inputSize,0);
592  std::vector<unsigned int > backwardTKStageIds (inputSize,0);
593 
594  std::vector<unsigned int > tkStageIdsLess1(inputSize,0);
595  std::vector<unsigned int > backwardTKStageIdsLess1(inputSize,0);
596 
597  for (unsigned int i = 0; i < inputSize; ++i) {
598  positionsData [i] = inputPositionsData[i];
599  backwardPositionsData [i] = inputPositionsData[inputSize-i-1];
600 
601  tkStageIds [i] = inputTKStageIds [i];
602  backwardTKStageIds [i] = inputTKStageIds [inputSize-i-1];
603 
604  tkStageIdsLess1[i] = inputTKStageIds [i];
605  backwardTKStageIdsLess1[i] = inputTKStageIds [inputSize-i-1];
606  }
607 
608  tkStageIdsLess1.pop_back();
609  backwardTKStageIdsLess1.pop_back();
610 
611  // Initialize cumulative variables
612  double logNumerator = 0.;
613  double logDenominator = 0.;
614  double alphasNumerator = 1.;
615  double alphasDenominator = 1.;
616 
617  // Compute cumulative variables
618  const P_V& _lastTKPosition = m_tk->preComputingPosition( tkStageIds[inputSize-1]);
619  const P_V& _lastBackwardTKPosition = m_tk->preComputingPosition(backwardTKStageIds[inputSize-1]);
620 
621 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
622  double numContrib = m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(_lastBackwardTKPosition,NULL,NULL,NULL,NULL);
623  double denContrib = m_tk->rv( tkStageIdsLess1).pdf().lnValue(_lastTKPosition ,NULL,NULL,NULL,NULL);
624 #else
625  double numContrib = -.5 * m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(_lastBackwardTKPosition,NULL,NULL,NULL,NULL);
626  double denContrib = -.5 * m_tk->rv( tkStageIdsLess1).pdf().lnValue(_lastTKPosition ,NULL,NULL,NULL,NULL);
627 #endif
628  if ((m_env.subDisplayFile() ) &&
629  (m_env.displayVerbosity() >= 10 ) &&
630  (m_optionsObj->m_ov.m_totallyMute == false)) {
631  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
632  << ", inputSize = " << inputSize
633  << ", before loop"
634  << ": numContrib = " << numContrib
635  << ", denContrib = " << denContrib
636  << std::endl;
637  }
638  logNumerator += numContrib;
639  logDenominator += denContrib;
640 
641  for (unsigned int i = 0; i < (inputSize-2); ++i) { // That is why size must be >= 2
642  positionsData.pop_back();
643  backwardPositionsData.pop_back();
644 
645  const P_V& lastTKPosition = m_tk->preComputingPosition( tkStageIds[inputSize-2-i]);
646  const P_V& lastBackwardTKPosition = m_tk->preComputingPosition(backwardTKStageIds[inputSize-2-i]);
647 
648  tkStageIds.pop_back();
649  backwardTKStageIds.pop_back();
650 
651  tkStageIdsLess1.pop_back();
652  backwardTKStageIdsLess1.pop_back();
653 
654 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
655  numContrib = m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(lastBackwardTKPosition,NULL,NULL,NULL,NULL);
656  denContrib = m_tk->rv( tkStageIdsLess1).pdf().lnValue(lastTKPosition ,NULL,NULL,NULL,NULL);
657 #else
658  numContrib = -.5 * m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(lastBackwardTKPosition,NULL,NULL,NULL,NULL);
659  denContrib = -.5 * m_tk->rv( tkStageIdsLess1).pdf().lnValue(lastTKPosition ,NULL,NULL,NULL,NULL);
660 #endif
661  if ((m_env.subDisplayFile() ) &&
662  (m_env.displayVerbosity() >= 10 ) &&
663  (m_optionsObj->m_ov.m_totallyMute == false)) {
664  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
665  << ", inputSize = " << inputSize
666  << ", in loop, i = " << i
667  << ": numContrib = " << numContrib
668  << ", denContrib = " << denContrib
669  << std::endl;
670  }
671  logNumerator += numContrib;
672  logDenominator += denContrib;
673 
674  alphasNumerator *= (1. - this->alpha(backwardPositionsData,backwardTKStageIds));
675  alphasDenominator *= (1. - this->alpha( positionsData, tkStageIds));
676  }
677 
678  double numeratorLogTargetToUse = backwardPositionsData[0]->logTarget();
679  numContrib = numeratorLogTargetToUse;
680  denContrib = positionsData[0]->logTarget();
681  if ((m_env.subDisplayFile() ) &&
682  (m_env.displayVerbosity() >= 10 ) &&
683  (m_optionsObj->m_ov.m_totallyMute == false)) {
684  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
685  << ", inputSize = " << inputSize
686  << ", after loop"
687  << ": numContrib = " << numContrib
688  << ", denContrib = " << denContrib
689  << std::endl;
690  }
691  logNumerator += numContrib;
692  logDenominator += denContrib;
693 
694  if ((m_env.subDisplayFile() ) &&
695  (m_env.displayVerbosity() >= 10 ) &&
696  (m_optionsObj->m_ov.m_totallyMute == false)) {
697  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
698  << ", inputSize = " << inputSize
699  << ": alphasNumerator = " << alphasNumerator
700  << ", alphasDenominator = " << alphasDenominator
701  << ", logNumerator = " << logNumerator
702  << ", logDenominator = " << logDenominator
703  << std::endl;
704  }
705 
706  // Return result
707  return std::min(1.,(alphasNumerator/alphasDenominator)*std::exp(logNumerator-logDenominator));
708 }
QUESO::BaseEnvironment::subRank
int subRank() const
Access function for sub-rank.
Definition: Environment.C:263
QUESO::BaseEnvironment::inter0Rank
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:289
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::BaseEnvironment::fullRank
int fullRank() const
Returns the process full rank.
Definition: Environment.C:241
QUESO::BaseTKGroup::rv
virtual const GaussianVectorRV< V, M > & rv(unsigned int stageId)=0
Gaussian increment property to construct a transition kernel. See template specialization.
QUESO::MetropolisHastingsSG::alpha
double alpha(const MarkovChainPositionData< P_V > &x, const MarkovChainPositionData< P_V > &y, unsigned int xStageId, unsigned int yStageId, double *alphaQuotientPtr=NULL)
Calculates acceptance ration.
Definition: MetropolisHastingsSG.C:393
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:436
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:329
QUESO::BaseTKGroup::preComputingPosition
const V & preComputingPosition(unsigned int stageId) const
Pre-computing position; access to protected attribute *m_preComputingPositions[stageId].
Definition: TKGroup.C:85
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 299 of file MetropolisHastingsSG.C.

References UQ_FATAL_TEST_MACRO.

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

300 {
302  // Instantiate the appropriate TK (transition kernel)
304  if ((m_env.subDisplayFile() ) &&
305  (m_optionsObj->m_ov.m_totallyMute == false)) {
306  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
307  << std::endl;
308  }
309 
310  if (m_optionsObj->m_ov.m_initialPositionDataInputFileName != ".") { // palms
311  std::set<unsigned int> tmpSet;
312  tmpSet.insert(m_env.subId());
313  m_initialPosition.subReadContents((m_optionsObj->m_ov.m_initialPositionDataInputFileName+"_sub"+m_env.subIdString()),
314  m_optionsObj->m_ov.m_initialPositionDataInputFileType,
315  tmpSet);
316  if ((m_env.subDisplayFile() ) &&
317  (m_optionsObj->m_ov.m_totallyMute == false)) {
318  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
319  << ": just read initial position contents = " << m_initialPosition
320  << std::endl;
321  }
322  }
323 
324  if (m_optionsObj->m_ov.m_parameterDisabledSet.size() > 0) { // gpmsa2
325  for (std::set<unsigned int>::iterator setIt = m_optionsObj->m_ov.m_parameterDisabledSet.begin(); setIt != m_optionsObj->m_ov.m_parameterDisabledSet.end(); ++setIt) {
326  unsigned int paramId = *setIt;
327  if (paramId < m_vectorSpace.dimLocal()) {
328  m_numDisabledParameters++;
329  m_parameterEnabledStatus[paramId] = false;
330  }
331  }
332  }
333 
334  std::vector<double> drScalesAll(m_optionsObj->m_ov.m_drScalesForExtraStages.size()+1,1.);
335  for (unsigned int i = 1; i < (m_optionsObj->m_ov.m_drScalesForExtraStages.size()+1); ++i) {
336  drScalesAll[i] = m_optionsObj->m_ov.m_drScalesForExtraStages[i-1];
337  }
338  if (m_optionsObj->m_ov.m_tkUseLocalHessian) { // sep2011
339  m_tk = new HessianCovMatricesTKGroup<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
340  m_vectorSpace,
341  drScalesAll,
342  *m_targetPdfSynchronizer);
343  if ((m_env.subDisplayFile() ) &&
344  (m_optionsObj->m_ov.m_totallyMute == false)) {
345  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
346  << ": just instantiated a 'HessianCovMatrices' TK class"
347  << std::endl;
348  }
349  }
350  else {
351  if (m_optionsObj->m_ov.m_initialProposalCovMatrixDataInputFileName != ".") { // palms
352  std::set<unsigned int> tmpSet;
353  tmpSet.insert(m_env.subId());
354  m_initialProposalCovMatrix.subReadContents((m_optionsObj->m_ov.m_initialProposalCovMatrixDataInputFileName+"_sub"+m_env.subIdString()),
355  m_optionsObj->m_ov.m_initialProposalCovMatrixDataInputFileType,
356  tmpSet);
357  if ((m_env.subDisplayFile() ) &&
358  (m_optionsObj->m_ov.m_totallyMute == false)) {
359  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
360  << ": just read initial proposal cov matrix contents = " << m_initialProposalCovMatrix
361  << std::endl;
362  }
363  }
364  else {
365  UQ_FATAL_TEST_MACRO(m_nullInputProposalCovMatrix,
366  m_env.worldRank(),
367  "MetropolisHastingsSG<P_V,P_M>::commonConstructor()",
368  "proposal cov matrix should have been passed by user, since, according to the input algorithm options, local Hessians will not be used in the proposal");
369  }
370 
371  m_tk = new ScaledCovMatrixTKGroup<P_V,P_M>(m_optionsObj->m_prefix.c_str(),
372  m_vectorSpace,
373  drScalesAll,
374  m_initialProposalCovMatrix);
375  if ((m_env.subDisplayFile() ) &&
376  (m_optionsObj->m_ov.m_totallyMute == false)) {
377  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
378  << ": just instantiated a 'ScaledCovMatrix' TK class"
379  << std::endl;
380  }
381  }
382 
383  if ((m_env.subDisplayFile() ) &&
384  (m_optionsObj->m_ov.m_totallyMute == false)) {
385  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
386  << std::endl;
387  }
388  return;
389 }
QUESO::BaseEnvironment::subIdString
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:335
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::MhOptionsValues::m_parameterDisabledSet
std::set< unsigned int > m_parameterDisabledSet
Definition: MetropolisHastingsSGOptions.h:144
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::MhOptionsValues::m_initialProposalCovMatrixDataInputFileName
std::string m_initialProposalCovMatrixDataInputFileName
Definition: MetropolisHastingsSGOptions.h:142
QUESO::MetropolisHastingsSG::m_targetPdfSynchronizer
const ScalarFunctionSynchronizer< P_V, P_M > * m_targetPdfSynchronizer
Definition: MetropolisHastingsSG.h:262
QUESO::MhOptionsValues::m_initialProposalCovMatrixDataInputFileType
std::string m_initialProposalCovMatrixDataInputFileType
Definition: MetropolisHastingsSGOptions.h:143
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::MhOptionsValues::m_drScalesForExtraStages
std::vector< double > m_drScalesForExtraStages
Definition: MetropolisHastingsSGOptions.h:176
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::MetropolisHastingsSG::m_parameterEnabledStatus
std::vector< bool > m_parameterEnabledStatus
Definition: MetropolisHastingsSG.h:261
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:329
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 1270 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::estimateConvBrooksGelman(), QUESO::BaseVectorSequence< V, M >::getPositionValues(), QUESO::MarkovChainPositionData< V >::logLikelihood(), QUESO::MarkovChainPositionData< V >::logTarget(), QUESO::MiscGetEllapsedSeconds(), QUESO::BaseVectorSequence< V, M >::name(), QUESO::SequenceOfVectors< V, M >::resizeSequence(), QUESO::ScalarSequence< T >::resizeSequence(), QUESO::BaseVectorSequence< V, M >::resizeSequence(), QUESO::SequenceOfVectors< V, M >::setPositionValues(), QUESO::BaseVectorSequence< V, M >::setPositionValues(), QUESO::BaseVectorSequence< V, M >::subSequenceSize(), QUESO::SequenceOfVectors< V, M >::subSequenceSize(), QUESO::BaseVectorSequence< V, M >::subWriteContents(), QUESO::ScalarSequence< T >::subWriteContents(), QUESO::ScaledCovMatrixTKGroup< V, M >::updateLawCovMatrix(), UQ_FATAL_RC_MACRO, UQ_FATAL_TEST_MACRO, QUESO::UQ_INCOMPLETE_IMPLEMENTATION_RC, QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC, QUESO::UQ_OK_RC, and QUESO::MarkovChainPositionData< V >::vecValues().

1276 {
1277  //m_env.syncPrintDebugMsg("Entering MetropolisHastingsSG<P_V,P_M>::generateFullChain()",3,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
1278 
1279  if ((m_env.subDisplayFile() ) &&
1280  (m_optionsObj->m_ov.m_totallyMute == false)) {
1281  *m_env.subDisplayFile() << "Starting the generation of Markov chain " << workingChain.name()
1282  << ", with " << chainSize
1283  << " positions..."
1284  << std::endl;
1285  }
1286 
1287  int iRC = UQ_OK_RC;
1288  struct timeval timevalChain;
1289  struct timeval timevalCandidate;
1290  struct timeval timevalTarget;
1291  struct timeval timevalMhAlpha;
1292  struct timeval timevalDrAlpha;
1293  struct timeval timevalDR;
1294  struct timeval timevalAM;
1295 
1296  m_positionIdForDebugging = 0;
1297  m_stageIdForDebugging = 0;
1298 
1299  m_rawChainInfo.reset();
1300 
1301  iRC = gettimeofday(&timevalChain, NULL);
1302 
1303  if ((m_env.subDisplayFile() ) &&
1304  (m_optionsObj->m_ov.m_totallyMute == false)) {
1305  *m_env.subDisplayFile() << "\nIn MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1306  << ": contents of initial position are:";
1307  *m_env.subDisplayFile() << valuesOf1stPosition; // FIX ME: might need parallelism
1308  *m_env.subDisplayFile() << "\nIn MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1309  << ": targetPdf.domaintSet() info is:"
1310  << m_targetPdf.domainSet();
1311  *m_env.subDisplayFile() << std::endl;
1312  }
1313 
1314  bool outOfTargetSupport = !m_targetPdf.domainSet().contains(valuesOf1stPosition);
1315  if ((m_env.subDisplayFile()) &&
1316  (outOfTargetSupport )) {
1317  *m_env.subDisplayFile() << "ERROR: In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1318  << ": contents of initial position are:\n";
1319  *m_env.subDisplayFile() << valuesOf1stPosition; // FIX ME: might need parallelism
1320  *m_env.subDisplayFile() << "\nERROR: In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1321  << ": targetPdf.domaintSet() info is:\n"
1322  << m_targetPdf.domainSet();
1323  *m_env.subDisplayFile() << std::endl;
1324  }
1325  UQ_FATAL_TEST_MACRO(outOfTargetSupport,
1326  m_env.worldRank(),
1327  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
1328  "initial position should not be out of target pdf support");
1329  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalTarget, NULL);
1330  double logPrior = 0.;
1331  double logLikelihood = 0.;
1332 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
1333  double logTarget = m_targetPdfSynchronizer->callFunction(&valuesOf1stPosition,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment // KEY
1334 #else
1335  double logTarget = -0.5 * m_targetPdfSynchronizer->callFunction(&valuesOf1stPosition,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1336 #endif
1337  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime += MiscGetEllapsedSeconds(&timevalTarget);
1338  m_rawChainInfo.numTargetCalls++;
1339  if ((m_env.subDisplayFile() ) &&
1340  (m_env.displayVerbosity() >= 3 ) &&
1341  (m_optionsObj->m_ov.m_totallyMute == false)) {
1342  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1343  << ": just returned from likelihood() for initial chain position"
1344  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1345  << ", logPrior = " << logPrior
1346  << ", logLikelihood = " << logLikelihood
1347  << ", logTarget = " << logTarget
1348  << std::endl;
1349  }
1350 
1351  //*m_env.subDisplayFile() << "AQUI 001" << std::endl;
1352  MarkovChainPositionData<P_V> currentPositionData(m_env,
1353  valuesOf1stPosition,
1354  outOfTargetSupport,
1355  logLikelihood,
1356  logTarget);
1357 
1358  P_V gaussianVector(m_vectorSpace.zeroVector());
1359  P_V tmpVecValues(m_vectorSpace.zeroVector());
1360  MarkovChainPositionData<P_V> currentCandidateData(m_env);
1361 
1362  //****************************************************
1363  // Set chain position with positionId = 0
1364  //****************************************************
1365  workingChain.resizeSequence(chainSize);
1366  m_numPositionsNotSubWritten = 0;
1367  if (workingLogLikelihoodValues) workingLogLikelihoodValues->resizeSequence(chainSize);
1368  if (workingLogTargetValues ) workingLogTargetValues->resizeSequence (chainSize);
1369  if (true/*m_uniqueChainGenerate*/) m_idsOfUniquePositions.resize(chainSize,0);
1370  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1371  m_logTargets.resize (chainSize,0.);
1372  m_alphaQuotients.resize(chainSize,0.);
1373  }
1374 
1375  unsigned int uniquePos = 0;
1376  workingChain.setPositionValues(0,currentPositionData.vecValues());
1377  m_numPositionsNotSubWritten++;
1378  if ((m_optionsObj->m_ov.m_rawChainDataOutputPeriod > 0 ) &&
1379  (((0+1) % m_optionsObj->m_ov.m_rawChainDataOutputPeriod) == 0 ) &&
1380  (m_optionsObj->m_ov.m_rawChainDataOutputFileName != ".")) {
1381  workingChain.subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1382  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1383  m_optionsObj->m_ov.m_rawChainDataOutputFileName,
1384  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1385  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1386  if ((m_env.subDisplayFile() ) &&
1387  (m_optionsObj->m_ov.m_totallyMute == false)) {
1388  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1389  << ": just wrote (per period request) " << m_numPositionsNotSubWritten << " chain positions "
1390  << ", " << 0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod << " <= pos <= " << 0
1391  << std::endl;
1392  }
1393 
1394  if (workingLogLikelihoodValues) {
1395  workingLogLikelihoodValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1396  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1397  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_likelihood",
1398  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1399  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1400  }
1401 
1402  if (workingLogTargetValues) {
1403  workingLogTargetValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1404  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1405  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_target",
1406  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1407  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1408  }
1409 
1410  m_numPositionsNotSubWritten = 0;
1411  }
1412 
1413  if (workingLogLikelihoodValues) (*workingLogLikelihoodValues)[0] = currentPositionData.logLikelihood();
1414  if (workingLogTargetValues ) (*workingLogTargetValues )[0] = currentPositionData.logTarget();
1415  if (true/*m_uniqueChainGenerate*/) m_idsOfUniquePositions[uniquePos++] = 0;
1416  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1417  m_logTargets [0] = currentPositionData.logTarget();
1418  m_alphaQuotients[0] = 1.;
1419  }
1420  //*m_env.subDisplayFile() << "AQUI 002" << std::endl;
1421 
1422  if ((m_env.subDisplayFile() ) &&
1423  (m_env.displayVerbosity() >= 10 ) &&
1424  (m_optionsObj->m_ov.m_totallyMute == false)) {
1425  *m_env.subDisplayFile() << "\n"
1426  << "\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
1427  << "\n"
1428  << std::endl;
1429  }
1430 
1431  //m_env.syncPrintDebugMsg("In MetropolisHastingsSG<P_V,P_M>::generateFullChain(), right before main loop",3,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
1432 
1433  //****************************************************
1434  // Begin chain loop from positionId = 1
1435  //****************************************************
1436  if ((m_env.numSubEnvironments() < (unsigned int) m_env.fullComm().NumProc()) &&
1437  (m_initialPosition.numOfProcsForStorage() == 1 ) &&
1438  (m_env.subRank() != 0 )) {
1439  // subRank != 0 --> Enter the barrier and wait for processor 0 to decide to call the targetPdf
1440  double aux = 0.;
1441  aux = m_targetPdfSynchronizer->callFunction(NULL,
1442  NULL,
1443  NULL,
1444  NULL,
1445  NULL,
1446  NULL,
1447  NULL);
1448  if (aux) {}; // just to remove compiler warning
1449  for (unsigned int positionId = 1; positionId < workingChain.subSequenceSize(); ++positionId) {
1450  // Multiply by position values by 'positionId' in order to avoid a constant sequence,
1451  // which would cause zero variance and eventually OVERFLOW flags raised
1452  workingChain.setPositionValues(positionId,((double) positionId) * currentPositionData.vecValues());
1453  m_rawChainInfo.numRejections++;
1454  }
1455  }
1456  else for (unsigned int positionId = 1; positionId < workingChain.subSequenceSize(); ++positionId) {
1457  //****************************************************
1458  // Point 1/6 of logic for new position
1459  // Loop: initialize variables and print some information
1460  //****************************************************
1461  m_positionIdForDebugging = positionId;
1462  if ((m_env.subDisplayFile() ) &&
1463  (m_env.displayVerbosity() >= 3 ) &&
1464  (m_optionsObj->m_ov.m_totallyMute == false)) {
1465  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1466  << ": beginning chain position of id = " << positionId
1467  << ", m_optionsObj->m_ov.m_drMaxNumExtraStages = " << m_optionsObj->m_ov.m_drMaxNumExtraStages
1468  << std::endl;
1469  }
1470  unsigned int stageId = 0;
1471  m_stageIdForDebugging = stageId;
1472 
1473  m_tk->clearPreComputingPositions();
1474 
1475  if ((m_env.subDisplayFile() ) &&
1476  (m_env.displayVerbosity() >= 5 ) &&
1477  (m_optionsObj->m_ov.m_totallyMute == false)) {
1478  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1479  << ": about to set TK pre computing position of local id " << 0
1480  << ", values = " << currentPositionData.vecValues()
1481  << std::endl;
1482  }
1483  bool validPreComputingPosition = m_tk->setPreComputingPosition(currentPositionData.vecValues(),0);
1484  if ((m_env.subDisplayFile() ) &&
1485  (m_env.displayVerbosity() >= 5 ) &&
1486  (m_optionsObj->m_ov.m_totallyMute == false)) {
1487  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1488  << ": returned from setting TK pre computing position of local id " << 0
1489  << ", values = " << currentPositionData.vecValues()
1490  << ", valid = " << validPreComputingPosition
1491  << std::endl;
1492  }
1493  UQ_FATAL_TEST_MACRO(validPreComputingPosition == false,
1494  m_env.worldRank(),
1495  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
1496  "initial position should not be an invalid pre computing position");
1497 
1498  //****************************************************
1499  // Point 2/6 of logic for new position
1500  // Loop: generate new position
1501  //****************************************************
1502  // sep2011
1503  bool keepGeneratingCandidates = true;
1504  while (keepGeneratingCandidates) {
1505  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalCandidate, NULL);
1506  m_tk->rv(0).realizer().realization(tmpVecValues);
1507  if (m_numDisabledParameters > 0) { // gpmsa2
1508  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
1509  if (m_parameterEnabledStatus[paramId] == false) {
1510  tmpVecValues[paramId] = m_initialPosition[paramId];
1511  }
1512  }
1513  }
1514  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.candidateRunTime += MiscGetEllapsedSeconds(&timevalCandidate);
1515 
1516  outOfTargetSupport = !m_targetPdf.domainSet().contains(tmpVecValues);
1517 
1518  bool displayDetail = (m_env.displayVerbosity() >= 10/*99*/) || m_optionsObj->m_ov.m_displayCandidates;
1519  if ((m_env.subDisplayFile() ) &&
1520  (displayDetail ) &&
1521  (m_optionsObj->m_ov.m_totallyMute == false)) {
1522  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1523  << ": for chain position of id = " << positionId
1524  << ", candidate = " << tmpVecValues // FIX ME: might need parallelism
1525  << ", outOfTargetSupport = " << outOfTargetSupport
1526  << std::endl;
1527  }
1528 
1529  if (m_optionsObj->m_ov.m_putOutOfBoundsInChain) keepGeneratingCandidates = false;
1530  else keepGeneratingCandidates = outOfTargetSupport;
1531  }
1532 
1533  if ((m_env.subDisplayFile() ) &&
1534  (m_env.displayVerbosity() >= 5 ) &&
1535  (m_optionsObj->m_ov.m_totallyMute == false)) {
1536  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1537  << ": about to set TK pre computing position of local id " << stageId+1
1538  << ", values = " << tmpVecValues
1539  << std::endl;
1540  }
1541  validPreComputingPosition = m_tk->setPreComputingPosition(tmpVecValues,stageId+1);
1542  if ((m_env.subDisplayFile() ) &&
1543  (m_env.displayVerbosity() >= 5 ) &&
1544  (m_optionsObj->m_ov.m_totallyMute == false)) {
1545  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1546  << ": returned from setting TK pre computing position of local id " << stageId+1
1547  << ", values = " << tmpVecValues
1548  << ", valid = " << validPreComputingPosition
1549  << std::endl;
1550  }
1551 
1552  if (outOfTargetSupport) {
1553  m_rawChainInfo.numOutOfTargetSupport++;
1554  logPrior = -INFINITY;
1555  logLikelihood = -INFINITY;
1556  logTarget = -INFINITY;
1557  }
1558  else {
1559  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalTarget, NULL);
1560 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
1561  logTarget = m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1562 #else
1563  logTarget = -0.5 * m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1564 #endif
1565  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime += MiscGetEllapsedSeconds(&timevalTarget);
1566  m_rawChainInfo.numTargetCalls++;
1567  if ((m_env.subDisplayFile() ) &&
1568  (m_env.displayVerbosity() >= 3 ) &&
1569  (m_optionsObj->m_ov.m_totallyMute == false)) {
1570  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1571  << ": just returned from likelihood() for chain position of id " << positionId
1572  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1573  << ", logPrior = " << logPrior
1574  << ", logLikelihood = " << logLikelihood
1575  << ", logTarget = " << logTarget
1576  << std::endl;
1577  }
1578  }
1579  currentCandidateData.set(tmpVecValues,
1580  outOfTargetSupport,
1581  logLikelihood,
1582  logTarget);
1583 
1584  if ((m_env.subDisplayFile() ) &&
1585  (m_env.displayVerbosity() >= 10 ) &&
1586  (m_optionsObj->m_ov.m_totallyMute == false)) {
1587  *m_env.subDisplayFile() << "\n"
1588  << "\n-----------------------------------------------------------\n"
1589  << "\n"
1590  << std::endl;
1591  }
1592  bool accept = false;
1593  double alphaFirstCandidate = 0.;
1594  if (outOfTargetSupport) {
1595  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1596  m_alphaQuotients[positionId] = 0.;
1597  }
1598  }
1599  else {
1600  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalMhAlpha, NULL);
1601  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1602  alphaFirstCandidate = this->alpha(currentPositionData,currentCandidateData,0,1,&m_alphaQuotients[positionId]);
1603  }
1604  else {
1605  alphaFirstCandidate = this->alpha(currentPositionData,currentCandidateData,0,1,NULL);
1606  }
1607  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.mhAlphaRunTime += MiscGetEllapsedSeconds(&timevalMhAlpha);
1608  if ((m_env.subDisplayFile() ) &&
1609  (m_env.displayVerbosity() >= 10 ) &&
1610  (m_optionsObj->m_ov.m_totallyMute == false)) {
1611  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1612  << ": for chain position of id = " << positionId
1613  << std::endl;
1614  }
1615  accept = acceptAlpha(alphaFirstCandidate);
1616  }
1617 
1618  bool displayDetail = (m_env.displayVerbosity() >= 10/*99*/) || m_optionsObj->m_ov.m_displayCandidates;
1619  if ((m_env.subDisplayFile() ) &&
1620  (displayDetail ) &&
1621  (m_optionsObj->m_ov.m_totallyMute == false)) {
1622  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1623  << ": for chain position of id = " << positionId
1624  << ", outOfTargetSupport = " << outOfTargetSupport
1625  << ", alpha = " << alphaFirstCandidate
1626  << ", accept = " << accept
1627  << ", currentCandidateData.vecValues() = ";
1628  *m_env.subDisplayFile() << currentCandidateData.vecValues(); // FIX ME: might need parallelism
1629  *m_env.subDisplayFile() << "\n"
1630  << "\n curLogTarget = " << currentPositionData.logTarget()
1631  << "\n"
1632  << "\n canLogTarget = " << currentCandidateData.logTarget()
1633  << "\n"
1634  << std::endl;
1635  }
1636  if ((m_env.subDisplayFile() ) &&
1637  (m_env.displayVerbosity() >= 10 ) &&
1638  (m_optionsObj->m_ov.m_totallyMute == false)) {
1639  *m_env.subDisplayFile() << "\n"
1640  << "\n-----------------------------------------------------------\n"
1641  << "\n"
1642  << std::endl;
1643  }
1644 
1645  //****************************************************
1646  // Point 3/6 of logic for new position
1647  // Loop: delayed rejection
1648  //****************************************************
1649  // sep2011
1650  std::vector<MarkovChainPositionData<P_V>*> drPositionsData(stageId+2,NULL);
1651  std::vector<unsigned int> tkStageIds (stageId+2,0);
1652  if ((accept == false) &&
1653  (outOfTargetSupport == false) && // IMPORTANT
1654  (m_optionsObj->m_ov.m_drMaxNumExtraStages > 0 )) {
1655  if ((m_optionsObj->m_ov.m_drDuringAmNonAdaptiveInt == false ) &&
1656  (m_optionsObj->m_ov.m_tkUseLocalHessian == false ) &&
1657  (m_optionsObj->m_ov.m_amInitialNonAdaptInterval > 0 ) &&
1658  (m_optionsObj->m_ov.m_amAdaptInterval > 0 ) &&
1659  (positionId <= m_optionsObj->m_ov.m_amInitialNonAdaptInterval)) {
1660  // Avoid DR now
1661  }
1662  else {
1663  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalDR, NULL);
1664 
1665  drPositionsData[0] = new MarkovChainPositionData<P_V>(currentPositionData );
1666  drPositionsData[1] = new MarkovChainPositionData<P_V>(currentCandidateData);
1667 
1668  tkStageIds[0] = 0;
1669  tkStageIds[1] = 1;
1670 
1671  while ((validPreComputingPosition == true ) &&
1672  (accept == false ) &&
1673  (stageId < m_optionsObj->m_ov.m_drMaxNumExtraStages)) {
1674  if ((m_env.subDisplayFile() ) &&
1675  (m_env.displayVerbosity() >= 10 ) &&
1676  (m_optionsObj->m_ov.m_totallyMute == false)) {
1677  *m_env.subDisplayFile() << "\n"
1678  << "\n+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"
1679  << "\n"
1680  << std::endl;
1681  }
1682  m_rawChainInfo.numDRs++;
1683  stageId++;
1684  m_stageIdForDebugging = stageId;
1685  if ((m_env.subDisplayFile() ) &&
1686  (m_env.displayVerbosity() >= 10 ) &&
1687  (m_optionsObj->m_ov.m_totallyMute == false)) {
1688  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1689  << ": for chain position of id = " << positionId
1690  << ", beginning stageId = " << stageId
1691  << std::endl;
1692  }
1693 
1694  keepGeneratingCandidates = true;
1695  while (keepGeneratingCandidates) {
1696  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalCandidate, NULL);
1697  m_tk->rv(tkStageIds).realizer().realization(tmpVecValues);
1698  if (m_numDisabledParameters > 0) { // gpmsa2
1699  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
1700  if (m_parameterEnabledStatus[paramId] == false) {
1701  tmpVecValues[paramId] = m_initialPosition[paramId];
1702  }
1703  }
1704  }
1705  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.candidateRunTime += MiscGetEllapsedSeconds(&timevalCandidate);
1706 
1707  outOfTargetSupport = !m_targetPdf.domainSet().contains(tmpVecValues);
1708 
1709  if (m_optionsObj->m_ov.m_putOutOfBoundsInChain) keepGeneratingCandidates = false;
1710  else keepGeneratingCandidates = outOfTargetSupport;
1711  }
1712 
1713  if ((m_env.subDisplayFile() ) &&
1714  (m_env.displayVerbosity() >= 5 ) &&
1715  (m_optionsObj->m_ov.m_totallyMute == false)) {
1716  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1717  << ": about to set TK pre computing position of local id " << stageId+1
1718  << ", values = " << tmpVecValues
1719  << std::endl;
1720  }
1721  validPreComputingPosition = m_tk->setPreComputingPosition(tmpVecValues,stageId+1);
1722  if ((m_env.subDisplayFile() ) &&
1723  (m_env.displayVerbosity() >= 5 ) &&
1724  (m_optionsObj->m_ov.m_totallyMute == false)) {
1725  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1726  << ": returned from setting TK pre computing position of local id " << stageId+1
1727  << ", values = " << tmpVecValues
1728  << ", valid = " << validPreComputingPosition
1729  << std::endl;
1730  }
1731 
1732  if (outOfTargetSupport) {
1733  m_rawChainInfo.numOutOfTargetSupportInDR++; // new 2010/May/12
1734  logPrior = -INFINITY;
1735  logLikelihood = -INFINITY;
1736  logTarget = -INFINITY;
1737  }
1738  else {
1739  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalTarget, NULL);
1740 #ifdef QUESO_EXPECTS_LN_LIKELIHOOD_INSTEAD_OF_MINUS_2_LN
1741  logTarget = m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1742 #else
1743  logTarget = -0.5 * m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood); // Might demand parallel environment
1744 #endif
1745  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime += MiscGetEllapsedSeconds(&timevalTarget);
1746  m_rawChainInfo.numTargetCalls++;
1747  if ((m_env.subDisplayFile() ) &&
1748  (m_env.displayVerbosity() >= 3 ) &&
1749  (m_optionsObj->m_ov.m_totallyMute == false)) {
1750  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1751  << ": just returned from likelihood() for chain position of id " << positionId
1752  << ", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1753  << ", stageId = " << stageId
1754  << ", logPrior = " << logPrior
1755  << ", logLikelihood = " << logLikelihood
1756  << ", logTarget = " << logTarget
1757  << std::endl;
1758  }
1759  }
1760  currentCandidateData.set(tmpVecValues,
1761  outOfTargetSupport,
1762  logLikelihood,
1763  logTarget);
1764 
1765  drPositionsData.push_back(new MarkovChainPositionData<P_V>(currentCandidateData));
1766  tkStageIds.push_back (stageId+1);
1767 
1768  double alphaDR = 0.;
1769  if (outOfTargetSupport == false) {
1770  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalDrAlpha, NULL);
1771  alphaDR = this->alpha(drPositionsData,tkStageIds);
1772  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.drAlphaRunTime += MiscGetEllapsedSeconds(&timevalDrAlpha);
1773  accept = acceptAlpha(alphaDR);
1774  }
1775 
1776  displayDetail = (m_env.displayVerbosity() >= 10/*99*/) || m_optionsObj->m_ov.m_displayCandidates;
1777  if ((m_env.subDisplayFile() ) &&
1778  (displayDetail ) &&
1779  (m_optionsObj->m_ov.m_totallyMute == false)) {
1780  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1781  << ": for chain position of id = " << positionId
1782  << " and stageId = " << stageId
1783  << ", outOfTargetSupport = " << outOfTargetSupport
1784  << ", alpha = " << alphaDR
1785  << ", accept = " << accept
1786  << ", currentCandidateData.vecValues() = ";
1787  *m_env.subDisplayFile() << currentCandidateData.vecValues(); // FIX ME: might need parallelism
1788  *m_env.subDisplayFile() << std::endl;
1789  }
1790  } // while
1791 
1792  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.drRunTime += MiscGetEllapsedSeconds(&timevalDR);
1793  } // if-else "Avoid DR now"
1794  } // end of 'delayed rejection' logic
1795 
1796  for (unsigned int i = 0; i < drPositionsData.size(); ++i) {
1797  if (drPositionsData[i]) delete drPositionsData[i];
1798  }
1799 
1800  //****************************************************
1801  // Point 4/6 of logic for new position
1802  // Loop: update chain
1803  //****************************************************
1804  if (accept) {
1805  workingChain.setPositionValues(positionId,currentCandidateData.vecValues());
1806  if (true/*m_uniqueChainGenerate*/) m_idsOfUniquePositions[uniquePos++] = positionId;
1807  currentPositionData = currentCandidateData;
1808  }
1809  else {
1810  workingChain.setPositionValues(positionId,currentPositionData.vecValues());
1811  m_rawChainInfo.numRejections++;
1812  }
1813  m_numPositionsNotSubWritten++;
1814  if ((m_optionsObj->m_ov.m_rawChainDataOutputPeriod > 0 ) &&
1815  (((positionId+1) % m_optionsObj->m_ov.m_rawChainDataOutputPeriod) == 0 ) &&
1816  (m_optionsObj->m_ov.m_rawChainDataOutputFileName != ".")) {
1817  if ((m_env.subDisplayFile() ) &&
1818  (m_env.displayVerbosity() >= 10 ) &&
1819  (m_optionsObj->m_ov.m_totallyMute == false)) {
1820  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1821  << ", for chain position of id = " << positionId
1822  << ": about to write (per period request) " << m_numPositionsNotSubWritten << " chain positions "
1823  << ", " << positionId + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod << " <= pos <= " << positionId
1824  << std::endl;
1825  }
1826  workingChain.subWriteContents(positionId + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1827  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1828  m_optionsObj->m_ov.m_rawChainDataOutputFileName,
1829  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1830  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1831  if ((m_env.subDisplayFile() ) &&
1832  (m_optionsObj->m_ov.m_totallyMute == false)) {
1833  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1834  << ", for chain position of id = " << positionId
1835  << ": just wrote (per period request) " << m_numPositionsNotSubWritten << " chain positions "
1836  << ", " << positionId + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod << " <= pos <= " << positionId
1837  << std::endl;
1838  }
1839 
1840  if (workingLogLikelihoodValues) {
1841  workingLogLikelihoodValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1842  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1843  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_likelihood",
1844  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1845  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1846  }
1847 
1848  if (workingLogTargetValues) {
1849  workingLogTargetValues->subWriteContents(0 + 1 - m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1850  m_optionsObj->m_ov.m_rawChainDataOutputPeriod,
1851  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_target",
1852  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
1853  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
1854  }
1855 
1856  m_numPositionsNotSubWritten = 0;
1857  }
1858 
1859 
1860  if (workingLogLikelihoodValues) (*workingLogLikelihoodValues)[positionId] = currentPositionData.logLikelihood();
1861  if (workingLogTargetValues ) (*workingLogTargetValues )[positionId] = currentPositionData.logTarget();
1862 
1863  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
1864  m_logTargets[positionId] = currentPositionData.logTarget();
1865  }
1866 
1867  if( m_optionsObj->m_ov.m_enableBrooksGelmanConvMonitor > 0 ) {
1868  if( positionId%m_optionsObj->m_ov.m_enableBrooksGelmanConvMonitor == 0 &&
1869  positionId > m_optionsObj->m_ov.m_BrooksGelmanLag+1 ) { //+1 to help ensure there are at least 2 samples to use
1870 
1871  double conv_est = workingChain.estimateConvBrooksGelman( m_optionsObj->m_ov.m_BrooksGelmanLag,
1872  positionId - m_optionsObj->m_ov.m_BrooksGelmanLag );
1873 
1874  if ( m_env.subDisplayFile() ) {
1875  *m_env.subDisplayFile() << "positionId = " << positionId
1876  << ", conv_est = " << conv_est << std::endl;
1877  (*m_env.subDisplayFile()).flush();
1878  }
1879  }
1880  }
1881 
1882  //****************************************************
1883  // Point 5/6 of logic for new position
1884  // Loop: adaptive Metropolis (adaptation of covariance matrix)
1885  //****************************************************
1886  // sep2011
1887  if ((m_optionsObj->m_ov.m_tkUseLocalHessian == false) && // IMPORTANT
1888  (m_optionsObj->m_ov.m_amInitialNonAdaptInterval > 0) &&
1889  (m_optionsObj->m_ov.m_amAdaptInterval > 0)) {
1890  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) iRC = gettimeofday(&timevalAM, NULL);
1891 
1892  // Now might be the moment to adapt
1893  unsigned int idOfFirstPositionInSubChain = 0;
1894  SequenceOfVectors<P_V,P_M> partialChain(m_vectorSpace,0,m_optionsObj->m_prefix+"partialChain");
1895 
1896  // Check if now is indeed the moment to adapt
1897  bool printAdaptedMatrix = false;
1898  if (positionId < m_optionsObj->m_ov.m_amInitialNonAdaptInterval) {
1899  // Do nothing
1900  }
1901  else if (positionId == m_optionsObj->m_ov.m_amInitialNonAdaptInterval) {
1902  idOfFirstPositionInSubChain = 0;
1903  partialChain.resizeSequence(m_optionsObj->m_ov.m_amInitialNonAdaptInterval+1);
1904  m_lastMean = m_vectorSpace.newVector();
1905  m_lastAdaptedCovMatrix = m_vectorSpace.newMatrix();
1906  printAdaptedMatrix = true;
1907  }
1908  else {
1909  unsigned int interval = positionId - m_optionsObj->m_ov.m_amInitialNonAdaptInterval;
1910  if ((interval % m_optionsObj->m_ov.m_amAdaptInterval) == 0) {
1911  idOfFirstPositionInSubChain = positionId - m_optionsObj->m_ov.m_amAdaptInterval;
1912  partialChain.resizeSequence(m_optionsObj->m_ov.m_amAdaptInterval);
1913 
1914  if (m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputPeriod > 0) {
1915  if ((interval % m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputPeriod) == 0) {
1916  printAdaptedMatrix = true;
1917  }
1918  }
1919  }
1920  }
1921 
1922  // If now is indeed the moment to adapt, then do it!
1923  if (partialChain.subSequenceSize() > 0) {
1924  P_V transporterVec(m_vectorSpace.zeroVector());
1925  for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
1926  workingChain.getPositionValues(idOfFirstPositionInSubChain+i,transporterVec);
1927  partialChain.setPositionValues(i,transporterVec);
1928  }
1929  updateAdaptedCovMatrix(partialChain,
1930  idOfFirstPositionInSubChain,
1931  m_lastChainSize,
1932  *m_lastMean,
1933  *m_lastAdaptedCovMatrix);
1934 
1935  if ((printAdaptedMatrix == true) &&
1936  (m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputFileName != "." )) { // palms
1937  char varNamePrefix[64];
1938  sprintf(varNamePrefix,"mat_am%d",positionId);
1939 
1940  char tmpChar[64];
1941  sprintf(tmpChar,"_am%d",positionId);
1942 
1943  std::set<unsigned int> tmpSet;
1944  tmpSet.insert(m_env.subId());
1945 
1946  m_lastAdaptedCovMatrix->subWriteContents(varNamePrefix,
1947  (m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputFileName+tmpChar),
1948  m_optionsObj->m_ov.m_amAdaptedMatricesDataOutputFileType,
1949  tmpSet);
1950  if ((m_env.subDisplayFile() ) &&
1951  (m_optionsObj->m_ov.m_totallyMute == false)) {
1952  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1953  << ": just wrote last adapted proposal cov matrix contents = " << *m_lastAdaptedCovMatrix
1954  << std::endl;
1955  }
1956  } // if (printAdaptedMatrix && ...)
1957 
1958  bool tmpCholIsPositiveDefinite = false;
1959  P_M tmpChol(*m_lastAdaptedCovMatrix);
1960  P_M attemptedMatrix(tmpChol);
1961  if ((m_env.subDisplayFile() ) &&
1962  (m_env.displayVerbosity() >= 10)) {
1963  //(m_optionsObj->m_ov.m_totallyMute == false)) {
1964  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1965  << ", positionId = " << positionId
1966  << ": 'am' calling first tmpChol.chol()"
1967  << std::endl;
1968  }
1969  iRC = tmpChol.chol();
1970  if (iRC) {
1971  std::cerr << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain(): first chol failed\n";
1972  }
1973  if ((m_env.subDisplayFile() ) &&
1974  (m_env.displayVerbosity() >= 10)) {
1975  //(m_optionsObj->m_ov.m_totallyMute == false)) {
1976  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1977  << ", positionId = " << positionId
1978  << ": 'am' got first tmpChol.chol() with iRC = " << iRC
1979  << std::endl;
1980  if (iRC == 0) {
1981  double diagMult = 1.;
1982  for (unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
1983  diagMult *= tmpChol(j,j);
1984  }
1985  *m_env.subDisplayFile() << "diagMult = " << diagMult
1986  << std::endl;
1987  }
1988  }
1989 #if 0 // tentative logic
1990  if (iRC == 0) {
1991  double diagMult = 1.;
1992  for (unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
1993  diagMult *= tmpChol(j,j);
1994  }
1995  if (diagMult < 1.e-40) {
1996  iRC = UQ_MATRIX_IS_NOT_POS_DEFINITE_RC;
1997  }
1998  }
1999 #endif
2000 
2001  if (iRC) {
2002  UQ_FATAL_TEST_MACRO(iRC != UQ_MATRIX_IS_NOT_POS_DEFINITE_RC,
2003  m_env.worldRank(),
2004  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
2005  "invalid iRC returned from first chol()");
2006  // Matrix is not positive definite
2007  P_M* tmpDiag = m_vectorSpace.newDiagMatrix(m_optionsObj->m_ov.m_amEpsilon);
2008  if (m_numDisabledParameters > 0) { // gpmsa2
2009  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2010  if (m_parameterEnabledStatus[paramId] == false) {
2011  (*tmpDiag)(paramId,paramId) = 0.;
2012  }
2013  }
2014  }
2015  tmpChol = *m_lastAdaptedCovMatrix + *tmpDiag;
2016  attemptedMatrix = tmpChol;
2017  delete tmpDiag;
2018  if ((m_env.subDisplayFile() ) &&
2019  (m_env.displayVerbosity() >= 10)) {
2020  //(m_optionsObj->m_ov.m_totallyMute == false)) {
2021  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2022  << ", positionId = " << positionId
2023  << ": 'am' calling second tmpChol.chol()"
2024  << std::endl;
2025  }
2026  iRC = tmpChol.chol();
2027  if (iRC) {
2028  std::cerr << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain(): second chol failed\n";
2029  }
2030  if ((m_env.subDisplayFile() ) &&
2031  (m_env.displayVerbosity() >= 10)) {
2032  //(m_optionsObj->m_ov.m_totallyMute == false)) {
2033  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2034  << ", positionId = " << positionId
2035  << ": 'am' got second tmpChol.chol() with iRC = " << iRC
2036  << std::endl;
2037  if (iRC == 0) {
2038  double diagMult = 1.;
2039  for (unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2040  diagMult *= tmpChol(j,j);
2041  }
2042  *m_env.subDisplayFile() << "diagMult = " << diagMult
2043  << std::endl;
2044  }
2045  else {
2046  *m_env.subDisplayFile() << "attemptedMatrix = " << attemptedMatrix // FIX ME: might demand parallelism
2047  << std::endl;
2048  }
2049  }
2050  if (iRC) {
2051  UQ_FATAL_TEST_MACRO(iRC != UQ_MATRIX_IS_NOT_POS_DEFINITE_RC,
2052  m_env.worldRank(),
2053  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
2054  "invalid iRC returned from second chol()");
2055  // Do nothing
2056  }
2057  else {
2058  tmpCholIsPositiveDefinite = true;
2059  }
2060  }
2061  else {
2062  tmpCholIsPositiveDefinite = true;
2063  }
2064  if (tmpCholIsPositiveDefinite) {
2065  ScaledCovMatrixTKGroup<P_V,P_M>* tempTK = dynamic_cast<ScaledCovMatrixTKGroup<P_V,P_M>* >(m_tk);
2066  P_M tmpMatrix(m_optionsObj->m_ov.m_amEta*attemptedMatrix);
2067  if (m_numDisabledParameters > 0) { // gpmsa2
2068  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2069  if (m_parameterEnabledStatus[paramId] == false) {
2070  for (unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2071  tmpMatrix(i,paramId) = 0.;
2072  }
2073  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2074  tmpMatrix(paramId,j) = 0.;
2075  }
2076  tmpMatrix(paramId,paramId) = 1.;
2077  }
2078  }
2079  }
2080  tempTK->updateLawCovMatrix(tmpMatrix);
2081 
2082 #ifdef UQ_DRAM_MCG_REQUIRES_INVERTED_COV_MATRICES
2083  UQ_FATAL_RC_MACRO(UQ_INCOMPLETE_IMPLEMENTATION_RC,
2084  m_env.worldRank(),
2085  "MetropolisHastingsSG<P_V,P_M>::generateFullChain()",
2086  "need to code the update of m_upperCholProposalPrecMatrices");
2087 #endif
2088  }
2089 
2090  //for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2091  // if (partialChain[i]) delete partialChain[i];
2092  //}
2093  } // if (partialChain.subSequenceSize() > 0)
2094 
2095  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) m_rawChainInfo.amRunTime += MiscGetEllapsedSeconds(&timevalAM);
2096  } // End of 'adaptive Metropolis' logic
2097 
2098  //****************************************************
2099  // Point 6/6 of logic for new position
2100  // Loop: print some information before going to the next chain position
2101  //****************************************************
2102  if ((m_env.subDisplayFile() ) &&
2103  (m_env.displayVerbosity() >= 3 ) &&
2104  (m_optionsObj->m_ov.m_totallyMute == false)) {
2105  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2106  << ": finishing chain position of id = " << positionId
2107  << ", accept = " << accept
2108  << ", curLogTarget = " << currentPositionData.logTarget()
2109  << ", canLogTarget = " << currentCandidateData.logTarget()
2110  << std::endl;
2111  }
2112 
2113  if ((m_optionsObj->m_ov.m_rawChainDisplayPeriod > 0) &&
2114  (((positionId+1) % m_optionsObj->m_ov.m_rawChainDisplayPeriod) == 0)) {
2115  if ((m_env.subDisplayFile() ) &&
2116  (m_optionsObj->m_ov.m_totallyMute == false)) {
2117  *m_env.subDisplayFile() << "Finished generating " << positionId+1
2118  << " positions" // root
2119  << ", curret rejection percentage = " << (100. * ((double) m_rawChainInfo.numRejections)/((double) (positionId+1)))
2120  << " %"
2121  << std::endl;
2122  }
2123  }
2124 
2125  if ((m_env.subDisplayFile() ) &&
2126  (m_env.displayVerbosity() >= 10 ) &&
2127  (m_optionsObj->m_ov.m_totallyMute == false)) {
2128  *m_env.subDisplayFile() << "\n"
2129  << "\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
2130  << "\n"
2131  << std::endl;
2132  }
2133  } // end chain loop [for (unsigned int positionId = 1; positionId < workingChain.subSequenceSize(); ++positionId) {]
2134 
2135  if ((m_env.numSubEnvironments() < (unsigned int) m_env.fullComm().NumProc()) &&
2136  (m_initialPosition.numOfProcsForStorage() == 1 ) &&
2137  (m_env.subRank() == 0 )) {
2138  // subRank == 0 --> Tell all other processors to exit barrier now that the chain has been fully generated
2139  double aux = 0.;
2140  aux = m_targetPdfSynchronizer->callFunction(NULL,
2141  NULL,
2142  NULL,
2143  NULL,
2144  NULL,
2145  NULL,
2146  NULL);
2147  if (aux) {}; // just to remove compiler warning
2148  }
2149 
2150  //****************************************************
2151  // Print basic information about the chain
2152  //****************************************************
2153  m_rawChainInfo.runTime += MiscGetEllapsedSeconds(&timevalChain);
2154  if ((m_env.subDisplayFile() ) &&
2155  (m_optionsObj->m_ov.m_totallyMute == false)) {
2156  *m_env.subDisplayFile() << "Finished the generation of Markov chain " << workingChain.name()
2157  << ", with sub " << workingChain.subSequenceSize()
2158  << " positions";
2159  *m_env.subDisplayFile() << "\nSome information about this chain:"
2160  << "\n Chain run time = " << m_rawChainInfo.runTime
2161  << " seconds";
2162  if (m_optionsObj->m_ov.m_rawChainMeasureRunTimes) {
2163  *m_env.subDisplayFile() << "\n\n Breaking of the chain run time:\n";
2164  *m_env.subDisplayFile() << "\n Candidate run time = " << m_rawChainInfo.candidateRunTime
2165  << " seconds (" << 100.*m_rawChainInfo.candidateRunTime/m_rawChainInfo.runTime
2166  << "%)";
2167  *m_env.subDisplayFile() << "\n Num target calls = " << m_rawChainInfo.numTargetCalls;
2168  *m_env.subDisplayFile() << "\n Target d. run time = " << m_rawChainInfo.targetRunTime
2169  << " seconds (" << 100.*m_rawChainInfo.targetRunTime/m_rawChainInfo.runTime
2170  << "%)";
2171  *m_env.subDisplayFile() << "\n Avg target run time = " << m_rawChainInfo.targetRunTime/((double) m_rawChainInfo.numTargetCalls)
2172  << " seconds";
2173  *m_env.subDisplayFile() << "\n Mh alpha run time = " << m_rawChainInfo.mhAlphaRunTime
2174  << " seconds (" << 100.*m_rawChainInfo.mhAlphaRunTime/m_rawChainInfo.runTime
2175  << "%)";
2176  *m_env.subDisplayFile() << "\n Dr alpha run time = " << m_rawChainInfo.drAlphaRunTime
2177  << " seconds (" << 100.*m_rawChainInfo.drAlphaRunTime/m_rawChainInfo.runTime
2178  << "%)";
2179  *m_env.subDisplayFile() << "\n---------------------- --------------";
2180  double sumRunTime = m_rawChainInfo.candidateRunTime + m_rawChainInfo.targetRunTime + m_rawChainInfo.mhAlphaRunTime + m_rawChainInfo.drAlphaRunTime;
2181  *m_env.subDisplayFile() << "\n Sum = " << sumRunTime
2182  << " seconds (" << 100.*sumRunTime/m_rawChainInfo.runTime
2183  << "%)";
2184  *m_env.subDisplayFile() << "\n\n Other run times:";
2185  *m_env.subDisplayFile() << "\n DR run time = " << m_rawChainInfo.drRunTime
2186  << " seconds (" << 100.*m_rawChainInfo.drRunTime/m_rawChainInfo.runTime
2187  << "%)";
2188  *m_env.subDisplayFile() << "\n AM run time = " << m_rawChainInfo.amRunTime
2189  << " seconds (" << 100.*m_rawChainInfo.amRunTime/m_rawChainInfo.runTime
2190  << "%)";
2191  }
2192  *m_env.subDisplayFile() << "\n Number of DRs = " << m_rawChainInfo.numDRs << "(num_DRs/chain_size = " << (double) m_rawChainInfo.numDRs/(double) workingChain.subSequenceSize()
2193  << ")";
2194  *m_env.subDisplayFile() << "\n Out of target support in DR = " << m_rawChainInfo.numOutOfTargetSupportInDR;
2195  *m_env.subDisplayFile() << "\n Rejection percentage = " << 100. * (double) m_rawChainInfo.numRejections/(double) workingChain.subSequenceSize()
2196  << " %";
2197  *m_env.subDisplayFile() << "\n Out of target support percentage = " << 100. * (double) m_rawChainInfo.numOutOfTargetSupport/(double) workingChain.subSequenceSize()
2198  << " %";
2199  *m_env.subDisplayFile() << std::endl;
2200  }
2201 
2202  //****************************************************
2203  // Release memory before leaving routine
2204  //****************************************************
2205  //m_env.syncPrintDebugMsg("Leaving MetropolisHastingsSG<P_V,P_M>::generateFullChain()",3,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
2206 
2207  return;
2208 }
QUESO::BaseVectorSequence::name
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
Definition: VectorSequence.C:121
QUESO::BaseEnvironment::subRank
int subRank() const
Access function for sub-rank.
Definition: Environment.C:263
QUESO::MetropolisHastingsSG::m_logTargets
std::vector< double > m_logTargets
Definition: MetropolisHastingsSG.h:268
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::BaseTKGroup::clearPreComputingPositions
virtual void clearPreComputingPositions()
Clears the pre-computing positions m_preComputingPositions[stageId].
Definition: TKGroup.C:121
QUESO::ScalarSequence::subWriteContents
void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Writes the sub-sequence to a file.
Definition: ScalarSequence.C:2742
QUESO::BaseVectorSequence::setPositionValues
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Definition: Defines.h:83
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::BaseVectorSequence::resizeSequence
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:218
QUESO::BaseEnvironment::fullComm
const MpiComm & fullComm() const
Access function for MpiComm full communicator.
Definition: Environment.C:247
QUESO::MetropolisHastingsSG::m_targetPdfSynchronizer
const ScalarFunctionSynchronizer< P_V, P_M > * m_targetPdfSynchronizer
Definition: MetropolisHastingsSG.h:262
QUESO::BaseVectorSequence::subWriteContents
virtual void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const =0
Writes info of the sub-sequence to a file. See template specialization.
QUESO::VectorSpace::newDiagMatrix
M * newDiagMatrix(const V &v) const
Creates a diagonal matrix with the elements and size of vector v.
Definition: VectorSpace.C:237
QUESO::UQ_INCOMPLETE_IMPLEMENTATION_RC
const int UQ_INCOMPLETE_IMPLEMENTATION_RC
Definition: Defines.h:76
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:75
QUESO::MHRawChainInfoStruct::reset
void reset()
Resets Metropolis-Hastings chain info.
Definition: MetropolisHastingsSG.C:77
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:181
QUESO::MhOptionsValues::m_amEta
double m_amEta
Proposal covariance scaling factor, usually 2.4 * 2.4 / d.
Definition: MetropolisHastingsSGOptions.h:196
QUESO::MhOptionsValues::m_rawChainDataOutputAllowedSet
std::set< unsigned int > m_rawChainDataOutputAllowedSet
Definition: MetropolisHastingsSGOptions.h:155
QUESO::MetropolisHastingsSG::m_rawChainInfo
MHRawChainInfoStruct m_rawChainInfo
Definition: MetropolisHastingsSG.h:275
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:287
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
QUESO::BaseTKGroup::setPreComputingPosition
virtual bool setPreComputingPosition(const V &position, unsigned int stageId)
Sets the pre-computing positions m_preComputingPositions[stageId] with a new vector of size position...
Definition: TKGroup.C:102
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::MetropolisHastingsSG::m_targetPdf
const BaseJointPdf< P_V, P_M > & m_targetPdf
Definition: MetropolisHastingsSG.h:256
QUESO::MetropolisHastingsSG::m_tk
BaseTKGroup< P_V, P_M > * m_tk
Definition: MetropolisHastingsSG.h:264
QUESO::ScalarFunctionSynchronizer::callFunction
double callFunction(const V *vecValues, const V *vecDirection, V *gradVector, M *hessianMatrix, V *hessianEffect, double *extraOutput1, double *extraOutput2) const
Calls the scalar function which will be synchronized.
Definition: ScalarFunctionSynchronizer.C:61
QUESO::MetropolisHastingsSG::updateAdaptedCovMatrix
void updateAdaptedCovMatrix(const BaseVectorSequence< P_V, P_M > &subChain, unsigned int idOfFirstPositionInSubChain, double &lastChainSize, P_V &lastMean, P_M &lastAdaptedCovMatrix)
This method updates the adapted covariance matrix.
Definition: MetropolisHastingsSG.C:2212
QUESO::VectorSpace::newMatrix
M * newMatrix() const
Creates an empty matrix of size given by Map&amp; map. See template specialization.
QUESO::BaseEnvironment::numSubEnvironments
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:319
QUESO::BaseTKGroup::rv
virtual const GaussianVectorRV< V, M > & rv(unsigned int stageId)=0
Gaussian increment property to construct a transition kernel. See template specialization.
QUESO::MetropolisHastingsSG::alpha
double alpha(const MarkovChainPositionData< P_V > &x, const MarkovChainPositionData< P_V > &y, unsigned int xStageId, unsigned int yStageId, double *alphaQuotientPtr=NULL)
Calculates acceptance ration.
Definition: MetropolisHastingsSG.C:393
QUESO::MetropolisHastingsSG::m_idsOfUniquePositions
std::vector< unsigned int > m_idsOfUniquePositions
Definition: MetropolisHastingsSG.h:267
QUESO::BaseVectorSequence::getPositionValues
virtual void getPositionValues(unsigned int posId, V &vec) const =0
Gets the values of the sequence at position posId and stores them at vec. See template specialization...
QUESO::BaseScalarFunction::domainSet
const VectorSet< V, M > & domainSet() const
Access to the protected attribute m_domainSet: domain set of the scalar function. ...
Definition: ScalarFunction.C:50
QUESO::MhOptionsValues::m_amEpsilon
double m_amEpsilon
Regularisation parameter for the DRAM covariance matrix.
Definition: MetropolisHastingsSGOptions.h:206
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:436
QUESO::MetropolisHastingsSG::m_parameterEnabledStatus
std::vector< bool > m_parameterEnabledStatus
Definition: MetropolisHastingsSG.h:261
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::BaseVectorSequence::estimateConvBrooksGelman
virtual double estimateConvBrooksGelman(unsigned int initialPos, unsigned int numPos) const =0
Estimates convergence rate using Brooks &amp; Gelman method. See template specialization.
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:329
QUESO::VectorSpace::newVector
V * newVector() const
Creates an empty vector of size given by Map&amp; map. See template specialization.
UQ_FATAL_RC_MACRO
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
Definition: Defines.h:206
QUESO::MpiComm::NumProc
int NumProc() const
Returns total number of processes.
Definition: MpiComm.C:121
QUESO::MetropolisHastingsSG::m_alphaQuotients
std::vector< double > m_alphaQuotients
Definition: MetropolisHastingsSG.h:269
QUESO::MetropolisHastingsSG::acceptAlpha
bool acceptAlpha(double alpha)
Decides whether or not to accept alpha.
Definition: MetropolisHastingsSG.C:712
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
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 811 of file MetropolisHastingsSG.C.

References QUESO::BaseVectorSequence< V, M >::computeFilterParams(), QUESO::BaseVectorSequence< V, M >::filter(), QUESO::ScalarSequence< T >::filter(), QUESO::SequenceOfVectors< V, M >::getPositionValues(), QUESO::BaseVectorSequence< V, M >::name(), QUESO::FilePtrSetStruct::ofsVar, 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(), UQ_FATAL_RC_MACRO, UQ_FATAL_TEST_MACRO, UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, UQ_MH_SG_FILENAME_FOR_NO_FILE, 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().

815 {
816  if ((m_env.subDisplayFile() ) &&
817  (m_env.displayVerbosity() >= 5 ) &&
818  (m_optionsObj->m_ov.m_totallyMute == false)) {
819  *m_env.subDisplayFile() << "Entering MetropolisHastingsSG<P_V,P_M>::generateSequence()..."
820  << std::endl;
821  }
822 
823  UQ_FATAL_TEST_MACRO(m_vectorSpace.dimLocal() != workingChain.vectorSizeLocal(),
824  m_env.worldRank(),
825  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
826  "'m_vectorSpace' and 'workingChain' are related to vector spaces of different dimensions");
827 
828  //m_env.syncPrintDebugMsg("Entering MetropolisHastingsSG<P_V,P_M>::generateSequence()",2,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
829  MiscCheckTheParallelEnvironment<P_V,P_V>(m_initialPosition,
830  m_initialPosition);
831 
832  P_V valuesOf1stPosition(m_initialPosition);
833  int iRC = UQ_OK_RC;
834 
835  workingChain.setName(m_optionsObj->m_prefix + "rawChain");
836 
837  //****************************************************
838  // Generate chain
839  //****************************************************
840  if (m_optionsObj->m_ov.m_rawChainDataInputFileName == UQ_MH_SG_FILENAME_FOR_NO_FILE) {
841  generateFullChain(valuesOf1stPosition,
842  m_optionsObj->m_ov.m_rawChainSize,
843  workingChain,
844  workingLogLikelihoodValues,
845  workingLogTargetValues);
846  }
847  else {
848  readFullChain(m_optionsObj->m_ov.m_rawChainDataInputFileName,
849  m_optionsObj->m_ov.m_rawChainDataInputFileType,
850  m_optionsObj->m_ov.m_rawChainSize,
851  workingChain);
852  }
853 
854  //****************************************************
855  // Open generic output file
856  //****************************************************
857  if ((m_env.subDisplayFile() ) &&
858  (m_optionsObj->m_ov.m_totallyMute == false)) {
859  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
860  << ", prefix = " << m_optionsObj->m_prefix
861  << ", chain name = " << workingChain.name()
862  << ": about to try to open generic output file '" << m_optionsObj->m_ov.m_dataOutputFileName
863  << "." << UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT // Yes, always ".m"
864  << "', subId = " << m_env.subId()
865  << ", subenv is allowed to write (1/true or 0/false) = " << (m_optionsObj->m_ov.m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_ov.m_dataOutputAllowedSet.end())
866  << "..."
867  << std::endl;
868  }
869 
870  FilePtrSetStruct genericFilePtrSet;
871  m_env.openOutputFile(m_optionsObj->m_ov.m_dataOutputFileName,
872  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
873  m_optionsObj->m_ov.m_dataOutputAllowedSet,
874  false,
875  genericFilePtrSet);
876 
877  if ((m_env.subDisplayFile() ) &&
878  (m_optionsObj->m_ov.m_totallyMute == false)) {
879  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
880  << ", prefix = " << m_optionsObj->m_prefix
881  << ", raw chain name = " << workingChain.name()
882  << ": returned from opening generic output file '" << m_optionsObj->m_ov.m_dataOutputFileName
883  << "." << UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT // Yes, always ".m"
884  << "', subId = " << m_env.subId()
885  << std::endl;
886  }
887 
888  //****************************************************************************************
889  // Eventually:
890  // --> write raw chain
891  // --> compute statistics on it
892  //****************************************************************************************
893  if ((m_optionsObj->m_ov.m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) &&
894  (m_optionsObj->m_ov.m_totallyMute == false )) {
895 
896  // Take "sub" care of raw chain
897  if ((m_env.subDisplayFile() ) &&
898  (m_optionsObj->m_ov.m_totallyMute == false)) {
899  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
900  << ", prefix = " << m_optionsObj->m_prefix
901  << ", raw chain name = " << workingChain.name()
902  << ": about to try to write raw sub chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
903  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
904  << "', subId = " << m_env.subId()
905  << ", subenv is allowed to write 1/true or 0/false) = " << (m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet.end())
906  << "..."
907  << std::endl;
908  }
909 
910  if ((m_numPositionsNotSubWritten > 0 ) &&
911  (m_optionsObj->m_ov.m_rawChainDataOutputFileName != ".")) {
912  workingChain.subWriteContents(m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten,
913  m_numPositionsNotSubWritten,
914  m_optionsObj->m_ov.m_rawChainDataOutputFileName,
915  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
916  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
917  if ((m_env.subDisplayFile() ) &&
918  (m_optionsObj->m_ov.m_totallyMute == false)) {
919  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
920  << ": just wrote (per period request) remaining " << m_numPositionsNotSubWritten << " chain positions "
921  << ", " << m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten << " <= pos <= " << m_optionsObj->m_ov.m_rawChainSize - 1
922  << std::endl;
923  }
924 
925  if (workingLogLikelihoodValues) {
926  workingLogLikelihoodValues->subWriteContents(m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten,
927  m_numPositionsNotSubWritten,
928  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_likelihood",
929  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
930  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
931  }
932 
933  if (workingLogTargetValues) {
934  workingLogTargetValues->subWriteContents(m_optionsObj->m_ov.m_rawChainSize - m_numPositionsNotSubWritten,
935  m_numPositionsNotSubWritten,
936  m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_target",
937  m_optionsObj->m_ov.m_rawChainDataOutputFileType,
938  m_optionsObj->m_ov.m_rawChainDataOutputAllowedSet);
939  }
940 
941  m_numPositionsNotSubWritten = 0;
942  }
943 
944  // Compute raw sub MLE
945  if (workingLogLikelihoodValues) {
946  SequenceOfVectors<P_V,P_M> rawSubMLEpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawSubMLEseq");
947  double rawSubMLEvalue = workingChain.subPositionsOfMaximum(*workingLogLikelihoodValues,
948  rawSubMLEpositions);
949  UQ_FATAL_TEST_MACRO(rawSubMLEpositions.subSequenceSize() == 0,
950  m_env.worldRank(),
951  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
952  "rawSubMLEpositions.subSequenceSize() = 0");
953 
954  if ((m_env.subDisplayFile() ) &&
955  (m_optionsObj->m_ov.m_totallyMute == false)) {
956  P_V tmpVec(m_vectorSpace.zeroVector());
957  rawSubMLEpositions.getPositionValues(0,tmpVec);
958  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
959  << ": just computed MLE"
960  << ", rawSubMLEvalue = " << rawSubMLEvalue
961  << ", rawSubMLEpositions.subSequenceSize() = " << rawSubMLEpositions.subSequenceSize()
962  << ", rawSubMLEpositions[0] = " << tmpVec
963  << std::endl;
964  }
965  }
966 
967  // Compute raw sub MAP
968  if (workingLogTargetValues) {
969  SequenceOfVectors<P_V,P_M> rawSubMAPpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawSubMAPseq");
970  double rawSubMAPvalue = workingChain.subPositionsOfMaximum(*workingLogTargetValues,
971  rawSubMAPpositions);
972  UQ_FATAL_TEST_MACRO(rawSubMAPpositions.subSequenceSize() == 0,
973  m_env.worldRank(),
974  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
975  "rawSubMAPpositions.subSequenceSize() = 0");
976 
977  if ((m_env.subDisplayFile() ) &&
978  (m_optionsObj->m_ov.m_totallyMute == false)) {
979  P_V tmpVec(m_vectorSpace.zeroVector());
980  rawSubMAPpositions.getPositionValues(0,tmpVec);
981  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
982  << ": just computed MAP"
983  << ", rawSubMAPvalue = " << rawSubMAPvalue
984  << ", rawSubMAPpositions.subSequenceSize() = " << rawSubMAPpositions.subSequenceSize()
985  << ", rawSubMAPpositions[0] = " << tmpVec
986  << std::endl;
987  }
988  }
989 
990  if ((m_env.subDisplayFile() ) &&
991  (m_optionsObj->m_ov.m_totallyMute == false)) {
992  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
993  << ", prefix = " << m_optionsObj->m_prefix
994  << ", raw chain name = " << workingChain.name()
995  << ": returned from writing raw sub chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
996  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
997  << "', subId = " << m_env.subId()
998  << std::endl;
999  }
1000 
1001  // Take "unified" care of raw chain
1002  if ((m_env.subDisplayFile() ) &&
1003  (m_optionsObj->m_ov.m_totallyMute == false)) {
1004  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1005  << ", prefix = " << m_optionsObj->m_prefix
1006  << ", raw chain name = " << workingChain.name()
1007  << ": about to try to write raw unified chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
1008  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
1009  << "', subId = " << m_env.subId()
1010  << "..."
1011  << std::endl;
1012  }
1013 
1014  workingChain.unifiedWriteContents(m_optionsObj->m_ov.m_rawChainDataOutputFileName,
1015  m_optionsObj->m_ov.m_rawChainDataOutputFileType);
1016  if ((m_env.subDisplayFile() ) &&
1017  (m_optionsObj->m_ov.m_totallyMute == false)) {
1018  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1019  << ", prefix = " << m_optionsObj->m_prefix
1020  << ", raw chain name = " << workingChain.name()
1021  << ": returned from writing raw unified chain output file '" << m_optionsObj->m_ov.m_rawChainDataOutputFileName
1022  << "." << m_optionsObj->m_ov.m_rawChainDataOutputFileType
1023  << "', subId = " << m_env.subId()
1024  << std::endl;
1025  }
1026 
1027  if (workingLogLikelihoodValues) {
1028  workingLogLikelihoodValues->unifiedWriteContents(m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_likelihood",
1029  m_optionsObj->m_ov.m_rawChainDataOutputFileType);
1030  }
1031 
1032  if (workingLogTargetValues) {
1033  workingLogTargetValues->unifiedWriteContents(m_optionsObj->m_ov.m_rawChainDataOutputFileName + "_target",
1034  m_optionsObj->m_ov.m_rawChainDataOutputFileType);
1035  }
1036 
1037  // Compute raw unified MLE
1038  if (workingLogLikelihoodValues) {
1039  SequenceOfVectors<P_V,P_M> rawUnifiedMLEpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawUnifiedMLEseq");
1040  double rawUnifiedMLEvalue = workingChain.unifiedPositionsOfMaximum(*workingLogLikelihoodValues,
1041  rawUnifiedMLEpositions);
1042  UQ_FATAL_TEST_MACRO(rawUnifiedMLEpositions.subSequenceSize() == 0,
1043  m_env.worldRank(),
1044  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
1045  "rawUnifiedMLEpositions.subSequenceSize() = 0");
1046 
1047  if ((m_env.subDisplayFile() ) &&
1048  (m_optionsObj->m_ov.m_totallyMute == false)) {
1049  P_V tmpVec(m_vectorSpace.zeroVector());
1050  rawUnifiedMLEpositions.getPositionValues(0,tmpVec);
1051  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1052  << ": just computed MLE"
1053  << ", rawUnifiedMLEvalue = " << rawUnifiedMLEvalue
1054  << ", rawUnifiedMLEpositions.subSequenceSize() = " << rawUnifiedMLEpositions.subSequenceSize()
1055  << ", rawUnifiedMLEpositions[0] = " << tmpVec
1056  << std::endl;
1057  }
1058  }
1059 
1060  // Compute raw unified MAP
1061  if (workingLogTargetValues) {
1062  SequenceOfVectors<P_V,P_M> rawUnifiedMAPpositions(m_vectorSpace,0,m_optionsObj->m_prefix+"rawUnifiedMAPseq");
1063  double rawUnifiedMAPvalue = workingChain.unifiedPositionsOfMaximum(*workingLogTargetValues,
1064  rawUnifiedMAPpositions);
1065 
1066  UQ_FATAL_TEST_MACRO(rawUnifiedMAPpositions.subSequenceSize() == 0,
1067  m_env.worldRank(),
1068  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
1069  "rawUnifiedMAPpositions.subSequenceSize() = 0");
1070 
1071  if ((m_env.subDisplayFile() ) &&
1072  (m_optionsObj->m_ov.m_totallyMute == false)) {
1073  P_V tmpVec(m_vectorSpace.zeroVector());
1074  rawUnifiedMAPpositions.getPositionValues(0,tmpVec);
1075  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1076  << ": just computed MAP"
1077  << ", rawUnifiedMAPvalue = " << rawUnifiedMAPvalue
1078  << ", rawUnifiedMAPpositions.subSequenceSize() = " << rawUnifiedMAPpositions.subSequenceSize()
1079  << ", rawUnifiedMAPpositions[0] = " << tmpVec
1080  << std::endl;
1081  }
1082  }
1083  }
1084 
1085  // Take care of other aspects of raw chain
1086  if ((genericFilePtrSet.ofsVar ) &&
1087  (m_optionsObj->m_ov.m_totallyMute == false)) {
1088  // Write likelihoodValues and alphaValues, if they were requested by user
1089  iRC = writeInfo(workingChain,
1090  *genericFilePtrSet.ofsVar);
1091  UQ_FATAL_RC_MACRO(iRC,
1092  m_env.worldRank(),
1093  "MetropolisHastingsSG<P_V,P_M>::generateSequence()",
1094  "improper writeInfo() return");
1095  }
1096 
1097 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
1098  if (m_optionsObj->m_ov.m_rawChainComputeStats) {
1099  workingChain.computeStatistics(*m_optionsObj->m_rawChainStatisticalOptionsObj,
1100  genericFilePtrSet.ofsVar);
1101  }
1102 #endif
1103 
1104  //****************************************************************************************
1105  // Eventually:
1106  // --> filter the raw chain
1107  // --> write it
1108  // --> compute statistics on it
1109  //****************************************************************************************
1110  if (m_optionsObj->m_ov.m_filteredChainGenerate) {
1111  // Compute filter parameters
1112  unsigned int filterInitialPos = (unsigned int) (m_optionsObj->m_ov.m_filteredChainDiscardedPortion * (double) workingChain.subSequenceSize());
1113  unsigned int filterSpacing = m_optionsObj->m_ov.m_filteredChainLag;
1114  if (filterSpacing == 0) {
1115  workingChain.computeFilterParams(genericFilePtrSet.ofsVar,
1116  filterInitialPos,
1117  filterSpacing);
1118  }
1119 
1120  // Filter positions from the converged portion of the chain
1121  workingChain.filter(filterInitialPos,
1122  filterSpacing);
1123  workingChain.setName(m_optionsObj->m_prefix + "filtChain");
1124 
1125  if (workingLogLikelihoodValues) workingLogLikelihoodValues->filter(filterInitialPos,
1126  filterSpacing);
1127 
1128  if (workingLogTargetValues) workingLogTargetValues->filter(filterInitialPos,
1129  filterSpacing);
1130 
1131  // Write filtered chain
1132  if ((m_env.subDisplayFile() ) &&
1133  (m_optionsObj->m_ov.m_totallyMute == false)) {
1134  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1135  << ", prefix = " << m_optionsObj->m_prefix
1136  << ": checking necessity of opening output files for filtered chain " << workingChain.name()
1137  << "..."
1138  << std::endl;
1139  }
1140 
1141  // Take "sub" care of filtered chain
1142  if ((m_optionsObj->m_ov.m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) &&
1143  (m_optionsObj->m_ov.m_totallyMute == false )) {
1144  workingChain.subWriteContents(0,
1145  workingChain.subSequenceSize(),
1146  m_optionsObj->m_ov.m_filteredChainDataOutputFileName,
1147  m_optionsObj->m_ov.m_filteredChainDataOutputFileType,
1148  m_optionsObj->m_ov.m_filteredChainDataOutputAllowedSet);
1149  if ((m_env.subDisplayFile() ) &&
1150  (m_optionsObj->m_ov.m_totallyMute == false)) {
1151  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1152  << ", prefix = " << m_optionsObj->m_prefix
1153  << ": closed sub output file '" << m_optionsObj->m_ov.m_filteredChainDataOutputFileName
1154  << "' for filtered chain " << workingChain.name()
1155  << std::endl;
1156  }
1157 
1158  if (workingLogLikelihoodValues) {
1159  workingLogLikelihoodValues->subWriteContents(0,
1160  workingChain.subSequenceSize(),
1161  m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_likelihood",
1162  m_optionsObj->m_ov.m_filteredChainDataOutputFileType,
1163  m_optionsObj->m_ov.m_filteredChainDataOutputAllowedSet);
1164  }
1165 
1166  if (workingLogTargetValues) {
1167  workingLogTargetValues->subWriteContents(0,
1168  workingChain.subSequenceSize(),
1169  m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_target",
1170  m_optionsObj->m_ov.m_filteredChainDataOutputFileType,
1171  m_optionsObj->m_ov.m_filteredChainDataOutputAllowedSet);
1172  }
1173  }
1174 
1175  // Compute sub filtered MLE and sub filtered MAP
1176 
1177  // Take "unified" care of filtered chain
1178  if ((m_optionsObj->m_ov.m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) &&
1179  (m_optionsObj->m_ov.m_totallyMute == false )) {
1180  workingChain.unifiedWriteContents(m_optionsObj->m_ov.m_filteredChainDataOutputFileName,
1181  m_optionsObj->m_ov.m_filteredChainDataOutputFileType);
1182  if ((m_env.subDisplayFile() ) &&
1183  (m_optionsObj->m_ov.m_totallyMute == false)) {
1184  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1185  << ", prefix = " << m_optionsObj->m_prefix
1186  << ": closed unified output file '" << m_optionsObj->m_ov.m_filteredChainDataOutputFileName
1187  << "' for filtered chain " << workingChain.name()
1188  << std::endl;
1189  }
1190 
1191  if (workingLogLikelihoodValues) {
1192  workingLogLikelihoodValues->unifiedWriteContents(m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_likelihood",
1193  m_optionsObj->m_ov.m_filteredChainDataOutputFileType);
1194  }
1195 
1196  if (workingLogTargetValues) {
1197  workingLogTargetValues->unifiedWriteContents(m_optionsObj->m_ov.m_filteredChainDataOutputFileName + "_target",
1198  m_optionsObj->m_ov.m_filteredChainDataOutputFileType);
1199  }
1200  }
1201 
1202  // Compute unified filtered MLE and unified filtered MAP
1203 
1204  // Compute statistics
1205 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
1206  if (m_optionsObj->m_ov.m_filteredChainComputeStats) {
1207  workingChain.computeStatistics(*m_optionsObj->m_filteredChainStatisticalOptionsObj,
1208  genericFilePtrSet.ofsVar);
1209  }
1210 #endif
1211  }
1212 
1213  //****************************************************
1214  // Close generic output file
1215  //****************************************************
1216  if (genericFilePtrSet.ofsVar) {
1217  //genericFilePtrSet.ofsVar->close();
1218  delete genericFilePtrSet.ofsVar;
1219  if ((m_env.subDisplayFile() ) &&
1220  (m_optionsObj->m_ov.m_totallyMute == false)) {
1221  *m_env.subDisplayFile() << "In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1222  << ", prefix = " << m_optionsObj->m_prefix
1223  << ": closed generic output file '" << m_optionsObj->m_ov.m_dataOutputFileName
1224  << "' (chain name is " << workingChain.name()
1225  << ")"
1226  << std::endl;
1227  }
1228  }
1229 
1230  if ((m_env.subDisplayFile() ) &&
1231  (m_optionsObj->m_ov.m_totallyMute == false)) {
1232  *m_env.subDisplayFile() << std::endl;
1233  }
1234 
1235  //m_env.syncPrintDebugMsg("Leaving MetropolisHastingsSG<P_V,P_M>::generateSequence()",2,3000000,m_env.fullComm()); // Dangerous to barrier on fullComm ... // KAUST
1236 
1237  if ((m_env.subDisplayFile() ) &&
1238  (m_env.displayVerbosity() >= 5 ) &&
1239  (m_optionsObj->m_ov.m_totallyMute == false)) {
1240  *m_env.subDisplayFile() << "Leaving MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1241  << std::endl;
1242  }
1243 
1244  return;
1245 }
QUESO::BaseVectorSequence::name
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
Definition: VectorSequence.C:121
QUESO::MetropolisHastingsSG::generateFullChain
void generateFullChain(const P_V &valuesOf1stPosition, unsigned int chainSize, BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
This method generates the chain.
Definition: MetropolisHastingsSG.C:1270
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::ScalarSequence::subWriteContents
void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Writes the sub-sequence to a file.
Definition: ScalarSequence.C:2742
QUESO::MetropolisHastingsSG::writeInfo
int writeInfo(const BaseVectorSequence< P_V, P_M > &workingChain, std::ofstream &ofsvar) const
Writes information about the Markov chain in a file.
Definition: MetropolisHastingsSG.C:726
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::BaseVectorSequence::setName
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
Definition: VectorSequence.C:128
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:218
QUESO::BaseVectorSequence::filter
virtual void filter(unsigned int initialPos, unsigned int spacing)=0
Filters positions in the sequence of vectors, starting at initialPos, and with spacing given by spaci...
QUESO::BaseVectorSequence::subWriteContents
virtual void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const =0
Writes info of the sub-sequence to a file. See template specialization.
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2575
UQ_MH_SG_FILENAME_FOR_NO_FILE
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
Definition: MetropolisHastingsSGOptions.h:35
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:75
QUESO::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2799
QUESO::BaseVectorSequence::vectorSizeLocal
unsigned int vectorSizeLocal() const
Local dimension (size) of the vector space.
Definition: VectorSequence.C:100
QUESO::MhOptionsValues::m_rawChainDataOutputAllowedSet
std::set< unsigned int > m_rawChainDataOutputAllowedSet
Definition: MetropolisHastingsSGOptions.h:155
QUESO::BaseVectorSequence::unifiedWriteContents
virtual void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const =0
Writes info of the unified sequence to a file. See template specialization.
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
QUESO::MhOptionsValues::m_filteredChainDataOutputAllowedSet
std::set< unsigned int > m_filteredChainDataOutputAllowedSet
Definition: MetropolisHastingsSGOptions.h:166
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::MetropolisHastingsSG::readFullChain
void readFullChain(const std::string &inputFileName, const std::string &inputFileType, unsigned int chainSize, BaseVectorSequence< P_V, P_M > &workingChain)
This method reads the chain contents.
Definition: MetropolisHastingsSG.C:1258
UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:88
QUESO::MhOptionsValues::m_dataOutputAllowedSet
std::set< unsigned int > m_dataOutputAllowedSet
Definition: MetropolisHastingsSGOptions.h:137
QUESO::BaseVectorSequence::computeFilterParams
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
Definition: VectorSequence.C:613
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::BaseEnvironment::displayVerbosity
unsigned int displayVerbosity() const
Definition: Environment.C:436
QUESO::BaseVectorSequence::subPositionsOfMaximum
double subPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
Definition: VectorSequence.C:386
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::BaseEnvironment::subId
unsigned int subId() const
Access function to the number of each sub-environment Id: m_subId.
Definition: Environment.C:329
QUESO::BaseEnvironment::openOutputFile
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:516
UQ_FATAL_RC_MACRO
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
Definition: Defines.h:206
QUESO::BaseVectorSequence::unifiedPositionsOfMaximum
double unifiedPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
Definition: VectorSequence.C:433
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
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 1250 of file MetropolisHastingsSG.C.

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

1251 {
1252  info = m_rawChainInfo;
1253  return;
1254 }
QUESO::MetropolisHastingsSG::m_rawChainInfo
MHRawChainInfoStruct m_rawChainInfo
Definition: MetropolisHastingsSG.h:275
template<class P_V, class P_M>
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 1258 of file MetropolisHastingsSG.C.

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

1263 {
1264  workingChain.unifiedReadContents(inputFileName,inputFileType,chainSize);
1265  return;
1266 }
QUESO::BaseVectorSequence::unifiedReadContents
virtual void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)=0
Reads info of the unified sequence from a file. See template specialization.
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 2212 of file MetropolisHastingsSG.C.

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

2218 {
2219  double doubleSubChainSize = (double) partialChain.subSequenceSize();
2220  if (lastChainSize == 0) {
2221  UQ_FATAL_TEST_MACRO(partialChain.subSequenceSize() < 2,
2222  m_env.worldRank(),
2223  "MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix()",
2224  "'partialChain.subSequenceSize()' should be >= 2");
2225 
2226 #if 1 // prudenci-2012-07-06
2227  lastMean = partialChain.subMeanPlain();
2228 #else
2229  partialChain.subMeanExtra(0,partialChain.subSequenceSize(),lastMean);
2230 #endif
2231 
2232  P_V tmpVec(m_vectorSpace.zeroVector());
2233  lastAdaptedCovMatrix = -doubleSubChainSize * matrixProduct(lastMean,lastMean);
2234  for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2235  partialChain.getPositionValues(i,tmpVec);
2236  lastAdaptedCovMatrix += matrixProduct(tmpVec,tmpVec);
2237  }
2238  lastAdaptedCovMatrix /= (doubleSubChainSize - 1.); // That is why partialChain size must be >= 2
2239  }
2240  else {
2241  UQ_FATAL_TEST_MACRO(partialChain.subSequenceSize() < 1,
2242  m_env.worldRank(),
2243  "MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix()",
2244  "'partialChain.subSequenceSize()' should be >= 1");
2245 
2246  UQ_FATAL_TEST_MACRO(idOfFirstPositionInSubChain < 1,
2247  m_env.worldRank(),
2248  "MetropolisHastingsSG<P_V,P_M>::updateAdaptedCovMatrix()",
2249  "'idOfFirstPositionInSubChain' should be >= 1");
2250 
2251  P_V tmpVec (m_vectorSpace.zeroVector());
2252  P_V diffVec(m_vectorSpace.zeroVector());
2253  for (unsigned int i = 0; i < partialChain.subSequenceSize(); ++i) {
2254  double doubleCurrentId = (double) (idOfFirstPositionInSubChain+i);
2255  partialChain.getPositionValues(i,tmpVec);
2256  diffVec = tmpVec - lastMean;
2257 
2258  double ratio1 = (1. - 1./doubleCurrentId); // That is why idOfFirstPositionInSubChain must be >= 1
2259  double ratio2 = (1./(1.+doubleCurrentId));
2260  lastAdaptedCovMatrix = ratio1 * lastAdaptedCovMatrix + ratio2 * matrixProduct(diffVec,diffVec);
2261  lastMean += ratio2 * diffVec;
2262  }
2263  }
2264  lastChainSize += doubleSubChainSize;
2265 
2266  if (m_numDisabledParameters > 0) { // gpmsa2
2267  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2268  if (m_parameterEnabledStatus[paramId] == false) {
2269  for (unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2270  lastAdaptedCovMatrix(i,paramId) = 0.;
2271  }
2272  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2273  lastAdaptedCovMatrix(paramId,j) = 0.;
2274  }
2275  lastAdaptedCovMatrix(paramId,paramId) = 1.;
2276  }
2277  }
2278  }
2279 
2280  return;
2281 }
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the process world rank.
Definition: Environment.C:235
QUESO::VectorSpace::zeroVector
const V & zeroVector() const
Returns a vector filled with zeros.
Definition: VectorSpace.C:218
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2409
QUESO::VectorSpace::dimLocal
unsigned int dimLocal() const
Definition: VectorSpace.C:199
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:222
QUESO::MetropolisHastingsSG::m_parameterEnabledStatus
std::vector< bool > m_parameterEnabledStatus
Definition: MetropolisHastingsSG.h:261
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
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 726 of file MetropolisHastingsSG.C.

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

729 {
730  if ((m_env.subDisplayFile() ) &&
731  (m_optionsObj->m_ov.m_totallyMute == false)) {
732  *m_env.subDisplayFile() << "\n"
733  << "\n-----------------------------------------------------"
734  << "\n Writing more information about the Markov chain " << workingChain.name() << " to output file ..."
735  << "\n-----------------------------------------------------"
736  << "\n"
737  << std::endl;
738  }
739 
740  int iRC = UQ_OK_RC;
741 
742  if (m_optionsObj->m_ov.m_rawChainGenerateExtra) {
743  // Write m_logTargets
744  ofsvar << m_optionsObj->m_prefix << "logTargets_sub" << m_env.subIdString() << " = zeros(" << m_logTargets.size()
745  << "," << 1
746  << ");"
747  << std::endl;
748  ofsvar << m_optionsObj->m_prefix << "logTargets_sub" << m_env.subIdString() << " = [";
749  for (unsigned int i = 0; i < m_logTargets.size(); ++i) {
750  ofsvar << m_logTargets[i]
751  << std::endl;
752  }
753  ofsvar << "];\n";
754 
755  // Write m_alphaQuotients
756  ofsvar << m_optionsObj->m_prefix << "alphaQuotients_sub" << m_env.subIdString() << " = zeros(" << m_alphaQuotients.size()
757  << "," << 1
758  << ");"
759  << std::endl;
760  ofsvar << m_optionsObj->m_prefix << "alphaQuotients_sub" << m_env.subIdString() << " = [";
761  for (unsigned int i = 0; i < m_alphaQuotients.size(); ++i) {
762  ofsvar << m_alphaQuotients[i]
763  << std::endl;
764  }
765  ofsvar << "];\n";
766  }
767 
768  // Write number of rejections
769  ofsvar << m_optionsObj->m_prefix << "rejected = " << (double) m_rawChainInfo.numRejections/(double) (workingChain.subSequenceSize()-1)
770  << ";\n"
771  << std::endl;
772 
773  if (false) { // Don't see need for code below. Let it there though, compiling, in case it is needed in the future.
774  // Write names of components
775  ofsvar << m_optionsObj->m_prefix << "componentNames = {";
776  m_vectorSpace.printComponentsNames(ofsvar,false);
777  ofsvar << "};\n";
778 
779  // Write number of out of target support
780  ofsvar << m_optionsObj->m_prefix << "outTargetSupport = " << (double) m_rawChainInfo.numOutOfTargetSupport/(double) (workingChain.subSequenceSize()-1)
781  << ";\n"
782  << std::endl;
783 
784  // Write chain run time
785  ofsvar << m_optionsObj->m_prefix << "runTime = " << m_rawChainInfo.runTime
786  << ";\n"
787  << std::endl;
788  }
789 
790  if ((m_env.subDisplayFile() ) &&
791  (m_optionsObj->m_ov.m_totallyMute == false)) {
792  *m_env.subDisplayFile() << "\n-----------------------------------------------------"
793  << "\n Finished writing more information about the Markov chain " << workingChain.name()
794  << "\n-----------------------------------------------------"
795  << "\n"
796  << std::endl;
797  }
798 
799  return iRC;
800 }
QUESO::BaseEnvironment::subIdString
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:335
QUESO::BaseVectorSequence::name
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
Definition: VectorSequence.C:121
QUESO::MetropolisHastingsSG::m_logTargets
std::vector< double > m_logTargets
Definition: MetropolisHastingsSG.h:268
QUESO::MetropolisHastingsSG::m_vectorSpace
const VectorSpace< P_V, P_M > & m_vectorSpace
Definition: MetropolisHastingsSG.h:255
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:305
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:75
QUESO::MetropolisHastingsSG::m_rawChainInfo
MHRawChainInfoStruct m_rawChainInfo
Definition: MetropolisHastingsSG.h:275
QUESO::MetropolisHastingsSG::m_optionsObj
MetropolisHastingsSGOptions * m_optionsObj
Definition: MetropolisHastingsSG.h:278
QUESO::VectorSpace::printComponentsNames
void printComponentsNames(std::ostream &os, bool printHorizontally) const
Prints the local component names.
Definition: VectorSpace.C:318
QUESO::MetropolisHastingsSGOptions::m_ov
MhOptionsValues m_ov
This class is where the actual options are stored.
Definition: MetropolisHastingsSGOptions.h:265
QUESO::MetropolisHastingsSG::m_env
const BaseEnvironment & m_env
Definition: MetropolisHastingsSG.h:254
QUESO::MetropolisHastingsSG::m_alphaQuotients
std::vector< double > m_alphaQuotients
Definition: MetropolisHastingsSG.h:269
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.

Friends And Related Function Documentation

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

Definition at line 186 of file MetropolisHastingsSG.h.

188  {
189  obj.print(os);
190 
191  return os;
192  }

Member Data Documentation

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

Definition at line 269 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
MhOptionsValues QUESO::MetropolisHastingsSG< P_V, P_M >::m_alternativeOptionsValues
private

Definition at line 277 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
const BaseEnvironment& QUESO::MetropolisHastingsSG< P_V, P_M >::m_env
private

Definition at line 254 of file MetropolisHastingsSG.h.

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

template<class P_V, class P_M>
std::vector<unsigned int> QUESO::MetropolisHastingsSG< P_V, P_M >::m_idsOfUniquePositions
private

Definition at line 267 of file MetropolisHastingsSG.h.

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

Definition at line 257 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
P_M QUESO::MetropolisHastingsSG< P_V, P_M >::m_initialProposalCovMatrix
private

Definition at line 258 of file MetropolisHastingsSG.h.

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

template<class P_V, class P_M>
P_M* QUESO::MetropolisHastingsSG< P_V, P_M >::m_lastAdaptedCovMatrix
private

Definition at line 272 of file MetropolisHastingsSG.h.

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

Definition at line 270 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
P_V* QUESO::MetropolisHastingsSG< P_V, P_M >::m_lastMean
private

Definition at line 271 of file MetropolisHastingsSG.h.

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

Definition at line 268 of file MetropolisHastingsSG.h.

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

Definition at line 259 of file MetropolisHastingsSG.h.

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

Definition at line 260 of file MetropolisHastingsSG.h.

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

Definition at line 273 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
MetropolisHastingsSGOptions* QUESO::MetropolisHastingsSG< P_V, P_M >::m_optionsObj
private

Definition at line 278 of file MetropolisHastingsSG.h.

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

template<class P_V, class P_M>
std::vector<bool> QUESO::MetropolisHastingsSG< P_V, P_M >::m_parameterEnabledStatus
private

Definition at line 261 of file MetropolisHastingsSG.h.

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

Definition at line 265 of file MetropolisHastingsSG.h.

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

Definition at line 275 of file MetropolisHastingsSG.h.

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

Definition at line 266 of file MetropolisHastingsSG.h.

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

Definition at line 256 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
const ScalarFunctionSynchronizer<P_V,P_M>* QUESO::MetropolisHastingsSG< P_V, P_M >::m_targetPdfSynchronizer
private

Definition at line 262 of file MetropolisHastingsSG.h.

template<class P_V, class P_M>
BaseTKGroup<P_V,P_M>* QUESO::MetropolisHastingsSG< P_V, P_M >::m_tk
private

Definition at line 264 of file MetropolisHastingsSG.h.

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

Definition at line 255 of file MetropolisHastingsSG.h.

The documentation for this class was generated from the following files:
Generated on Thu Apr 23 2015 19:18:36 for queso-0.50.1 by  doxygen 1.8.5