25 #include <queso/asserts.h>
26 #include <queso/MetropolisHastingsSG.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
30 #include <queso/HessianCovMatricesTKGroup.h>
31 #include <queso/ScaledCovMatrixTKGroup.h>
32 #include <queso/TransformedScaledCovMatrixTKGroup.h>
34 #include <queso/InvLogitGaussianJointPdf.h>
36 #include <queso/GaussianJointPdf.h>
38 #include <queso/TransitionKernelFactory.h>
39 #include <queso/AlgorithmFactory.h>
127 "MHRawChainInfoStruct::mpiSum()",
128 "failed MPI.Allreduce() for sum of doubles");
131 "MHRawChainInfoStruct::mpiSum()",
132 "failed MPI.Allreduce() for sum of unsigned ints");
138 template<
class P_V,
class P_M>
146 m_env (sourceRv.env()),
147 m_vectorSpace (sourceRv.imageSet().vectorSpace()),
148 m_targetPdf (sourceRv.pdf()),
149 m_initialPosition (initialPosition),
150 m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
151 m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
152 m_numDisabledParameters (0),
153 m_parameterEnabledStatus (m_vectorSpace.dimLocal(),
true),
157 m_positionIdForDebugging (0),
158 m_stageIdForDebugging (0),
159 m_idsOfUniquePositions (0),
161 m_alphaQuotients (0),
164 m_lastAdaptedCovMatrix (NULL),
165 m_numPositionsNotSubWritten (0),
166 m_optionsObj (alternativeOptionsValues),
167 m_computeInitialPriorAndLikelihoodValues(
true),
168 m_initialLogPriorValue (0.),
169 m_initialLogLikelihoodValue (0.),
170 m_userDidNotProvideOptions(
false)
172 if (inputProposalCovMatrix != NULL) {
173 m_initialProposalCovMatrix = *inputProposalCovMatrix;
176 if (m_optionsObj == NULL) {
181 m_optionsObj = tempOptions;
184 m_userDidNotProvideOptions =
true;
187 if (m_optionsObj->m_help !=
"") {
188 if (m_env.subDisplayFile() && !m_optionsObj->m_totallyMute) {
189 *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
193 if ((m_env.subDisplayFile() ) &&
194 (m_optionsObj->m_totallyMute ==
false)) {
195 *m_env.subDisplayFile() <<
"Entering MetropolisHastingsSG<P_V,P_M>::constructor(1)"
196 <<
": prefix = " << prefix
197 <<
", alternativeOptionsValues = " << alternativeOptionsValues
198 <<
", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
199 <<
", m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
203 queso_require_equal_to_msg(sourceRv.imageSet().vectorSpace().dimLocal(), initialPosition.sizeLocal(),
"'sourceRv' and 'initialPosition' should have equal dimensions");
205 if (inputProposalCovMatrix) {
206 queso_require_equal_to_msg(sourceRv.imageSet().vectorSpace().dimLocal(), inputProposalCovMatrix->numRowsLocal(),
"'sourceRv' and 'inputProposalCovMatrix' should have equal dimensions");
207 queso_require_equal_to_msg(inputProposalCovMatrix->numCols(), inputProposalCovMatrix->numRowsGlobal(),
"'inputProposalCovMatrix' should be a square matrix");
212 if ((m_env.subDisplayFile() ) &&
213 (m_optionsObj->m_totallyMute ==
false)) {
214 *m_env.subDisplayFile() <<
"Leaving MetropolisHastingsSG<P_V,P_M>::constructor(1)"
218 template<
class P_V,
class P_M>
228 m_env (sourceRv.env()),
229 m_vectorSpace (sourceRv.imageSet().vectorSpace()),
230 m_targetPdf (sourceRv.pdf()),
231 m_initialPosition (initialPosition),
232 m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
233 m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
234 m_numDisabledParameters (0),
235 m_parameterEnabledStatus (m_vectorSpace.dimLocal(),
true),
239 m_positionIdForDebugging (0),
240 m_stageIdForDebugging (0),
241 m_idsOfUniquePositions (0),
243 m_alphaQuotients (0),
246 m_lastAdaptedCovMatrix (NULL),
247 m_numPositionsNotSubWritten (0),
248 m_optionsObj (alternativeOptionsValues),
249 m_computeInitialPriorAndLikelihoodValues(
false),
250 m_initialLogPriorValue (initialLogPrior),
251 m_initialLogLikelihoodValue (initialLogLikelihood),
252 m_userDidNotProvideOptions(
false)
254 if (inputProposalCovMatrix != NULL) {
255 m_initialProposalCovMatrix = *inputProposalCovMatrix;
259 if (m_optionsObj == NULL) {
264 m_optionsObj = tempOptions;
266 m_userDidNotProvideOptions =
true;
269 if (m_optionsObj->m_help !=
"") {
270 if (m_env.subDisplayFile() && !m_optionsObj->m_totallyMute) {
271 *m_env.subDisplayFile() << (*m_optionsObj) << std::endl;
275 if ((m_env.subDisplayFile() ) &&
276 (m_optionsObj->m_totallyMute ==
false)) {
277 *m_env.subDisplayFile() <<
"Entering MetropolisHastingsSG<P_V,P_M>::constructor(2)"
278 <<
": prefix = " << prefix
279 <<
", alternativeOptionsValues = " << alternativeOptionsValues
280 <<
", m_env.optionsInputFileName() = " << m_env.optionsInputFileName()
281 <<
", m_initialProposalCovMatrix = " << m_initialProposalCovMatrix
285 queso_require_equal_to_msg(sourceRv.imageSet().vectorSpace().dimLocal(), initialPosition.sizeLocal(),
"'sourceRv' and 'initialPosition' should have equal dimensions");
287 if (inputProposalCovMatrix) {
288 queso_require_equal_to_msg(sourceRv.imageSet().vectorSpace().dimLocal(), inputProposalCovMatrix->numRowsLocal(),
"'sourceRv' and 'inputProposalCovMatrix' should have equal dimensions");
289 queso_require_equal_to_msg(inputProposalCovMatrix->numCols(), inputProposalCovMatrix->numRowsGlobal(),
"'inputProposalCovMatrix' should be a square matrix");
294 if ((m_env.subDisplayFile() ) &&
295 (m_optionsObj->m_totallyMute ==
false)) {
296 *m_env.subDisplayFile() <<
"Leaving MetropolisHastingsSG<P_V,P_M>::constructor(2)"
301 template<
class P_V,
class P_M>
305 const P_V& initialPosition,
306 const P_M* inputProposalCovMatrix)
308 m_env (sourceRv.env()),
309 m_vectorSpace (sourceRv.imageSet().vectorSpace()),
310 m_targetPdf (sourceRv.pdf()),
311 m_initialPosition (initialPosition),
312 m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
313 m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
314 m_numDisabledParameters (0),
315 m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true),
319 m_positionIdForDebugging (0),
320 m_stageIdForDebugging (0),
321 m_idsOfUniquePositions (0),
323 m_alphaQuotients (0),
326 m_lastAdaptedCovMatrix (NULL),
327 m_computeInitialPriorAndLikelihoodValues(true),
328 m_initialLogPriorValue (0.),
329 m_initialLogLikelihoodValue (0.),
330 m_userDidNotProvideOptions(true)
342 if (inputProposalCovMatrix != NULL) {
366 template<
class P_V,
class P_M>
370 const P_V& initialPosition,
371 double initialLogPrior,
372 double initialLogLikelihood,
373 const P_M* inputProposalCovMatrix)
375 m_env (sourceRv.env()),
376 m_vectorSpace (sourceRv.imageSet().vectorSpace()),
377 m_targetPdf (sourceRv.pdf()),
378 m_initialPosition (initialPosition),
379 m_initialProposalCovMatrix (m_vectorSpace.zeroVector()),
380 m_nullInputProposalCovMatrix(inputProposalCovMatrix == NULL),
381 m_numDisabledParameters (0),
382 m_parameterEnabledStatus (m_vectorSpace.dimLocal(),true),
386 m_positionIdForDebugging (0),
387 m_stageIdForDebugging (0),
388 m_idsOfUniquePositions (0),
390 m_alphaQuotients (0),
393 m_lastAdaptedCovMatrix (NULL),
394 m_computeInitialPriorAndLikelihoodValues(false),
395 m_initialLogPriorValue (initialLogPrior),
396 m_initialLogLikelihoodValue (initialLogLikelihood),
397 m_userDidNotProvideOptions(true)
409 if (inputProposalCovMatrix != NULL) {
433 template<
class P_V,
class P_M>
441 if (m_lastAdaptedCovMatrix)
delete m_lastAdaptedCovMatrix;
442 if (m_lastMean)
delete m_lastMean;
444 m_rawChainInfo.reset();
445 m_alphaQuotients.clear();
446 m_logTargets.clear();
447 m_numDisabledParameters = 0;
448 m_parameterEnabledStatus.clear();
449 m_positionIdForDebugging = 0;
450 m_stageIdForDebugging = 0;
451 m_idsOfUniquePositions.clear();
453 if (m_targetPdfSynchronizer)
delete m_targetPdfSynchronizer;
458 if (m_optionsObj && m_userDidNotProvideOptions) {
469 template<
class P_V,
class P_M>
476 if ((m_env.subDisplayFile() ) &&
477 (m_optionsObj->m_totallyMute ==
false)) {
478 *m_env.subDisplayFile() <<
"Entering MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
482 if (m_optionsObj->m_initialPositionDataInputFileName !=
".") {
483 std::set<unsigned int> tmpSet;
484 tmpSet.insert(m_env.subId());
485 m_initialPosition.subReadContents((m_optionsObj->m_initialPositionDataInputFileName+
"_sub"+m_env.subIdString()),
486 m_optionsObj->m_initialPositionDataInputFileType,
488 if ((m_env.subDisplayFile() ) &&
489 (m_optionsObj->m_totallyMute ==
false)) {
490 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
491 <<
": just read initial position contents = " << m_initialPosition
496 if (m_optionsObj->m_parameterDisabledSet.size() > 0) {
497 for (std::set<unsigned int>::iterator setIt = m_optionsObj->m_parameterDisabledSet.begin(); setIt != m_optionsObj->m_parameterDisabledSet.end(); ++setIt) {
498 unsigned int paramId = *setIt;
499 if (paramId < m_vectorSpace.dimLocal()) {
500 m_numDisabledParameters++;
501 m_parameterEnabledStatus[paramId] =
false;
506 std::vector<double> drScalesAll(m_optionsObj->m_drScalesForExtraStages.size()+1,1.);
507 for (
unsigned int i = 1; i < (m_optionsObj->m_drScalesForExtraStages.size()+1); ++i) {
508 drScalesAll[i] = m_optionsObj->m_drScalesForExtraStages[i-1];
514 msg =
"The doLogitTransform option is deprecated. ";
515 msg +=
"Use both ip_mh_algorithm and ip_mh_tk instead.";
519 if (m_optionsObj->m_initialProposalCovMatrixDataInputFileName !=
".") {
520 std::set<unsigned int> tmpSet;
521 tmpSet.insert(m_env.subId());
522 m_initialProposalCovMatrix.subReadContents((m_optionsObj->m_initialProposalCovMatrixDataInputFileName+
"_sub"+m_env.subIdString()),
523 m_optionsObj->m_initialProposalCovMatrixDataInputFileType,
525 if ((m_env.subDisplayFile() ) &&
526 (m_optionsObj->m_totallyMute ==
false)) {
527 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::commonConstructor()"
528 <<
": just read initial proposal cov matrix contents = " << m_initialProposalCovMatrix
533 queso_require_msg(!(m_nullInputProposalCovMatrix),
"proposal cov matrix should have been passed by user, since, according to the input algorithm options, local Hessians will not be used in the proposal");
539 if ((m_optionsObj->m_algorithm ==
"logit_random_walk") &&
540 (m_optionsObj->m_tk ==
"logit_random_walk") &&
541 m_optionsObj->m_doLogitTransform) {
543 transformInitialCovMatrixToGaussianSpace(
560 template<
class P_V,
class P_M>
564 const std::vector<unsigned int >& inputTKStageIds)
572 unsigned int inputSize = inputPositionsData.size();
573 if ((m_env.subDisplayFile() ) &&
574 (m_env.displayVerbosity() >= 10 ) &&
575 (m_optionsObj->m_totallyMute ==
false)) {
576 *m_env.subDisplayFile() <<
"Entering MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
577 <<
", inputSize = " << inputSize
581 queso_require_equal_to_msg(inputSize, inputPositionsData.size(),
"inputPositionsData and inputTKStageIds have different lengths");
584 if (inputPositionsData[0 ]->outOfTargetSupport())
return 0.;
585 if (inputPositionsData[inputSize-1]->outOfTargetSupport())
return 0.;
587 if ((inputPositionsData[0]->logTarget() == -INFINITY ) ||
588 (inputPositionsData[0]->logTarget() == INFINITY ) ||
589 (
queso_isnan(inputPositionsData[0]->logTarget()) )) {
590 std::cerr <<
"WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
591 <<
", worldRank " << m_env.worldRank()
592 <<
", fullRank " << m_env.fullRank()
593 <<
", subEnvironment " << m_env.subId()
594 <<
", subRank " << m_env.subRank()
595 <<
", inter0Rank " << m_env.inter0Rank()
596 <<
", positionId = " << m_positionIdForDebugging
597 <<
", stageId = " << m_stageIdForDebugging
598 <<
": inputSize = " << inputSize
599 <<
", inputPositionsData[0]->logTarget() = " << inputPositionsData[0]->logTarget()
600 <<
", [0]->values() = " << inputPositionsData[0]->vecValues()
601 <<
", [inputSize - 1]->values() = " << inputPositionsData[inputSize-1]->vecValues()
605 else if ((inputPositionsData[inputSize - 1]->logTarget() == -INFINITY ) ||
606 (inputPositionsData[inputSize - 1]->logTarget() == INFINITY ) ||
607 (
queso_isnan(inputPositionsData[inputSize - 1]->logTarget()) )) {
608 std::cerr <<
"WARNING In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
609 <<
", worldRank " << m_env.worldRank()
610 <<
", fullRank " << m_env.fullRank()
611 <<
", subEnvironment " << m_env.subId()
612 <<
", subRank " << m_env.subRank()
613 <<
", inter0Rank " << m_env.inter0Rank()
614 <<
", positionId = " << m_positionIdForDebugging
615 <<
", stageId = " << m_stageIdForDebugging
616 <<
": inputSize = " << inputSize
617 <<
", inputPositionsData[inputSize - 1]->logTarget() = " << inputPositionsData[inputSize-1]->logTarget()
618 <<
", [0]->values() = " << inputPositionsData[0]->vecValues()
619 <<
", [inputSize - 1]->values() = " << inputPositionsData[inputSize-1]->vecValues()
625 if (inputSize == 2) {
626 const P_V & tk_pos_x = m_tk->preComputingPosition(inputTKStageIds[inputSize-1]);
627 const P_V & tk_pos_y = m_tk->preComputingPosition(inputTKStageIds[0]);
628 return this->m_algorithm->acceptance_ratio(
629 *(inputPositionsData[0]),
630 *(inputPositionsData[inputSize - 1]),
636 std::vector<MarkovChainPositionData<P_V>*> positionsData (inputSize,NULL);
637 std::vector<MarkovChainPositionData<P_V>*> backwardPositionsData (inputSize,NULL);
639 std::vector<unsigned int > tkStageIds (inputSize,0);
640 std::vector<unsigned int > backwardTKStageIds (inputSize,0);
642 std::vector<unsigned int > tkStageIdsLess1(inputSize,0);
643 std::vector<unsigned int > backwardTKStageIdsLess1(inputSize,0);
645 for (
unsigned int i = 0; i < inputSize; ++i) {
646 positionsData [i] = inputPositionsData[i];
647 backwardPositionsData [i] = inputPositionsData[inputSize-i-1];
649 tkStageIds [i] = inputTKStageIds [i];
650 backwardTKStageIds [i] = inputTKStageIds [inputSize-i-1];
652 tkStageIdsLess1[i] = inputTKStageIds [i];
653 backwardTKStageIdsLess1[i] = inputTKStageIds [inputSize-i-1];
656 tkStageIdsLess1.pop_back();
657 backwardTKStageIdsLess1.pop_back();
660 double logNumerator = 0.;
661 double logDenominator = 0.;
662 double alphasNumerator = 1.;
663 double alphasDenominator = 1.;
666 const P_V& _lastTKPosition = m_tk->preComputingPosition( tkStageIds[inputSize-1]);
667 const P_V& _lastBackwardTKPosition = m_tk->preComputingPosition(backwardTKStageIds[inputSize-1]);
669 double numContrib = m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(_lastBackwardTKPosition,NULL,NULL,NULL,NULL);
671 double denContrib = m_tk->rv(tkStageIdsLess1).pdf().lnValue(_lastTKPosition,NULL,NULL,NULL,NULL);
672 if ((m_env.subDisplayFile() ) &&
673 (m_env.displayVerbosity() >= 10 ) &&
674 (m_optionsObj->m_totallyMute ==
false)) {
675 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
676 <<
", inputSize = " << inputSize
678 <<
": numContrib = " << numContrib
679 <<
", denContrib = " << denContrib
683 logNumerator += numContrib;
684 logDenominator += denContrib;
686 for (
unsigned int i = 0; i < (inputSize-2); ++i) {
687 positionsData.pop_back();
688 backwardPositionsData.pop_back();
690 const P_V& lastTKPosition = m_tk->preComputingPosition( tkStageIds[inputSize-2-i]);
691 const P_V& lastBackwardTKPosition = m_tk->preComputingPosition(backwardTKStageIds[inputSize-2-i]);
693 tkStageIds.pop_back();
694 backwardTKStageIds.pop_back();
696 tkStageIdsLess1.pop_back();
697 backwardTKStageIdsLess1.pop_back();
699 numContrib = m_tk->rv(backwardTKStageIdsLess1).pdf().lnValue(lastBackwardTKPosition,NULL,NULL,NULL,NULL);
701 denContrib = m_tk->rv(tkStageIdsLess1).pdf().lnValue(lastTKPosition,NULL,NULL,NULL,NULL);
702 if ((m_env.subDisplayFile() ) &&
703 (m_env.displayVerbosity() >= 10 ) &&
704 (m_optionsObj->m_totallyMute ==
false)) {
705 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
706 <<
", inputSize = " << inputSize
707 <<
", in loop, i = " << i
708 <<
": numContrib = " << numContrib
709 <<
", denContrib = " << denContrib
713 logNumerator += numContrib;
714 logDenominator += denContrib;
716 alphasNumerator *= (1. - this->alpha(backwardPositionsData,backwardTKStageIds));
717 alphasDenominator *= (1. - this->alpha( positionsData, tkStageIds));
720 double numeratorLogTargetToUse = backwardPositionsData[0]->logTarget();
721 numContrib = numeratorLogTargetToUse;
722 denContrib = positionsData[0]->logTarget();
723 if ((m_env.subDisplayFile() ) &&
724 (m_env.displayVerbosity() >= 10 ) &&
725 (m_optionsObj->m_totallyMute ==
false)) {
726 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
727 <<
", inputSize = " << inputSize
729 <<
": numContrib = " << numContrib
730 <<
", denContrib = " << denContrib
733 logNumerator += numContrib;
734 logDenominator += denContrib;
736 if ((m_env.subDisplayFile() ) &&
737 (m_env.displayVerbosity() >= 10 ) &&
738 (m_optionsObj->m_totallyMute ==
false)) {
739 *m_env.subDisplayFile() <<
"Leaving MetropolisHastingsSG<P_V,P_M>::alpha(vec)"
740 <<
", inputSize = " << inputSize
741 <<
": alphasNumerator = " << alphasNumerator
742 <<
", alphasDenominator = " << alphasDenominator
743 <<
", logNumerator = " << logNumerator
744 <<
", logDenominator = " << logDenominator
749 return std::min(1.,(alphasNumerator/alphasDenominator)*std::exp(logNumerator-logDenominator));
752 template<
class P_V,
class P_M>
758 if (alpha <= 0. ) result =
false;
759 else if (alpha >= 1. ) result =
true;
760 else if (alpha >= m_env.rngObject()->uniformSample()) result =
true;
766 template<
class P_V,
class P_M>
770 std::ofstream& ofsvar)
const
772 if ((m_env.subDisplayFile() ) &&
773 (m_optionsObj->m_totallyMute ==
false)) {
774 *m_env.subDisplayFile() <<
"\n"
775 <<
"\n-----------------------------------------------------"
776 <<
"\n Writing more information about the Markov chain " << workingChain.
name() <<
" to output file ..."
777 <<
"\n-----------------------------------------------------"
784 if (m_optionsObj->m_rawChainGenerateExtra) {
786 ofsvar << m_optionsObj->m_prefix <<
"logTargets_sub" << m_env.subIdString() <<
" = zeros(" << m_logTargets.size()
790 ofsvar << m_optionsObj->m_prefix <<
"logTargets_sub" << m_env.subIdString() <<
" = [";
791 for (
unsigned int i = 0; i < m_logTargets.size(); ++i) {
792 ofsvar << m_logTargets[i]
798 ofsvar << m_optionsObj->m_prefix <<
"alphaQuotients_sub" << m_env.subIdString() <<
" = zeros(" << m_alphaQuotients.size()
802 ofsvar << m_optionsObj->m_prefix <<
"alphaQuotients_sub" << m_env.subIdString() <<
" = [";
803 for (
unsigned int i = 0; i < m_alphaQuotients.size(); ++i) {
804 ofsvar << m_alphaQuotients[i]
811 ofsvar << m_optionsObj->m_prefix <<
"rejected = " << (double) m_rawChainInfo.numRejections/(
double) (workingChain.
subSequenceSize()-1)
817 ofsvar << m_optionsObj->m_prefix <<
"componentNames = {";
818 m_vectorSpace.printComponentsNames(ofsvar,
false);
822 ofsvar << m_optionsObj->m_prefix <<
"outTargetSupport = " << (double) m_rawChainInfo.numOutOfTargetSupport/(
double) (workingChain.
subSequenceSize()-1)
827 ofsvar << m_optionsObj->m_prefix <<
"runTime = " << m_rawChainInfo.runTime
832 if ((m_env.subDisplayFile() ) &&
833 (m_optionsObj->m_totallyMute ==
false)) {
834 *m_env.subDisplayFile() <<
"\n-----------------------------------------------------"
835 <<
"\n Finished writing more information about the Markov chain " << workingChain.
name()
836 <<
"\n-----------------------------------------------------"
844 template <
class P_V,
class P_M>
858 template <
class P_V,
class P_M>
865 if ((m_env.subDisplayFile() ) &&
866 (m_env.displayVerbosity() >= 5 ) &&
867 (m_optionsObj->m_totallyMute ==
false)) {
868 *m_env.subDisplayFile() <<
"Entering MetropolisHastingsSG<P_V,P_M>::generateSequence()..."
873 std::cerr <<
"'m_vectorSpace' and 'workingChain' are related to vector"
874 <<
"spaces of different dimensions"
880 bool writeLogLikelihood;
881 if ((workingLogLikelihoodValues != NULL) &&
882 (m_optionsObj->m_outputLogLikelihood)) {
883 writeLogLikelihood =
true;
886 writeLogLikelihood =
false;
891 if ((workingLogTargetValues != NULL) &&
892 (m_optionsObj->m_outputLogTarget)) {
893 writeLogTarget =
true;
896 writeLogTarget =
false;
899 MiscCheckTheParallelEnvironment<P_V,P_V>(m_initialPosition,
902 P_V valuesOf1stPosition(m_initialPosition);
905 workingChain.
setName(m_optionsObj->m_prefix +
"rawChain");
911 generateFullChain(valuesOf1stPosition,
912 m_optionsObj->m_rawChainSize,
914 workingLogLikelihoodValues,
915 workingLogTargetValues);
918 readFullChain(m_optionsObj->m_rawChainDataInputFileName,
919 m_optionsObj->m_rawChainDataInputFileType,
920 m_optionsObj->m_rawChainSize,
927 if ((m_env.subDisplayFile() ) &&
928 (m_optionsObj->m_totallyMute ==
false)) {
929 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
930 <<
", prefix = " << m_optionsObj->m_prefix
931 <<
", chain name = " << workingChain.
name()
932 <<
": about to try to open generic output file '" << m_optionsObj->m_dataOutputFileName
934 <<
"', subId = " << m_env.subId()
935 <<
", subenv is allowed to write (1/true or 0/false) = " << (m_optionsObj->m_dataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_dataOutputAllowedSet.end())
941 m_env.openOutputFile(m_optionsObj->m_dataOutputFileName,
943 m_optionsObj->m_dataOutputAllowedSet,
947 if ((m_env.subDisplayFile() ) &&
948 (m_optionsObj->m_totallyMute ==
false)) {
949 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
950 <<
", prefix = " << m_optionsObj->m_prefix
951 <<
", raw chain name = " << workingChain.
name()
952 <<
": returned from opening generic output file '" << m_optionsObj->m_dataOutputFileName
954 <<
"', subId = " << m_env.subId()
964 (m_optionsObj->m_totallyMute ==
false )) {
967 if ((m_env.subDisplayFile() ) &&
968 (m_optionsObj->m_totallyMute ==
false)) {
969 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
970 <<
", prefix = " << m_optionsObj->m_prefix
971 <<
", raw chain name = " << workingChain.
name()
972 <<
": about to try to write raw sub chain output file '" << m_optionsObj->m_rawChainDataOutputFileName
973 <<
"." << m_optionsObj->m_rawChainDataOutputFileType
974 <<
"', subId = " << m_env.subId()
975 <<
", subenv is allowed to write 1/true or 0/false) = " << (m_optionsObj->m_rawChainDataOutputAllowedSet.find(m_env.subId()) != m_optionsObj->m_rawChainDataOutputAllowedSet.end())
980 if ((m_numPositionsNotSubWritten > 0 ) &&
981 (m_optionsObj->m_rawChainDataOutputFileName !=
".")) {
982 workingChain.
subWriteContents(m_optionsObj->m_rawChainSize - m_numPositionsNotSubWritten,
983 m_numPositionsNotSubWritten,
984 m_optionsObj->m_rawChainDataOutputFileName,
985 m_optionsObj->m_rawChainDataOutputFileType,
986 m_optionsObj->m_rawChainDataOutputAllowedSet);
987 if ((m_env.subDisplayFile() ) &&
988 (m_optionsObj->m_totallyMute ==
false)) {
989 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
990 <<
": just wrote (per period request) remaining " << m_numPositionsNotSubWritten <<
" chain positions "
991 <<
", " << m_optionsObj->m_rawChainSize - m_numPositionsNotSubWritten <<
" <= pos <= " << m_optionsObj->m_rawChainSize - 1
995 if (writeLogLikelihood) {
996 workingLogLikelihoodValues->
subWriteContents(m_optionsObj->m_rawChainSize - m_numPositionsNotSubWritten,
997 m_numPositionsNotSubWritten,
998 m_optionsObj->m_rawChainDataOutputFileName +
"_loglikelihood",
999 m_optionsObj->m_rawChainDataOutputFileType,
1000 m_optionsObj->m_rawChainDataOutputAllowedSet);
1003 if (writeLogTarget) {
1004 workingLogTargetValues->
subWriteContents(m_optionsObj->m_rawChainSize - m_numPositionsNotSubWritten,
1005 m_numPositionsNotSubWritten,
1006 m_optionsObj->m_rawChainDataOutputFileName +
"_logtarget",
1007 m_optionsObj->m_rawChainDataOutputFileType,
1008 m_optionsObj->m_rawChainDataOutputAllowedSet);
1011 m_numPositionsNotSubWritten = 0;
1015 if (workingLogLikelihoodValues) {
1018 rawSubMLEpositions);
1021 if ((m_env.subDisplayFile() ) &&
1022 (m_optionsObj->m_totallyMute ==
false)) {
1023 P_V tmpVec(m_vectorSpace.zeroVector());
1025 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1026 <<
": just computed MLE"
1027 <<
", rawSubMLEvalue = " << rawSubMLEvalue
1028 <<
", rawSubMLEpositions.subSequenceSize() = " << rawSubMLEpositions.
subSequenceSize()
1029 <<
", rawSubMLEpositions[0] = " << tmpVec
1035 if (workingLogTargetValues) {
1038 rawSubMAPpositions);
1041 if ((m_env.subDisplayFile() ) &&
1042 (m_optionsObj->m_totallyMute ==
false)) {
1043 P_V tmpVec(m_vectorSpace.zeroVector());
1045 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1046 <<
": just computed MAP"
1047 <<
", rawSubMAPvalue = " << rawSubMAPvalue
1048 <<
", rawSubMAPpositions.subSequenceSize() = " << rawSubMAPpositions.
subSequenceSize()
1049 <<
", rawSubMAPpositions[0] = " << tmpVec
1054 if ((m_env.subDisplayFile() ) &&
1055 (m_optionsObj->m_totallyMute ==
false)) {
1056 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1057 <<
", prefix = " << m_optionsObj->m_prefix
1058 <<
", raw chain name = " << workingChain.
name()
1059 <<
": returned from writing raw sub chain output file '" << m_optionsObj->m_rawChainDataOutputFileName
1060 <<
"." << m_optionsObj->m_rawChainDataOutputFileType
1061 <<
"', subId = " << m_env.subId()
1066 if ((m_env.subDisplayFile() ) &&
1067 (m_optionsObj->m_totallyMute ==
false)) {
1068 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1069 <<
", prefix = " << m_optionsObj->m_prefix
1070 <<
", raw chain name = " << workingChain.
name()
1071 <<
": about to try to write raw unified chain output file '" << m_optionsObj->m_rawChainDataOutputFileName
1072 <<
"." << m_optionsObj->m_rawChainDataOutputFileType
1073 <<
"', subId = " << m_env.subId()
1079 m_optionsObj->m_rawChainDataOutputFileType);
1080 if ((m_env.subDisplayFile() ) &&
1081 (m_optionsObj->m_totallyMute ==
false)) {
1082 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1083 <<
", prefix = " << m_optionsObj->m_prefix
1084 <<
", raw chain name = " << workingChain.
name()
1085 <<
": returned from writing raw unified chain output file '" << m_optionsObj->m_rawChainDataOutputFileName
1086 <<
"." << m_optionsObj->m_rawChainDataOutputFileType
1087 <<
"', subId = " << m_env.subId()
1091 if (writeLogLikelihood) {
1092 workingLogLikelihoodValues->
unifiedWriteContents(m_optionsObj->m_rawChainDataOutputFileName +
"_loglikelihood",
1093 m_optionsObj->m_rawChainDataOutputFileType);
1096 if (writeLogTarget) {
1097 workingLogTargetValues->
unifiedWriteContents(m_optionsObj->m_rawChainDataOutputFileName +
"_logtarget",
1098 m_optionsObj->m_rawChainDataOutputFileType);
1102 if (workingLogLikelihoodValues && (m_env.subRank() == 0)) {
1106 rawUnifiedMLEpositions);
1108 if ((m_env.subDisplayFile() ) &&
1109 (m_optionsObj->m_totallyMute ==
false)) {
1110 P_V tmpVec(m_vectorSpace.zeroVector());
1116 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1117 <<
": just computed MLE"
1118 <<
", rawUnifiedMLEvalue = " << rawUnifiedMLEvalue
1119 <<
", rawUnifiedMLEpositions.subSequenceSize() = " << rawUnifiedMLEpositions.
subSequenceSize()
1120 <<
", rawUnifiedMLEpositions[0] = " << tmpVec
1127 if (workingLogTargetValues && (m_env.subRank() == 0)) {
1130 rawUnifiedMAPpositions);
1132 if ((m_env.subDisplayFile() ) &&
1133 (m_optionsObj->m_totallyMute ==
false)) {
1134 P_V tmpVec(m_vectorSpace.zeroVector());
1140 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1141 <<
": just computed MAP"
1142 <<
", rawUnifiedMAPvalue = " << rawUnifiedMAPvalue
1143 <<
", rawUnifiedMAPpositions.subSequenceSize() = " << rawUnifiedMAPpositions.
subSequenceSize()
1144 <<
", rawUnifiedMAPpositions[0] = " << tmpVec
1152 if ((genericFilePtrSet.
ofsVar ) &&
1153 (m_optionsObj->m_totallyMute ==
false)) {
1155 iRC = writeInfo(workingChain,
1156 *genericFilePtrSet.
ofsVar);
1160 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
1161 if (m_optionsObj->m_rawChainComputeStats) {
1162 workingChain.computeStatistics(*m_optionsObj->m_rawChainStatisticalOptionsObj,
1163 genericFilePtrSet.
ofsVar);
1173 if (m_optionsObj->m_filteredChainGenerate) {
1175 unsigned int filterInitialPos = (
unsigned int) (m_optionsObj->m_filteredChainDiscardedPortion * (
double) workingChain.
subSequenceSize());
1176 unsigned int filterSpacing = m_optionsObj->m_filteredChainLag;
1177 if (filterSpacing == 0) {
1184 workingChain.
filter(filterInitialPos,
1186 workingChain.
setName(m_optionsObj->m_prefix +
"filtChain");
1188 if (workingLogLikelihoodValues) workingLogLikelihoodValues->
filter(filterInitialPos,
1191 if (workingLogTargetValues) workingLogTargetValues->
filter(filterInitialPos,
1195 if ((m_env.subDisplayFile() ) &&
1196 (m_optionsObj->m_totallyMute ==
false)) {
1197 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1198 <<
", prefix = " << m_optionsObj->m_prefix
1199 <<
": checking necessity of opening output files for filtered chain " << workingChain.
name()
1206 (m_optionsObj->m_totallyMute ==
false )) {
1209 m_optionsObj->m_filteredChainDataOutputFileName,
1210 m_optionsObj->m_filteredChainDataOutputFileType,
1211 m_optionsObj->m_filteredChainDataOutputAllowedSet);
1212 if ((m_env.subDisplayFile() ) &&
1213 (m_optionsObj->m_totallyMute ==
false)) {
1214 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1215 <<
", prefix = " << m_optionsObj->m_prefix
1216 <<
": closed sub output file '" << m_optionsObj->m_filteredChainDataOutputFileName
1217 <<
"' for filtered chain " << workingChain.
name()
1221 if (writeLogLikelihood) {
1224 m_optionsObj->m_filteredChainDataOutputFileName +
"_loglikelihood",
1225 m_optionsObj->m_filteredChainDataOutputFileType,
1226 m_optionsObj->m_filteredChainDataOutputAllowedSet);
1229 if (writeLogTarget) {
1232 m_optionsObj->m_filteredChainDataOutputFileName +
"_logtarget",
1233 m_optionsObj->m_filteredChainDataOutputFileType,
1234 m_optionsObj->m_filteredChainDataOutputAllowedSet);
1242 (m_optionsObj->m_totallyMute == false )) {
1244 m_optionsObj->m_filteredChainDataOutputFileType);
1245 if ((m_env.subDisplayFile() ) &&
1246 (m_optionsObj->m_totallyMute ==
false)) {
1247 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1248 <<
", prefix = " << m_optionsObj->m_prefix
1249 <<
": closed unified output file '" << m_optionsObj->m_filteredChainDataOutputFileName
1250 <<
"' for filtered chain " << workingChain.
name()
1254 if (writeLogLikelihood) {
1255 workingLogLikelihoodValues->
unifiedWriteContents(m_optionsObj->m_filteredChainDataOutputFileName +
"_loglikelihood",
1256 m_optionsObj->m_filteredChainDataOutputFileType);
1259 if (writeLogTarget) {
1260 workingLogTargetValues->
unifiedWriteContents(m_optionsObj->m_filteredChainDataOutputFileName +
"_logtarget",
1261 m_optionsObj->m_filteredChainDataOutputFileType);
1268 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
1269 if (m_optionsObj->m_filteredChainComputeStats) {
1270 workingChain.computeStatistics(*m_optionsObj->m_filteredChainStatisticalOptionsObj,
1271 genericFilePtrSet.
ofsVar);
1279 if (genericFilePtrSet.
ofsVar) {
1281 delete genericFilePtrSet.
ofsVar;
1282 if ((m_env.subDisplayFile() ) &&
1283 (m_optionsObj->m_totallyMute ==
false)) {
1284 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1285 <<
", prefix = " << m_optionsObj->m_prefix
1286 <<
": closed generic output file '" << m_optionsObj->m_dataOutputFileName
1287 <<
"' (chain name is " << workingChain.
name()
1293 if ((m_env.subDisplayFile() ) &&
1294 (m_optionsObj->m_totallyMute ==
false)) {
1295 *m_env.subDisplayFile() << std::endl;
1300 if ((m_env.subDisplayFile() ) &&
1301 (m_env.displayVerbosity() >= 5 ) &&
1302 (m_optionsObj->m_totallyMute ==
false)) {
1303 *m_env.subDisplayFile() <<
"Leaving MetropolisHastingsSG<P_V,P_M>::generateSequence()"
1311 template<
class P_V,
class P_M>
1315 info = m_rawChainInfo;
1319 template <
class P_V,
class P_M>
1322 const std::string& inputFileName,
1323 const std::string& inputFileType,
1324 unsigned int chainSize,
1331 template <
class P_V,
class P_M>
1334 const P_V& valuesOf1stPosition,
1335 unsigned int chainSize,
1342 if ((m_env.subDisplayFile() ) &&
1343 (m_optionsObj->m_totallyMute ==
false)) {
1344 *m_env.subDisplayFile() <<
"Starting the generation of Markov chain " << workingChain.
name()
1345 <<
", with " << chainSize
1351 struct timeval timevalChain;
1352 struct timeval timevalCandidate;
1353 struct timeval timevalTarget;
1354 struct timeval timevalMhAlpha;
1356 m_positionIdForDebugging = 0;
1357 m_stageIdForDebugging = 0;
1359 m_rawChainInfo.reset();
1361 iRC = gettimeofday(&timevalChain, NULL);
1364 if ((m_env.subDisplayFile() ) &&
1365 (m_optionsObj->m_totallyMute ==
false)) {
1366 *m_env.subDisplayFile() <<
"\nIn MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1367 <<
": contents of initial position are:";
1368 *m_env.subDisplayFile() << valuesOf1stPosition;
1369 *m_env.subDisplayFile() <<
"\nIn MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1370 <<
": targetPdf.domaintSet() info is:"
1371 << m_targetPdf.domainSet();
1372 *m_env.subDisplayFile() << std::endl;
1376 bool writeLogLikelihood;
1377 if ((workingLogLikelihoodValues != NULL) &&
1378 (m_optionsObj->m_outputLogLikelihood)) {
1379 writeLogLikelihood =
true;
1382 writeLogLikelihood =
false;
1386 bool writeLogTarget;
1387 if ((workingLogTargetValues != NULL) &&
1388 (m_optionsObj->m_outputLogTarget)) {
1389 writeLogTarget =
true;
1392 writeLogTarget =
false;
1395 bool outOfTargetSupport = !m_targetPdf.domainSet().contains(valuesOf1stPosition);
1396 if ((m_env.subDisplayFile()) &&
1397 (outOfTargetSupport )) {
1398 *m_env.subDisplayFile() <<
"ERROR: In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1399 <<
": contents of initial position are:\n";
1400 *m_env.subDisplayFile() << valuesOf1stPosition;
1401 *m_env.subDisplayFile() <<
"\nERROR: In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1402 <<
": targetPdf.domaintSet() info is:\n"
1403 << m_targetPdf.domainSet();
1404 *m_env.subDisplayFile() << std::endl;
1406 queso_require_msg(!(outOfTargetSupport),
"initial position should not be out of target pdf support");
1408 double logPrior = 0.;
1409 double logLikelihood = 0.;
1410 double logTarget = 0.;
1411 if (m_computeInitialPriorAndLikelihoodValues) {
1412 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1413 iRC = gettimeofday(&timevalTarget, NULL);
1416 logTarget = m_targetPdfSynchronizer->callFunction(&valuesOf1stPosition,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood);
1417 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1420 m_rawChainInfo.numTargetCalls++;
1421 if ((m_env.subDisplayFile() ) &&
1422 (m_env.displayVerbosity() >= 3 ) &&
1423 (m_optionsObj->m_totallyMute ==
false)) {
1424 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1425 <<
": just returned from likelihood() for initial chain position"
1426 <<
", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1427 <<
", logPrior = " << logPrior
1428 <<
", logLikelihood = " << logLikelihood
1429 <<
", logTarget = " << logTarget
1434 logPrior = m_initialLogPriorValue;
1435 logLikelihood = m_initialLogLikelihoodValue;
1436 logTarget = logPrior + logLikelihood;
1437 if ((m_env.subDisplayFile() ) &&
1438 (m_env.displayVerbosity() >= 3 ) &&
1439 (m_optionsObj->m_totallyMute ==
false)) {
1440 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1441 <<
": used input prior and likelihood for initial chain position"
1442 <<
", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1443 <<
", logPrior = " << logPrior
1444 <<
", logLikelihood = " << logLikelihood
1445 <<
", logTarget = " << logTarget
1452 valuesOf1stPosition,
1457 P_V gaussianVector(m_vectorSpace.zeroVector());
1458 P_V tmpVecValues(m_vectorSpace.zeroVector());
1465 m_numPositionsNotSubWritten = 0;
1466 if (workingLogLikelihoodValues) workingLogLikelihoodValues->
resizeSequence(chainSize);
1467 if (workingLogTargetValues ) workingLogTargetValues->
resizeSequence (chainSize);
1468 if (
true) m_idsOfUniquePositions.resize(chainSize,0);
1469 if (m_optionsObj->m_rawChainGenerateExtra) {
1470 m_logTargets.resize (chainSize,0.);
1471 m_alphaQuotients.resize(chainSize,0.);
1474 unsigned int uniquePos = 0;
1476 m_numPositionsNotSubWritten++;
1477 if ((m_optionsObj->m_rawChainDataOutputPeriod > 0 ) &&
1478 (((0+1) % m_optionsObj->m_rawChainDataOutputPeriod) == 0 ) &&
1479 (m_optionsObj->m_rawChainDataOutputFileName !=
".")) {
1480 workingChain.
subWriteContents(0 + 1 - m_optionsObj->m_rawChainDataOutputPeriod,
1481 m_optionsObj->m_rawChainDataOutputPeriod,
1482 m_optionsObj->m_rawChainDataOutputFileName,
1483 m_optionsObj->m_rawChainDataOutputFileType,
1484 m_optionsObj->m_rawChainDataOutputAllowedSet);
1485 if ((m_env.subDisplayFile() ) &&
1486 (m_optionsObj->m_totallyMute ==
false)) {
1487 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1488 <<
": just wrote (per period request) " << m_numPositionsNotSubWritten <<
" chain positions "
1489 <<
", " << 0 + 1 - m_optionsObj->m_rawChainDataOutputPeriod <<
" <= pos <= " << 0
1493 if (writeLogLikelihood) {
1494 workingLogLikelihoodValues->
subWriteContents(0 + 1 - m_optionsObj->m_rawChainDataOutputPeriod,
1495 m_optionsObj->m_rawChainDataOutputPeriod,
1496 m_optionsObj->m_rawChainDataOutputFileName +
"_loglikelihood",
1497 m_optionsObj->m_rawChainDataOutputFileType,
1498 m_optionsObj->m_rawChainDataOutputAllowedSet);
1501 if (writeLogTarget) {
1502 workingLogTargetValues->
subWriteContents(0 + 1 - m_optionsObj->m_rawChainDataOutputPeriod,
1503 m_optionsObj->m_rawChainDataOutputPeriod,
1504 m_optionsObj->m_rawChainDataOutputFileName +
"_logtarget",
1505 m_optionsObj->m_rawChainDataOutputFileType,
1506 m_optionsObj->m_rawChainDataOutputAllowedSet);
1509 m_numPositionsNotSubWritten = 0;
1512 if (workingLogLikelihoodValues) (*workingLogLikelihoodValues)[0] = currentPositionData.
logLikelihood();
1513 if (workingLogTargetValues ) (*workingLogTargetValues )[0] = currentPositionData.
logTarget();
1514 if (
true) m_idsOfUniquePositions[uniquePos++] = 0;
1515 if (m_optionsObj->m_rawChainGenerateExtra) {
1516 m_logTargets [0] = currentPositionData.
logTarget();
1517 m_alphaQuotients[0] = 1.;
1521 if ((m_env.subDisplayFile() ) &&
1522 (m_env.displayVerbosity() >= 10 ) &&
1523 (m_optionsObj->m_totallyMute ==
false)) {
1524 *m_env.subDisplayFile() <<
"\n"
1525 <<
"\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
1535 if ((m_env.numSubEnvironments() < (
unsigned int) m_env.fullComm().NumProc()) &&
1536 (m_initialPosition.numOfProcsForStorage() == 1 ) &&
1537 (m_env.subRank() != 0 )) {
1540 aux = m_targetPdfSynchronizer->callFunction(NULL,
1548 for (
unsigned int positionId = 1; positionId < workingChain.
subSequenceSize(); ++positionId) {
1552 m_rawChainInfo.numRejections++;
1555 else for (
unsigned int positionId = 1; positionId < workingChain.
subSequenceSize(); ++positionId) {
1560 m_positionIdForDebugging = positionId;
1561 if ((m_env.subDisplayFile() ) &&
1562 (m_env.displayVerbosity() >= 3 ) &&
1563 (m_optionsObj->m_totallyMute ==
false)) {
1564 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1565 <<
": beginning chain position of id = " << positionId
1566 <<
", m_optionsObj->m_drMaxNumExtraStages = " << m_optionsObj->m_drMaxNumExtraStages
1569 unsigned int stageId = 0;
1570 m_stageIdForDebugging = stageId;
1572 m_tk->clearPreComputingPositions();
1574 if ((m_env.subDisplayFile() ) &&
1575 (m_env.displayVerbosity() >= 5 ) &&
1576 (m_optionsObj->m_totallyMute ==
false)) {
1577 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1578 <<
": about to set TK pre computing position of local id " << 0
1579 <<
", values = " << currentPositionData.
vecValues()
1582 bool validPreComputingPosition = m_tk->setPreComputingPosition(currentPositionData.
vecValues(),0);
1583 if ((m_env.subDisplayFile() ) &&
1584 (m_env.displayVerbosity() >= 5 ) &&
1585 (m_optionsObj->m_totallyMute ==
false)) {
1586 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1587 <<
": returned from setting TK pre computing position of local id " << 0
1588 <<
", values = " << currentPositionData.
vecValues()
1589 <<
", valid = " << validPreComputingPosition
1592 queso_require_msg(validPreComputingPosition,
"initial position should not be an invalid pre computing position");
1598 bool keepGeneratingCandidates =
true;
1599 while (keepGeneratingCandidates) {
1600 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1601 iRC = gettimeofday(&timevalCandidate, NULL);
1605 m_tk->rv(currentPositionData.
vecValues()).realizer().realization(tmpVecValues);
1607 if (m_numDisabledParameters > 0) {
1608 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
1609 if (m_parameterEnabledStatus[paramId] ==
false) {
1610 tmpVecValues[paramId] = m_initialPosition[paramId];
1614 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.candidateRunTime +=
MiscGetEllapsedSeconds(&timevalCandidate);
1616 outOfTargetSupport = !m_targetPdf.domainSet().contains(tmpVecValues);
1618 bool displayDetail = (m_env.displayVerbosity() >= 10) || m_optionsObj->m_displayCandidates;
1619 if ((m_env.subDisplayFile() ) &&
1621 (m_optionsObj->m_totallyMute ==
false)) {
1622 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1623 <<
": for chain position of id = " << positionId
1624 <<
", candidate = " << tmpVecValues
1625 <<
", outOfTargetSupport = " << outOfTargetSupport
1629 if (m_optionsObj->m_putOutOfBoundsInChain) keepGeneratingCandidates =
false;
1630 else keepGeneratingCandidates = outOfTargetSupport;
1633 if ((m_env.subDisplayFile() ) &&
1634 (m_env.displayVerbosity() >= 5 ) &&
1635 (m_optionsObj->m_totallyMute ==
false)) {
1636 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1637 <<
": about to set TK pre computing position of local id " << stageId+1
1638 <<
", values = " << tmpVecValues
1641 validPreComputingPosition = m_tk->setPreComputingPosition(tmpVecValues,stageId+1);
1642 if ((m_env.subDisplayFile() ) &&
1643 (m_env.displayVerbosity() >= 5 ) &&
1644 (m_optionsObj->m_totallyMute ==
false)) {
1645 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1646 <<
": returned from setting TK pre computing position of local id " << stageId+1
1647 <<
", values = " << tmpVecValues
1648 <<
", valid = " << validPreComputingPosition
1652 if (outOfTargetSupport) {
1653 m_rawChainInfo.numOutOfTargetSupport++;
1654 logPrior = -INFINITY;
1655 logLikelihood = -INFINITY;
1656 logTarget = -INFINITY;
1659 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1660 iRC = gettimeofday(&timevalTarget, NULL);
1663 logTarget = m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood);
1664 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime +=
MiscGetEllapsedSeconds(&timevalTarget);
1665 m_rawChainInfo.numTargetCalls++;
1666 if ((m_env.subDisplayFile() ) &&
1667 (m_env.displayVerbosity() >= 3 ) &&
1668 (m_optionsObj->m_totallyMute ==
false)) {
1669 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1670 <<
": just returned from likelihood() for chain position of id " << positionId
1671 <<
", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
1672 <<
", logPrior = " << logPrior
1673 <<
", logLikelihood = " << logLikelihood
1674 <<
", logTarget = " << logTarget
1678 currentCandidateData.set(tmpVecValues,
1683 if ((m_env.subDisplayFile() ) &&
1684 (m_env.displayVerbosity() >= 10 ) &&
1685 (m_optionsObj->m_totallyMute ==
false)) {
1686 *m_env.subDisplayFile() <<
"\n"
1687 <<
"\n-----------------------------------------------------------\n"
1691 bool accept =
false;
1692 double alphaFirstCandidate = 0.;
1693 if (outOfTargetSupport) {
1694 if (m_optionsObj->m_rawChainGenerateExtra) {
1695 m_alphaQuotients[positionId] = 0.;
1699 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1700 iRC = gettimeofday(&timevalMhAlpha, NULL);
1703 if (m_optionsObj->m_rawChainGenerateExtra) {
1704 alphaFirstCandidate = m_algorithm->acceptance_ratio(
1705 currentPositionData,
1706 currentCandidateData,
1711 alphaFirstCandidate = m_algorithm->acceptance_ratio(
1712 currentPositionData,
1713 currentCandidateData,
1717 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.mhAlphaRunTime +=
MiscGetEllapsedSeconds(&timevalMhAlpha);
1718 if ((m_env.subDisplayFile() ) &&
1719 (m_env.displayVerbosity() >= 10 ) &&
1720 (m_optionsObj->m_totallyMute ==
false)) {
1721 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1722 <<
": for chain position of id = " << positionId
1725 accept = acceptAlpha(alphaFirstCandidate);
1728 bool displayDetail = (m_env.displayVerbosity() >= 10) || m_optionsObj->m_displayCandidates;
1729 if ((m_env.subDisplayFile() ) &&
1731 (m_optionsObj->m_totallyMute ==
false)) {
1732 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1733 <<
": for chain position of id = " << positionId
1734 <<
", outOfTargetSupport = " << outOfTargetSupport
1735 <<
", alpha = " << alphaFirstCandidate
1736 <<
", accept = " << accept
1737 <<
", currentCandidateData.vecValues() = ";
1738 *m_env.subDisplayFile() << currentCandidateData.vecValues();
1739 *m_env.subDisplayFile() <<
"\n"
1740 <<
"\n curLogTarget = " << currentPositionData.
logTarget()
1742 <<
"\n canLogTarget = " << currentCandidateData.logTarget()
1746 if ((m_env.subDisplayFile() ) &&
1747 (m_env.displayVerbosity() >= 10 ) &&
1748 (m_optionsObj->m_totallyMute ==
false)) {
1749 *m_env.subDisplayFile() <<
"\n"
1750 <<
"\n-----------------------------------------------------------\n"
1759 if ((accept ==
false) &&
1760 (outOfTargetSupport ==
false) &&
1761 (m_optionsObj->m_drMaxNumExtraStages > 0)) {
1762 accept = delayedRejection(positionId,
1763 currentPositionData,
1764 currentCandidateData);
1773 if (
true) m_idsOfUniquePositions[uniquePos++] = positionId;
1774 currentPositionData = currentCandidateData;
1778 m_rawChainInfo.numRejections++;
1780 m_numPositionsNotSubWritten++;
1781 if ((m_optionsObj->m_rawChainDataOutputPeriod > 0 ) &&
1782 (((positionId+1) % m_optionsObj->m_rawChainDataOutputPeriod) == 0 ) &&
1783 (m_optionsObj->m_rawChainDataOutputFileName !=
".")) {
1784 if ((m_env.subDisplayFile() ) &&
1785 (m_env.displayVerbosity() >= 10 ) &&
1786 (m_optionsObj->m_totallyMute ==
false)) {
1787 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1788 <<
", for chain position of id = " << positionId
1789 <<
": about to write (per period request) " << m_numPositionsNotSubWritten <<
" chain positions "
1790 <<
", " << positionId + 1 - m_optionsObj->m_rawChainDataOutputPeriod <<
" <= pos <= " << positionId
1793 workingChain.
subWriteContents(positionId + 1 - m_optionsObj->m_rawChainDataOutputPeriod,
1794 m_optionsObj->m_rawChainDataOutputPeriod,
1795 m_optionsObj->m_rawChainDataOutputFileName,
1796 m_optionsObj->m_rawChainDataOutputFileType,
1797 m_optionsObj->m_rawChainDataOutputAllowedSet);
1798 if ((m_env.subDisplayFile() ) &&
1799 (m_optionsObj->m_totallyMute ==
false)) {
1800 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1801 <<
", for chain position of id = " << positionId
1802 <<
": just wrote (per period request) " << m_numPositionsNotSubWritten <<
" chain positions "
1803 <<
", " << positionId + 1 - m_optionsObj->m_rawChainDataOutputPeriod <<
" <= pos <= " << positionId
1807 if (writeLogLikelihood) {
1808 workingLogLikelihoodValues->
subWriteContents(0 + 1 - m_optionsObj->m_rawChainDataOutputPeriod,
1809 m_optionsObj->m_rawChainDataOutputPeriod,
1810 m_optionsObj->m_rawChainDataOutputFileName +
"_loglikelihood",
1811 m_optionsObj->m_rawChainDataOutputFileType,
1812 m_optionsObj->m_rawChainDataOutputAllowedSet);
1815 if (writeLogTarget) {
1816 workingLogTargetValues->
subWriteContents(0 + 1 - m_optionsObj->m_rawChainDataOutputPeriod,
1817 m_optionsObj->m_rawChainDataOutputPeriod,
1818 m_optionsObj->m_rawChainDataOutputFileName +
"_logtarget",
1819 m_optionsObj->m_rawChainDataOutputFileType,
1820 m_optionsObj->m_rawChainDataOutputAllowedSet);
1823 m_numPositionsNotSubWritten = 0;
1827 if (workingLogLikelihoodValues) (*workingLogLikelihoodValues)[positionId] = currentPositionData.
logLikelihood();
1828 if (workingLogTargetValues ) (*workingLogTargetValues )[positionId] = currentPositionData.
logTarget();
1830 if (m_optionsObj->m_rawChainGenerateExtra) {
1831 m_logTargets[positionId] = currentPositionData.
logTarget();
1834 if (m_optionsObj->m_enableBrooksGelmanConvMonitor > 0) {
1835 if (positionId % m_optionsObj->m_enableBrooksGelmanConvMonitor == 0 &&
1836 positionId > m_optionsObj->m_BrooksGelmanLag + 1) {
1839 m_optionsObj->m_BrooksGelmanLag,
1840 positionId - m_optionsObj->m_BrooksGelmanLag);
1842 if (m_env.subDisplayFile()) {
1843 *m_env.subDisplayFile() <<
"positionId = " << positionId
1844 <<
", conv_est = " << conv_est
1846 (*m_env.subDisplayFile()).flush();
1855 this->adapt(positionId, workingChain);
1861 if ((m_env.subDisplayFile() ) &&
1862 (m_env.displayVerbosity() >= 3 ) &&
1863 (m_optionsObj->m_totallyMute ==
false)) {
1864 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
1865 <<
": finishing chain position of id = " << positionId
1866 <<
", accept = " << accept
1867 <<
", curLogTarget = " << currentPositionData.
logTarget()
1868 <<
", canLogTarget = " << currentCandidateData.logTarget()
1872 if ((m_optionsObj->m_rawChainDisplayPeriod > 0) &&
1873 (((positionId+1) % m_optionsObj->m_rawChainDisplayPeriod) == 0)) {
1874 if ((m_env.subDisplayFile() ) &&
1875 (m_optionsObj->m_totallyMute ==
false)) {
1876 *m_env.subDisplayFile() <<
"Finished generating " << positionId+1
1878 <<
", current rejection percentage = " << (100. * ((double) m_rawChainInfo.numRejections)/((double) (positionId+1)))
1884 if ((m_env.subDisplayFile() ) &&
1885 (m_env.displayVerbosity() >= 10 ) &&
1886 (m_optionsObj->m_totallyMute ==
false)) {
1887 *m_env.subDisplayFile() <<
"\n"
1888 <<
"\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
1894 if ((m_env.numSubEnvironments() < (
unsigned int) m_env.fullComm().NumProc()) &&
1895 (m_initialPosition.numOfProcsForStorage() == 1 ) &&
1896 (m_env.subRank() == 0 )) {
1899 aux = m_targetPdfSynchronizer->callFunction(NULL,
1913 if ((m_env.subDisplayFile() ) &&
1914 (m_optionsObj->m_totallyMute ==
false)) {
1915 *m_env.subDisplayFile() <<
"Finished the generation of Markov chain " << workingChain.
name()
1918 *m_env.subDisplayFile() <<
"\nSome information about this chain:"
1919 <<
"\n Chain run time = " << m_rawChainInfo.runTime
1921 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1922 *m_env.subDisplayFile() <<
"\n\n Breaking of the chain run time:\n";
1923 *m_env.subDisplayFile() <<
"\n Candidate run time = " << m_rawChainInfo.candidateRunTime
1924 <<
" seconds (" << 100.*m_rawChainInfo.candidateRunTime/m_rawChainInfo.runTime
1926 *m_env.subDisplayFile() <<
"\n Num target calls = " << m_rawChainInfo.numTargetCalls;
1927 *m_env.subDisplayFile() <<
"\n Target d. run time = " << m_rawChainInfo.targetRunTime
1928 <<
" seconds (" << 100.*m_rawChainInfo.targetRunTime/m_rawChainInfo.runTime
1930 *m_env.subDisplayFile() <<
"\n Avg target run time = " << m_rawChainInfo.targetRunTime/((double) m_rawChainInfo.numTargetCalls)
1932 *m_env.subDisplayFile() <<
"\n Mh alpha run time = " << m_rawChainInfo.mhAlphaRunTime
1933 <<
" seconds (" << 100.*m_rawChainInfo.mhAlphaRunTime/m_rawChainInfo.runTime
1935 *m_env.subDisplayFile() <<
"\n Dr alpha run time = " << m_rawChainInfo.drAlphaRunTime
1936 <<
" seconds (" << 100.*m_rawChainInfo.drAlphaRunTime/m_rawChainInfo.runTime
1938 *m_env.subDisplayFile() <<
"\n---------------------- --------------";
1939 double sumRunTime = m_rawChainInfo.candidateRunTime + m_rawChainInfo.targetRunTime + m_rawChainInfo.mhAlphaRunTime + m_rawChainInfo.drAlphaRunTime;
1940 *m_env.subDisplayFile() <<
"\n Sum = " << sumRunTime
1941 <<
" seconds (" << 100.*sumRunTime/m_rawChainInfo.runTime
1943 *m_env.subDisplayFile() <<
"\n\n Other run times:";
1944 *m_env.subDisplayFile() <<
"\n DR run time = " << m_rawChainInfo.drRunTime
1945 <<
" seconds (" << 100.*m_rawChainInfo.drRunTime/m_rawChainInfo.runTime
1947 *m_env.subDisplayFile() <<
"\n AM run time = " << m_rawChainInfo.amRunTime
1948 <<
" seconds (" << 100.*m_rawChainInfo.amRunTime/m_rawChainInfo.runTime
1951 *m_env.subDisplayFile() <<
"\n Number of DRs = " << m_rawChainInfo.numDRs <<
"(num_DRs/chain_size = " << (double) m_rawChainInfo.numDRs/(
double) workingChain.
subSequenceSize()
1953 *m_env.subDisplayFile() <<
"\n Out of target support in DR = " << m_rawChainInfo.numOutOfTargetSupportInDR;
1954 *m_env.subDisplayFile() <<
"\n Rejection percentage = " << 100. * (double) m_rawChainInfo.numRejections/(
double) workingChain.
subSequenceSize()
1956 *m_env.subDisplayFile() <<
"\n Out of target support percentage = " << 100. * (double) m_rawChainInfo.numOutOfTargetSupport/(
double) workingChain.
subSequenceSize()
1958 *m_env.subDisplayFile() << std::endl;
1969 template <
class P_V,
class P_M>
1975 struct timeval timevalAM;
1978 if ((m_optionsObj->m_tkUseLocalHessian ==
true) ||
1979 (m_optionsObj->m_amInitialNonAdaptInterval == 0) ||
1980 (m_optionsObj->m_amAdaptInterval == 0)) {
1985 if (m_optionsObj->m_rawChainMeasureRunTimes) {
1986 iRC = gettimeofday(&timevalAM, NULL);
1990 unsigned int idOfFirstPositionInSubChain = 0;
1994 bool printAdaptedMatrix =
false;
1995 if (positionId < m_optionsObj->m_amInitialNonAdaptInterval) {
1998 else if (positionId == m_optionsObj->m_amInitialNonAdaptInterval) {
1999 idOfFirstPositionInSubChain = 0;
2000 partialChain.
resizeSequence(m_optionsObj->m_amInitialNonAdaptInterval+1);
2001 m_lastMean = m_vectorSpace.newVector();
2002 m_lastAdaptedCovMatrix = m_vectorSpace.newMatrix();
2003 printAdaptedMatrix =
true;
2006 unsigned int interval = positionId - m_optionsObj->m_amInitialNonAdaptInterval;
2007 if ((interval % m_optionsObj->m_amAdaptInterval) == 0) {
2008 idOfFirstPositionInSubChain = positionId - m_optionsObj->m_amAdaptInterval;
2011 if (m_optionsObj->m_amAdaptedMatricesDataOutputPeriod > 0) {
2012 if ((interval % m_optionsObj->m_amAdaptedMatricesDataOutputPeriod) == 0) {
2013 printAdaptedMatrix =
true;
2022 if (m_optionsObj->m_rawChainMeasureRunTimes) {
2030 P_V transporterVec(m_vectorSpace.zeroVector());
2036 if (this->m_optionsObj->m_algorithm ==
"logit_random_walk") {
2040 P_V transformedTransporterVec(m_vectorSpace.zeroVector());
2042 m_tk.get())->transformToGaussianSpace(transporterVec,
2043 transformedTransporterVec);
2051 updateAdaptedCovMatrix(partialChain,
2052 idOfFirstPositionInSubChain,
2055 *m_lastAdaptedCovMatrix);
2058 if ((printAdaptedMatrix ==
true) &&
2059 (m_optionsObj->m_amAdaptedMatricesDataOutputFileName !=
"." )) {
2060 char varNamePrefix[64];
2061 sprintf(varNamePrefix,
"mat_am%d",positionId);
2064 sprintf(tmpChar,
"_am%d",positionId);
2066 std::set<unsigned int> tmpSet;
2067 tmpSet.insert(m_env.subId());
2069 m_lastAdaptedCovMatrix->subWriteContents(varNamePrefix,
2070 (m_optionsObj->m_amAdaptedMatricesDataOutputFileName+tmpChar),
2071 m_optionsObj->m_amAdaptedMatricesDataOutputFileType,
2073 if ((m_env.subDisplayFile() ) &&
2074 (m_optionsObj->m_totallyMute ==
false)) {
2075 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2076 <<
": just wrote last adapted proposal cov matrix contents = " << *m_lastAdaptedCovMatrix
2082 bool tmpCholIsPositiveDefinite =
false;
2083 P_M tmpChol(*m_lastAdaptedCovMatrix);
2084 P_M attemptedMatrix(tmpChol);
2085 if ((m_env.subDisplayFile() ) &&
2086 (m_env.displayVerbosity() >= 10)) {
2088 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2089 <<
", positionId = " << positionId
2090 <<
": 'am' calling first tmpChol.chol()"
2093 iRC = tmpChol.chol();
2095 std::string err1 =
"In MetropolisHastingsSG<P_V,P_M>::adapt(): first ";
2096 err1 +=
"Cholesky factorisation of proposal covariance matrix ";
2097 err1 +=
"failed. QUESO will now attempt to regularise the ";
2098 err1 +=
"matrix before proceeding. This is not a fatal error.";
2099 std::cerr << err1 << std::endl;
2103 if ((m_env.subDisplayFile() ) &&
2104 (m_env.displayVerbosity() >= 10)) {
2106 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2107 <<
", positionId = " << positionId
2108 <<
": 'am' got first tmpChol.chol() with iRC = " << iRC
2111 double diagMult = 1.;
2112 for (
unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2113 diagMult *= tmpChol(j,j);
2115 *m_env.subDisplayFile() <<
"diagMult = " << diagMult
2120 #if 0 // tentative logic
2122 double diagMult = 1.;
2123 for (
unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2124 diagMult *= tmpChol(j,j);
2126 if (diagMult < 1.e-40) {
2137 P_M* tmpDiag = m_vectorSpace.newDiagMatrix(m_optionsObj->m_amEpsilon);
2138 if (m_numDisabledParameters > 0) {
2139 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2140 if (m_parameterEnabledStatus[paramId] ==
false) {
2141 (*tmpDiag)(paramId,paramId) = 0.;
2145 tmpChol = *m_lastAdaptedCovMatrix + *tmpDiag;
2146 attemptedMatrix = tmpChol;
2149 if ((m_env.subDisplayFile() ) &&
2150 (m_env.displayVerbosity() >= 10)) {
2151 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2152 <<
", positionId = " << positionId
2153 <<
": 'am' calling second tmpChol.chol()"
2158 iRC = tmpChol.chol();
2160 std::string err2 =
"In MetropolisHastingsSG<P_V,P_M>::adapt(): second ";
2161 err2 +=
"Cholesky factorisation of (regularised) proposal ";
2162 err2 +=
"covariance matrix failed. QUESO is falling back to ";
2163 err2 +=
"the last factorisable proposal covariance matrix. ";
2164 err2 +=
"This is not a fatal error.";
2165 std::cerr << err2 << std::endl;
2169 if ((m_env.subDisplayFile() ) &&
2170 (m_env.displayVerbosity() >= 10)) {
2171 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2172 <<
", positionId = " << positionId
2173 <<
": 'am' got second tmpChol.chol() with iRC = " << iRC
2176 double diagMult = 1.;
2177 for (
unsigned int j = 0; j < tmpChol.numRowsLocal(); ++j) {
2178 diagMult *= tmpChol(j,j);
2180 *m_env.subDisplayFile() <<
"diagMult = " << diagMult
2184 *m_env.subDisplayFile() <<
"attemptedMatrix = " << attemptedMatrix
2195 tmpCholIsPositiveDefinite =
true;
2199 tmpCholIsPositiveDefinite =
true;
2203 if (tmpCholIsPositiveDefinite) {
2204 P_M tmpMatrix(m_optionsObj->m_amEta*attemptedMatrix);
2205 if (m_numDisabledParameters > 0) {
2206 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2207 if (m_parameterEnabledStatus[paramId] ==
false) {
2208 for (
unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2209 tmpMatrix(i,paramId) = 0.;
2211 for (
unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2212 tmpMatrix(paramId,j) = 0.;
2214 tmpMatrix(paramId,paramId) = 1.;
2221 if (this->m_optionsObj->m_algorithm ==
"logit_random_walk") {
2223 ->updateLawCovMatrix(tmpMatrix);
2225 else if (this->m_optionsObj->m_algorithm ==
"random_walk") {
2227 ->updateLawCovMatrix(tmpMatrix);
2230 #ifdef UQ_DRAM_MCG_REQUIRES_INVERTED_COV_MATRICES
2241 if (m_optionsObj->m_rawChainMeasureRunTimes) {
2246 template <
class P_V,
class P_M>
2252 if ((m_optionsObj->m_drDuringAmNonAdaptiveInt ==
false ) &&
2253 (m_optionsObj->m_tkUseLocalHessian ==
false ) &&
2254 (m_optionsObj->m_amInitialNonAdaptInterval > 0 ) &&
2255 (m_optionsObj->m_amAdaptInterval > 0 ) &&
2256 (positionId <= m_optionsObj->m_amInitialNonAdaptInterval)) {
2260 unsigned int stageId = 0;
2262 bool validPreComputingPosition;
2264 m_tk->clearPreComputingPositions();
2266 validPreComputingPosition = m_tk->setPreComputingPosition(
2269 validPreComputingPosition = m_tk->setPreComputingPosition(
2270 currentCandidateData.
vecValues(), stageId + 1);
2272 std::vector<MarkovChainPositionData<P_V>*> drPositionsData(stageId+2,NULL);
2273 std::vector<unsigned int> tkStageIds (stageId+2,0);
2276 struct timeval timevalDR;
2277 struct timeval timevalDrAlpha;
2278 struct timeval timevalCandidate;
2279 struct timeval timevalTarget;
2281 if (m_optionsObj->m_rawChainMeasureRunTimes) {
2282 iRC = gettimeofday(&timevalDR, NULL);
2292 bool accept =
false;
2293 while ((validPreComputingPosition ==
true ) &&
2294 (accept == false ) &&
2295 (stageId < m_optionsObj->m_drMaxNumExtraStages)) {
2296 if ((m_env.subDisplayFile() ) &&
2297 (m_env.displayVerbosity() >= 10 ) &&
2298 (m_optionsObj->m_totallyMute ==
false)) {
2299 *m_env.subDisplayFile() <<
"\n"
2300 <<
"\n+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"
2304 m_rawChainInfo.numDRs++;
2306 m_stageIdForDebugging = stageId;
2307 if ((m_env.subDisplayFile() ) &&
2308 (m_env.displayVerbosity() >= 10 ) &&
2309 (m_optionsObj->m_totallyMute ==
false)) {
2310 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2311 <<
": for chain position of id = " << positionId
2312 <<
", beginning stageId = " << stageId
2316 P_V tmpVecValues(currentCandidateData.
vecValues());
2317 bool keepGeneratingCandidates =
true;
2318 bool outOfTargetSupport =
false;
2319 while (keepGeneratingCandidates) {
2320 if (m_optionsObj->m_rawChainMeasureRunTimes) {
2321 iRC = gettimeofday(&timevalCandidate, NULL);
2324 m_tk->rv(tkStageIds).realizer().realization(tmpVecValues);
2325 if (m_numDisabledParameters > 0) {
2326 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2327 if (m_parameterEnabledStatus[paramId] ==
false) {
2328 tmpVecValues[paramId] = m_initialPosition[paramId];
2332 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.candidateRunTime +=
MiscGetEllapsedSeconds(&timevalCandidate);
2334 outOfTargetSupport = !m_targetPdf.domainSet().contains(tmpVecValues);
2336 if (m_optionsObj->m_putOutOfBoundsInChain) keepGeneratingCandidates =
false;
2337 else keepGeneratingCandidates = outOfTargetSupport;
2340 if ((m_env.subDisplayFile() ) &&
2341 (m_env.displayVerbosity() >= 5 ) &&
2342 (m_optionsObj->m_totallyMute ==
false)) {
2343 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2344 <<
": about to set TK pre computing position of local id " << stageId+1
2345 <<
", values = " << tmpVecValues
2348 validPreComputingPosition = m_tk->setPreComputingPosition(tmpVecValues,stageId+1);
2349 if ((m_env.subDisplayFile() ) &&
2350 (m_env.displayVerbosity() >= 5 ) &&
2351 (m_optionsObj->m_totallyMute ==
false)) {
2352 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2353 <<
": returned from setting TK pre computing position of local id " << stageId+1
2354 <<
", values = " << tmpVecValues
2355 <<
", valid = " << validPreComputingPosition
2360 double logLikelihood;
2362 if (outOfTargetSupport) {
2363 m_rawChainInfo.numOutOfTargetSupportInDR++;
2364 logPrior = -INFINITY;
2365 logLikelihood = -INFINITY;
2366 logTarget = -INFINITY;
2369 if (m_optionsObj->m_rawChainMeasureRunTimes) {
2370 iRC = gettimeofday(&timevalTarget, NULL);
2373 logTarget = m_targetPdfSynchronizer->callFunction(&tmpVecValues,NULL,NULL,NULL,NULL,&logPrior,&logLikelihood);
2374 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.targetRunTime +=
MiscGetEllapsedSeconds(&timevalTarget);
2375 m_rawChainInfo.numTargetCalls++;
2376 if ((m_env.subDisplayFile() ) &&
2377 (m_env.displayVerbosity() >= 3 ) &&
2378 (m_optionsObj->m_totallyMute ==
false)) {
2379 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2380 <<
": just returned from likelihood() for chain position of id " << positionId
2381 <<
", m_rawChainInfo.numTargetCalls = " << m_rawChainInfo.numTargetCalls
2382 <<
", stageId = " << stageId
2383 <<
", logPrior = " << logPrior
2384 <<
", logLikelihood = " << logLikelihood
2385 <<
", logTarget = " << logTarget
2389 currentCandidateData.
set(tmpVecValues,
2397 tkStageIds.push_back (stageId+1);
2399 double alphaDR = 0.;
2400 if (outOfTargetSupport ==
false) {
2401 if (m_optionsObj->m_rawChainMeasureRunTimes) {
2402 iRC = gettimeofday(&timevalDrAlpha, NULL);
2405 alphaDR = this->alpha(drPositionsData,tkStageIds);
2406 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.drAlphaRunTime +=
MiscGetEllapsedSeconds(&timevalDrAlpha);
2407 accept = acceptAlpha(alphaDR);
2410 bool displayDetail = (m_env.displayVerbosity() >= 10) || m_optionsObj->m_displayCandidates;
2411 if ((m_env.subDisplayFile() ) &&
2413 (m_optionsObj->m_totallyMute ==
false)) {
2414 *m_env.subDisplayFile() <<
"In MetropolisHastingsSG<P_V,P_M>::generateFullChain()"
2415 <<
": for chain position of id = " << positionId
2416 <<
" and stageId = " << stageId
2417 <<
", outOfTargetSupport = " << outOfTargetSupport
2418 <<
", alpha = " << alphaDR
2419 <<
", accept = " << accept
2420 <<
", currentCandidateData.vecValues() = ";
2421 *m_env.subDisplayFile() << currentCandidateData.
vecValues();
2422 *m_env.subDisplayFile() << std::endl;
2426 if (m_optionsObj->m_rawChainMeasureRunTimes) m_rawChainInfo.drRunTime +=
MiscGetEllapsedSeconds(&timevalDR);
2428 for (
unsigned int i = 0; i < drPositionsData.size(); ++i) {
2429 if (drPositionsData[i])
delete drPositionsData[i];
2436 template <
class P_V,
class P_M>
2440 unsigned int idOfFirstPositionInSubChain,
2441 double& lastChainSize,
2443 P_M& lastAdaptedCovMatrix)
2446 if (lastChainSize == 0) {
2449 #if 1 // prudenci-2012-07-06
2455 P_V tmpVec(m_vectorSpace.zeroVector());
2456 lastAdaptedCovMatrix = -doubleSubChainSize *
matrixProduct(lastMean,lastMean);
2461 lastAdaptedCovMatrix /= (doubleSubChainSize - 1.);
2467 P_V tmpVec (m_vectorSpace.zeroVector());
2468 P_V diffVec(m_vectorSpace.zeroVector());
2470 double doubleCurrentId = (double) (idOfFirstPositionInSubChain+i);
2472 diffVec = tmpVec - lastMean;
2474 double ratio1 = (1. - 1./doubleCurrentId);
2475 double ratio2 = (1./(1.+doubleCurrentId));
2476 lastAdaptedCovMatrix = ratio1 * lastAdaptedCovMatrix + ratio2 *
matrixProduct(diffVec,diffVec);
2477 lastMean += ratio2 * diffVec;
2480 lastChainSize += doubleSubChainSize;
2482 if (m_numDisabledParameters > 0) {
2483 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2484 if (m_parameterEnabledStatus[paramId] ==
false) {
2485 for (
unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2486 lastAdaptedCovMatrix(i,paramId) = 0.;
2488 for (
unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2489 lastAdaptedCovMatrix(paramId,j) = 0.;
2491 lastAdaptedCovMatrix(paramId,paramId) = 1.;
2499 template <
class P_V,
class P_M>
2504 P_V min_domain_bounds(boxSubset.
minValues());
2505 P_V max_domain_bounds(boxSubset.
maxValues());
2507 for (
unsigned int i = 0; i < min_domain_bounds.sizeLocal(); i++) {
2508 double min_val = min_domain_bounds[i];
2509 double max_val = max_domain_bounds[i];
2512 if (m_initialProposalCovMatrix(i, i) >= max_val - min_val) {
2515 std::cerr <<
"Proposal variance element "
2518 << m_initialProposalCovMatrix(i, i)
2519 <<
" but domain is of size "
2520 << max_val - min_val
2522 std::cerr <<
"QUESO does not support uniform-like proposal "
2523 <<
"distributions. Try making the proposal variance smaller"
2528 double transformJacobian = 4.0 / (max_val - min_val);
2532 for (
unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
2534 m_initialProposalCovMatrix(j, i) *= transformJacobian;
2536 for (
unsigned int j = 0; j < min_domain_bounds.sizeLocal(); j++) {
2538 m_initialProposalCovMatrix(i, j) *= transformJacobian;
bool m_totallyMute
If true, zero output is written to files. Default is false.
MhOptionsValues m_ov
This class is where the actual options are stored.
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
double logLikelihood() const
Logarithm of the value of the likelihood; access to private attribute m_logLikelihood.
static void set_dr_scales(const std::vector< double > &scales)
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
#define queso_require_greater_equal_msg(expr1, expr2, msg)
bool acceptAlpha(double alpha)
Decides whether or not to accept alpha.
#define queso_warning(message)
static void set_tk(const BaseTKGroup< GslVector, GslMatrix > &tk)
bool delayedRejection(unsigned int positionId, const MarkovChainPositionData< P_V > ¤tPositionData, MarkovChainPositionData< P_V > ¤tCandidateData)
Does delayed rejection.
void commonConstructor()
Reads the options values from the options input file.
unsigned int numOutOfTargetSupport
This class reads the options for the Metropolis-Hastings generator of samples from an input file...
~MetropolisHastingsSG()
Destructor.
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
MHRawChainInfoStruct & operator=(const MHRawChainInfoStruct &rhs)
Assignment operator.
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
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.
static void set_options(const MhOptionsValues &options)
void copy(const MHRawChainInfoStruct &src)
Copies Metropolis-Hastings chain info from src to this.
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.
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
static void set_target_pdf(const BaseJointPdf< GslVector, GslMatrix > &target_pdf)
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
unsigned int vectorSizeLocal() const
Local dimension (size) of the vector space.
#define queso_require_equal_to_msg(expr1, expr2, msg)
MHRawChainInfoStruct & operator+=(const MHRawChainInfoStruct &rhs)
Addition assignment operator.
unsigned int numRejections
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.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
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...
const V & subMeanPlain() const
Finds the mean value of the sub-sequence.
void adapt(unsigned int positionId, BaseVectorSequence< P_V, P_M > &workingChain)
Adaptive Metropolis method that deals with adapting the proposal covariance matrix.
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.
MHRawChainInfoStruct()
Constructor.
Struct for handling data input and output from files.
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
const int UQ_INCOMPLETE_IMPLEMENTATION_RC
bool queso_isfinite(T arg)
A templated class that represents a Metropolis-Hastings generator of samples.
virtual void subMeanExtra(unsigned int initialPos, unsigned int numPos, V &meanVec) const =0
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
void set(const V &vecValues, bool outOfTargetSupport, double logLikelihood, double logTarget)
Sets the values of the chain.
MetropolisHastingsSGOptions * m_oldOptions
#define UQ_MH_SG_DO_LOGIT_TRANSFORM
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.
void reset()
Resets Metropolis-Hastings chain info.
const MhOptionsValues * m_optionsObj
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.
This class provides options for each level of the Multilevel sequence generator if no input file is a...
void mpiSum(const MpiComm &comm, MHRawChainInfoStruct &sumInfo)
Calculates the MPI sum of this.
unsigned int numOutOfTargetSupportInDR
void transformInitialCovMatrixToGaussianSpace(const BoxSubset< P_V, P_M > &boxSubset)
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
A struct that represents a Metropolis-Hastings sample.
static void set_pdf_synchronizer(const ScalarFunctionSynchronizer< GslVector, GslMatrix > &synchronizer)
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
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...
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.
const V & vecValues() const
Values of the chain (vector); access to private attribute m_vecValues.
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
MetropolisHastingsSG(const char *prefix, const MhOptionsValues *alternativeOptionsValues, const BaseVectorRV< P_V, P_M > &sourceRv, const P_V &initialPosition, const P_M *inputProposalCovMatrix)
Constructor.
double MiscGetEllapsedSeconds(struct timeval *timeval0)
static void set_vectorspace(const VectorSpace< GslVector, GslMatrix > &v)
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define queso_require_msg(asserted, msg)
The QUESO MPI Communicator Class.
P_M m_initialProposalCovMatrix
This base class allows the representation of a transition kernel.
double logTarget() const
Logarithm of the value of the target; access to private attribute m_logTarget.
static SharedPtr< BaseTKGroup< GslVector, GslMatrix > >::Type build(const std::string &name)
static void set_environment(const BaseEnvironment &env)
This class allows the representation of a transition kernel with a scaled covariance matrix...
double alpha(const std::vector< MarkovChainPositionData< P_V > * > &inputPositions, const std::vector< unsigned int > &inputTKStageIds)
Calculates acceptance ratio.
This class represents a transition kernel with a scaled covariance matrix on hybrid bounded/unbounded...
const V & minValues() const
Vector of the minimum values of the box subset.
unsigned int numTargetCalls
This class provides options for the Metropolis-Hastings generator of samples if no input file is avai...
~MHRawChainInfoStruct()
Destructor.
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
double unifiedPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
static void set_initial_cov_matrix(GslMatrix &cov_matrix)
const BaseEnvironment & m_env
const std::string & name() const
Access to protected attribute m_name: name of the sequence of vectors.
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
double subPositionsOfMaximum(const ScalarSequence< double > &subCorrespondingScalarValues, BaseVectorSequence< V, M > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
virtual double estimateConvBrooksGelman(unsigned int initialPos, unsigned int numPos) const =0
Estimates convergence rate using Brooks & Gelman method. See template specialization.
const V & maxValues() const
Vector of the maximum values of the box subset.
int writeInfo(const BaseVectorSequence< P_V, P_M > &workingChain, std::ofstream &ofsvar) const
Writes information about the Markov chain in a file.
const BaseTKGroup< P_V, P_M > & transitionKernel() const
Returns the underlying transition kernel for this sequence generator.
A templated class that represents a Markov Chain.