25 #include <queso/MLSampling.h>
26 #include <queso/InstantiateIntersection.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
34 void BIP_routine(glp_tree *tree,
void *info)
36 const BaseEnvironment& env = *(((BIP_routine_struct *) info)->env);
37 unsigned int currLevel = ((BIP_routine_struct *) info)->currLevel;
39 int reason = glp_ios_reason(tree);
41 if ((env.subDisplayFile()) && (env.displayVerbosity() >= 1)) {
42 *env.subDisplayFile() <<
"In BIP_routine()"
44 <<
": glp_ios_reason() = " << reason
47 std::cout <<
"In BIP_routine: reason = " << reason << std::endl;
82 "invalid glp_ios_readon");
89 #endif // QUESO_HAS_GLPK
91 template <
class P_V,
class P_M>
94 unsigned int unifiedRequestedNumSamples,
95 const std::vector<double>& unifiedWeightStdVectorAtProc0Only,
96 std::vector<unsigned int>& unifiedIndexCountersAtProc0Only)
98 if (m_env.inter0Rank() != 0)
return;
100 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
101 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::sampleIndexes_proc0()"
103 <<
", step " << m_currStep
104 <<
": unifiedRequestedNumSamples = " << unifiedRequestedNumSamples
105 <<
", unifiedWeightStdVectorAtProc0Only.size() = " << unifiedWeightStdVectorAtProc0Only.size()
109 #if 0 // For debug only
110 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
111 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::sampleIndexes_proc0()"
113 <<
", step " << m_currStep
116 unsigned int numZeros = 0;
117 for (
unsigned int i = 0; i < unifiedWeightStdVectorAtProc0Only.size(); ++i) {
118 *m_env.subDisplayFile() <<
" unifiedWeightStdVectorAtProc0Only[" << i
119 <<
"] = " << unifiedWeightStdVectorAtProc0Only[i]
121 if (unifiedWeightStdVectorAtProc0Only[i] == 0.) numZeros++;
123 *m_env.subDisplayFile() <<
"Number of zeros in unifiedWeightStdVectorAtProc0Only = " << numZeros
128 if (m_env.inter0Rank() == 0) {
129 unsigned int resizeSize = unifiedWeightStdVectorAtProc0Only.size();
130 unifiedIndexCountersAtProc0Only.resize(resizeSize,0);
135 unifiedWeightStdVectorAtProc0Only);
136 for (
unsigned int i = 0; i < unifiedRequestedNumSamples; ++i) {
137 unsigned int index = tmpFd.
sample();
138 unifiedIndexCountersAtProc0Only[index] += 1;
145 template <
class P_V,
class P_M>
149 unsigned int indexOfFirstWeight,
150 unsigned int indexOfLastWeight,
151 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
152 std::vector<ExchangeInfoStruct>& exchangeStdVec)
156 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
157 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
159 <<
", step " << m_currStep
160 <<
": indexOfFirstWeight = " << indexOfFirstWeight
161 <<
", indexOfLastWeight = " << indexOfLastWeight
167 if (m_env.inter0Rank() >= 0) {
168 Np = (
unsigned int) m_env.inter0Comm().NumProc();
170 std::vector<unsigned int> allFirstIndexes(Np,0);
171 std::vector<unsigned int> allLastIndexes(Np,0);
173 if (m_env.inter0Rank() >= 0) {
177 unsigned int auxUInt = indexOfFirstWeight;
179 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
180 "failed MPI.Gather() for first indexes");
182 if (m_env.inter0Rank() == 0) {
185 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
186 "failed MPI.Gather() result for first indexes, at proc 0");
189 auxUInt = indexOfLastWeight;
191 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
192 "failed MPI.Gather() for last indexes");
194 if (m_env.inter0Rank() == 0) {
198 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
199 "failed MPI.Gather() result for last indexes, at proc 0");
206 if (m_env.inter0Rank() == 0) {
207 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
208 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
210 <<
", step " << m_currStep
211 <<
": original distribution of unified indexes in 'inter0Comm' is as follows"
213 for (
unsigned int r = 0; r < Np; ++r) {
214 *m_env.subDisplayFile() <<
" allFirstIndexes[" << r <<
"] = " << allFirstIndexes[r]
215 <<
" allLastIndexes[" << r <<
"] = " << allLastIndexes[r]
219 for (
unsigned int r = 0; r < (Np-1); ++r) {
222 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
226 for (
unsigned int r = 0; r < (Np-1); ++r) {
229 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
233 std::vector<unsigned int> origNumChainsPerNode (Np,0);
234 std::vector<unsigned int> origNumPositionsPerNode(Np,0);
236 for (
unsigned int i = 0; i < unifiedIndexCountersAtProc0Only.size(); ++i) {
237 if ((allFirstIndexes[r] <= i) &&
238 (i <= allLastIndexes[r] )) {
243 if ((r < (
int) Np ) &&
244 (allFirstIndexes[r] <= i) &&
245 (i <= allLastIndexes[r] )) {
249 std::cerr <<
"In MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
252 <<
", allFirstIndexes[r] = " << allFirstIndexes[r]
253 <<
", allLastIndexes[r] = " << allLastIndexes[r]
257 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
258 "wrong indexes or 'r' got too large");
261 if (unifiedIndexCountersAtProc0Only[i] != 0) {
262 origNumChainsPerNode [r] += 1;
263 origNumPositionsPerNode[r] += unifiedIndexCountersAtProc0Only[i];
270 exchangeStdVec.push_back(auxInfo);
276 unsigned int totalNumberOfChains = 0;
277 for (
unsigned int r = 0; r < Np; ++r) {
278 totalNumberOfChains += origNumChainsPerNode[r];
280 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
281 *m_env.subDisplayFile() <<
" KEY"
283 <<
", step " << m_currStep
285 <<
", totalNumberOfChains = " << totalNumberOfChains
290 unsigned int origMinPosPerNode = *std::min_element(origNumPositionsPerNode.begin(), origNumPositionsPerNode.end());
291 unsigned int origMaxPosPerNode = *std::max_element(origNumPositionsPerNode.begin(), origNumPositionsPerNode.end());
292 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
293 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
294 *m_env.subDisplayFile() <<
" KEY"
296 <<
", step " << m_currStep
297 <<
", origNumChainsPerNode[" << nodeId <<
"] = " << origNumChainsPerNode[nodeId]
298 <<
", origNumPositionsPerNode[" << nodeId <<
"] = " << origNumPositionsPerNode[nodeId]
302 double origRatioOfPosPerNode = ((double) origMaxPosPerNode ) / ((double) origMinPosPerNode);
303 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
304 *m_env.subDisplayFile() <<
" KEY"
306 <<
", step " << m_currStep
307 <<
", origRatioOfPosPerNode = " << origRatioOfPosPerNode
315 (m_env.numSubEnvironments() > 1 ) &&
316 (Np < totalNumberOfChains ) &&
323 m_env.fullComm().Barrier();
324 unsigned int tmpValue = result;
326 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
327 "failed MPI.Bcast() for 'result'");
328 if (m_env.inter0Rank() != 0) result = tmpValue;
330 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
331 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
333 <<
", step " << m_currStep
334 <<
": result = " << result
341 template <
class P_V,
class P_M>
350 std::vector<ExchangeInfoStruct>& exchangeStdVec,
353 if (m_env.inter0Rank() < 0)
return;
355 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
356 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
358 <<
", step " << m_currStep
362 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
363 if (m_env.inter0Rank() == 0) {
366 justBalance_proc0(currOptions,
372 #ifdef QUESO_HAS_GLPK
374 solveBIP_proc0(exchangeStdVec);
376 if (m_env.subDisplayFile()) {
377 *m_env.subDisplayFile() <<
"WARNING in MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
379 <<
"' has been requested, but this QUESO library has not been built with 'hdf5'"
380 <<
". Code will therefore process the algorithm id '" << 2
384 if (m_env.subRank() == 0) {
385 std::cerr <<
"WARNING in MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
387 <<
"' has been requested, but this QUESO library has not been built with 'hdf5'"
388 <<
". Code will therefore process the algorithm id '" << 2
392 justBalance_proc0(currOptions,
399 m_env.inter0Comm().Barrier();
404 unsigned int exchangeStdVecSize = exchangeStdVec.size();
406 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
407 "failed MPI.Bcast() for exchangeStdVec size");
408 if (m_env.inter0Rank() > 0) exchangeStdVec.resize(exchangeStdVecSize);
411 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
412 "failed MPI.Bcast() for exchangeStdVec data");
417 std::vector<unsigned int> finalNumChainsPerNode (Np,0);
418 std::vector<unsigned int> finalNumPositionsPerNode(Np,0);
419 unsigned int Nc = exchangeStdVec.size();
420 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
421 unsigned int nodeId = exchangeStdVec[chainId].finalNodeOfInitialPosition;
422 finalNumChainsPerNode [nodeId] += 1;
423 finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
429 unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
430 unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
431 double finalRatioOfPosPerNode = ((double) finalMaxPosPerNode) / ((double)finalMinPosPerNode);
434 std::vector<double> auxBuf(1,0.);
435 double minRatio = 0.;
436 auxBuf[0] = finalRatioOfPosPerNode;
438 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
439 "failed MPI.Allreduce() for min");
443 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
444 "failed minRatio sanity check");
446 double maxRatio = 0.;
447 auxBuf[0] = finalRatioOfPosPerNode;
449 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
450 "failed MPI.Allreduce() for max");
454 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
455 "failed maxRatio sanity check");
460 unsigned int finalNumChainsPerNodeSize = finalNumChainsPerNode.size();
462 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
463 "failed MPI.Bcast() for finalNumChainsPerNode size");
464 if (m_env.inter0Rank() > 0) finalNumChainsPerNode.resize(finalNumChainsPerNodeSize);
466 m_env.inter0Comm().Bcast((
void *) &finalNumChainsPerNode[0], (
int) finalNumChainsPerNodeSize,
RawValue_MPI_UNSIGNED, 0,
467 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
468 "failed MPI.Bcast() for finalNumChainsPerNode data");
474 mpiExchangePositions_inter0(prevChain,
477 prevLogLikelihoodValues,
480 finalNumChainsPerNode,
481 finalNumPositionsPerNode,
482 balancedLinkControl);
487 template <
class P_V,
class P_M>
490 unsigned int indexOfFirstWeight,
491 unsigned int indexOfLastWeight,
492 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
495 if (m_env.inter0Rank() < 0)
return;
497 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
498 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
500 <<
", step " << m_currStep
501 <<
": indexOfFirstWeight = " << indexOfFirstWeight
502 <<
", indexOfLastWeight = " << indexOfLastWeight
506 unsigned int subNumSamples = 0;
507 std::vector<unsigned int> unifiedIndexCountersAtAllProcs(0);
510 unsigned int resizeSize = unifiedIndexCountersAtProc0Only.size();
512 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
513 "failed MPI.Bcast() for resizeSize");
514 unifiedIndexCountersAtAllProcs.resize(resizeSize,0);
516 if (m_env.inter0Rank() == 0) unifiedIndexCountersAtAllProcs = unifiedIndexCountersAtProc0Only;
519 m_env.inter0Comm().Bcast((
void *) &unifiedIndexCountersAtAllProcs[0], (
int) unifiedIndexCountersAtAllProcs.size(),
RawValue_MPI_UNSIGNED, 0,
520 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
521 "failed MPI.Bcast() for unified index counters");
522 #if 0 // Use allgatherv ??? for subNumSamples instead
523 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
524 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
526 <<
", step " << m_currStep
529 for (
int r = 0; r < m_env.inter0Comm().NumProc(); ++r) {
530 *m_env.subDisplayFile() <<
" unifiedIndexCountersAtAllProcs[" << r <<
"] = " << unifiedIndexCountersAtAllProcs[r]
544 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
545 "invalid indexOfFirstWeight");
548 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
549 "invalid indexOfLastWeight");
551 for (
unsigned int i = indexOfFirstWeight; i <= indexOfLastWeight; ++i) {
552 subNumSamples += unifiedIndexCountersAtAllProcs[i];
555 std::vector<unsigned int> auxBuf(1,0);
557 unsigned int minModifiedSubNumSamples = 0;
558 auxBuf[0] = subNumSamples;
560 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
561 "failed MPI.Allreduce() for min");
563 unsigned int maxModifiedSubNumSamples = 0;
564 auxBuf[0] = subNumSamples;
566 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
567 "failed MPI.Allreduce() for max");
569 unsigned int sumModifiedSubNumSamples = 0;
570 auxBuf[0] = subNumSamples;
572 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
573 "failed MPI.Allreduce() for sum");
580 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
581 *m_env.subDisplayFile() <<
"KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
583 <<
", step " << m_currStep
584 <<
": subNumSamples = " << subNumSamples
585 <<
", unifiedIndexCountersAtAllProcs.size() = " << unifiedIndexCountersAtAllProcs.size()
587 *m_env.subDisplayFile() <<
"KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
589 <<
", step " << m_currStep
590 <<
": minModifiedSubNumSamples = " << minModifiedSubNumSamples
591 <<
", avgModifiedSubNumSamples = " << ((double) sumModifiedSubNumSamples)/((double) m_env.inter0Comm().NumProc())
592 <<
", maxModifiedSubNumSamples = " << maxModifiedSubNumSamples
596 unsigned int numberOfPositionsToGuaranteeForNode = subNumSamples;
597 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
598 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
600 <<
", step " << m_currStep
601 <<
": numberOfPositionsToGuaranteeForNode = " << numberOfPositionsToGuaranteeForNode
604 for (
unsigned int i = indexOfFirstWeight; i <= indexOfLastWeight; ++i) {
606 while (unifiedIndexCountersAtAllProcs[i] != 0) {
607 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 30)) {
608 *m_env.subDisplayFile() <<
", numberOfPositionsToGuaranteeForNode = " << numberOfPositionsToGuaranteeForNode
609 <<
", unifiedIndexCountersAtAllProcs[" << i
610 <<
"] = " << unifiedIndexCountersAtAllProcs[i]
613 if (unifiedIndexCountersAtAllProcs[i] < numberOfPositionsToGuaranteeForNode) {
619 numberOfPositionsToGuaranteeForNode -= unifiedIndexCountersAtAllProcs[i];
620 unifiedIndexCountersAtAllProcs[i] = 0;
622 else if ((unifiedIndexCountersAtAllProcs[i] == numberOfPositionsToGuaranteeForNode) &&
623 (unifiedIndexCountersAtAllProcs[i] > 0 )) {
630 unifiedIndexCountersAtAllProcs[i] -= numberOfPositionsToGuaranteeForNode;
631 numberOfPositionsToGuaranteeForNode = 0;
633 else if ((unifiedIndexCountersAtAllProcs[i] == numberOfPositionsToGuaranteeForNode) &&
634 (unifiedIndexCountersAtAllProcs[i] == 0 )) {
640 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
641 "should never get here");
647 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
648 "numberOfPositionsToGuaranteeForNode exited loop with wrong value");
651 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
652 *m_env.subDisplayFile() <<
"KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
654 <<
", step " << m_currStep
655 <<
": unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.
unbLinkedChains.size()
662 template <
class P_V,
class P_M>
666 const P_M& unifiedCovMatrix,
670 double& cumulativeRunTime,
671 unsigned int& cumulativeRejections,
675 m_env.fullComm().Barrier();
677 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
678 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
679 <<
": balancedLinkControl.balLinkedChains.size() = " << balancedLinkControl.
balLinkedChains.size()
683 P_V auxInitialPosition(m_vectorSpace.zeroVector());
684 double auxInitialLogPrior;
685 double auxInitialLogLikelihood;
687 unsigned int chainIdMax = 0;
688 if (m_env.inter0Rank() >= 0) {
693 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
694 "failed MPI.Bcast() for chainIdMax");
696 struct timeval timevalEntering;
698 iRC = gettimeofday(&timevalEntering, NULL);
701 if (m_env.inter0Rank() >= 0) {
702 unsigned int numberOfPositions = 0;
703 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
704 numberOfPositions += balancedLinkControl.
balLinkedChains[chainId].numberOfPositions;
707 std::vector<unsigned int> auxBuf(1,0);
709 unsigned int minNumberOfPositions = 0;
710 auxBuf[0] = numberOfPositions;
712 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
713 "failed MPI.Allreduce() for min");
715 unsigned int maxNumberOfPositions = 0;
716 auxBuf[0] = numberOfPositions;
718 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
719 "failed MPI.Allreduce() for max");
721 unsigned int sumNumberOfPositions = 0;
722 auxBuf[0] = numberOfPositions;
724 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
725 "failed MPI.Allreduce() for sum");
727 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
728 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
730 <<
", step " << m_currStep
731 <<
": chainIdMax = " << chainIdMax
732 <<
", numberOfPositions = " << numberOfPositions
733 <<
", at " << ctime(&timevalEntering.tv_sec)
735 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
737 <<
", step " << m_currStep
738 <<
": minNumberOfPositions = " << minNumberOfPositions
739 <<
", avgNumberOfPositions = " << ((double) sumNumberOfPositions)/((double) m_env.inter0Comm().NumProc())
740 <<
", maxNumberOfPositions = " << maxNumberOfPositions
746 if ((m_debugExponent == 1.) &&
747 (m_currStep == 10)) {
750 unsigned int cumulativeNumPositions = 0;
751 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
752 unsigned int tmpChainSize = 0;
753 if (m_env.inter0Rank() >= 0) {
755 auxInitialPosition = *(balancedLinkControl.
balLinkedChains[chainId].initialPosition);
756 auxInitialLogPrior = balancedLinkControl.
balLinkedChains[chainId].initialLogPrior;
757 auxInitialLogLikelihood = balancedLinkControl.
balLinkedChains[chainId].initialLogLikelihood;
758 tmpChainSize = balancedLinkControl.
balLinkedChains[chainId].numberOfPositions+1;
759 if ((m_env.subDisplayFile() ) &&
760 (m_env.displayVerbosity() >= 3)) {
761 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
763 <<
", step " << m_currStep
764 <<
", chainId = " << chainId
765 <<
" < " << chainIdMax
766 <<
": begin generating " << tmpChainSize
767 <<
" chain positions"
771 auxInitialPosition.mpiBcast(0, m_env.subComm());
773 #if 0 // For debug only
774 for (
int r = 0; r < m_env.subComm().NumProc(); ++r) {
775 if (r == m_env.subComm().MyPID()) {
776 std::cout <<
"Vector 'auxInitialPosition at rank " << r
777 <<
" has contents " << auxInitialPosition
780 m_env.subComm().Barrier();
787 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
788 "failed MPI.Bcast() for tmpChainSize");
793 m_options.m_prefix+
"tmp_chain");
800 "MLSamplingClass<P_V,P_M>::generateBalLinkedChains_all()",
801 "failed MPI.Bcast() for auxInitialLogPrior");
803 "MLSamplingClass<P_V,P_M>::generateBalLinkedChains_all()",
804 "failed MPI.Bcast() for auxInitialLogLikelihood");
810 auxInitialLogLikelihood,
813 &tmpLogLikelihoodValues,
814 &tmpLogTargetValues);
823 &tmpLogLikelihoodValues,
824 &tmpLogTargetValues);
828 cumulativeRunTime += mcRawInfo.
runTime;
831 if (m_env.inter0Rank() >= 0) {
832 if (m_env.exceptionalCircumstance()) {
833 if ((m_env.subDisplayFile() ) &&
834 (m_env.displayVerbosity() >= 0)) {
835 P_V tmpVec(m_vectorSpace.zeroVector());
836 for (
unsigned int i = 0; i < tmpLogLikelihoodValues.
subSequenceSize(); ++i) {
838 *m_env.subDisplayFile() <<
"DEBUG finalChain[" << cumulativeNumPositions+i <<
"] "
839 <<
"= tmpChain[" << i <<
"] = " << tmpVec
840 <<
", tmpLogLikelihoodValues[" << i <<
"] = " << tmpLogLikelihoodValues[i]
841 <<
", tmpLogTargetValues[" << i <<
"] = " << tmpLogTargetValues[i]
847 cumulativeNumPositions += tmpChainSize;
848 if (cumulativeNumPositions > 100) m_env.setExceptionalCircumstance(
false);
850 if ((m_env.subDisplayFile() ) &&
851 (m_env.displayVerbosity() >= 3)) {
852 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
854 <<
", step " << m_currStep
855 <<
", chainId = " << chainId
856 <<
" < " << chainIdMax
858 <<
" chain positions"
864 if (currLogLikelihoodValues) {
865 currLogLikelihoodValues->
append(tmpLogLikelihoodValues,1,tmpLogLikelihoodValues.
subSequenceSize()-1);
866 if ((m_env.subDisplayFile() ) &&
867 (m_env.displayVerbosity() >= 99) &&
869 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
871 <<
", step " << m_currStep
872 <<
", chainId = " << chainId
873 <<
", tmpLogLikelihoodValues.subSequenceSize() = " << tmpLogLikelihoodValues.
subSequenceSize()
874 <<
", tmpLogLikelihoodValues[0] = " << tmpLogLikelihoodValues[0]
875 <<
", tmpLogLikelihoodValues[1] = " << tmpLogLikelihoodValues[1]
876 <<
", currLogLikelihoodValues[0] = " << (*currLogLikelihoodValues)[0]
880 if (currLogTargetValues) {
889 struct timeval timevalBarrier;
890 iRC = gettimeofday(&timevalBarrier, NULL);
893 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
894 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
896 <<
", step " << m_currStep
897 <<
": ended chain loop after " << loopTime <<
" seconds"
898 <<
", calling fullComm().Barrier() at " << ctime(&timevalBarrier.tv_sec)
902 m_env.fullComm().Barrier();
904 struct timeval timevalLeaving;
905 iRC = gettimeofday(&timevalLeaving, NULL);
908 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
909 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
911 <<
", step " << m_currStep
912 <<
": after " << barrierTime <<
" seconds in fullComm().Barrier()"
913 <<
", at " << ctime(&timevalLeaving.tv_sec)
920 template <
class P_V,
class P_M>
924 const P_M& unifiedCovMatrix,
927 unsigned int indexOfFirstWeight,
934 double& cumulativeRunTime,
935 unsigned int& cumulativeRejections,
939 m_env.fullComm().Barrier();
941 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
942 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
943 <<
": unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.
unbLinkedChains.size()
944 <<
", indexOfFirstWeight = " << indexOfFirstWeight
948 P_V auxInitialPosition(m_vectorSpace.zeroVector());
949 double auxInitialLogPrior;
950 double auxInitialLogLikelihood;
952 unsigned int chainIdMax = 0;
953 if (m_env.inter0Rank() >= 0) {
958 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
959 "failed MPI.Bcast() for chainIdMax");
961 struct timeval timevalEntering;
963 iRC = gettimeofday(&timevalEntering, NULL);
966 if (m_env.inter0Rank() >= 0) {
967 unsigned int numberOfPositions = 0;
968 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
969 numberOfPositions += unbalancedLinkControl.
unbLinkedChains[chainId].numberOfPositions;
972 std::vector<unsigned int> auxBuf(1,0);
974 unsigned int minNumberOfPositions = 0;
975 auxBuf[0] = numberOfPositions;
977 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
978 "failed MPI.Allreduce() for min");
980 unsigned int maxNumberOfPositions = 0;
981 auxBuf[0] = numberOfPositions;
983 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
984 "failed MPI.Allreduce() for max");
986 unsigned int sumNumberOfPositions = 0;
987 auxBuf[0] = numberOfPositions;
989 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
990 "failed MPI.Allreduce() for sum");
992 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
993 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
995 <<
", step " << m_currStep
996 <<
": chainIdMax = " << chainIdMax
997 <<
", numberOfPositions = " << numberOfPositions
998 <<
", at " << ctime(&timevalEntering.tv_sec)
1000 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1002 <<
", step " << m_currStep
1003 <<
": minNumberOfPositions = " << minNumberOfPositions
1004 <<
", avgNumberOfPositions = " << ((double) sumNumberOfPositions)/((double) m_env.inter0Comm().NumProc())
1005 <<
", maxNumberOfPositions = " << maxNumberOfPositions
1009 if ((m_debugExponent == 1.) &&
1010 (m_currStep == 10)) {
1013 double expRatio = currExponent;
1014 if (prevExponent > 0.0) {
1015 expRatio /= prevExponent;
1017 unsigned int cumulativeNumPositions = 0;
1018 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
1019 unsigned int tmpChainSize = 0;
1020 if (m_env.inter0Rank() >= 0) {
1021 unsigned int auxIndex = unbalancedLinkControl.
unbLinkedChains[chainId].initialPositionIndexInPreviousChain - indexOfFirstWeight;
1023 auxInitialLogPrior = prevLogTargetValues[auxIndex] - prevLogLikelihoodValues[auxIndex];
1024 auxInitialLogLikelihood = expRatio * prevLogLikelihoodValues[auxIndex];
1025 tmpChainSize = unbalancedLinkControl.
unbLinkedChains[chainId].numberOfPositions+1;
1026 if ((m_env.subDisplayFile() ) &&
1027 (m_env.displayVerbosity() >= 3)) {
1028 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1030 <<
", step " << m_currStep
1031 <<
", chainId = " << chainId
1032 <<
" < " << chainIdMax
1033 <<
": begin generating " << tmpChainSize
1034 <<
" chain positions"
1038 auxInitialPosition.mpiBcast(0, m_env.subComm());
1040 #if 0 // For debug only
1041 for (
int r = 0; r < m_env.subComm().NumProc(); ++r) {
1042 if (r == m_env.subComm().MyPID()) {
1043 std::cout <<
"Vector 'auxInitialPosition at rank " << r
1044 <<
" has contents " << auxInitialPosition
1047 m_env.subComm().Barrier();
1054 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
1055 "failed MPI.Bcast() for tmpChainSize");
1060 m_options.m_prefix+
"tmp_chain");
1068 "MLSamplingClass<P_V,P_M>::generateUnbLinkedChains_all()",
1069 "failed MPI.Bcast() for auxInitialLogPrior");
1071 "MLSamplingClass<P_V,P_M>::generateUnbLinkedChains_all",
1072 "failed MPI.Bcast() for auxInitialLogLikelihood");
1077 auxInitialLogLikelihood,
1080 &tmpLogLikelihoodValues,
1081 &tmpLogTargetValues);
1090 &tmpLogLikelihoodValues,
1091 &tmpLogTargetValues);
1095 cumulativeRunTime += mcRawInfo.
runTime;
1098 if (m_env.inter0Rank() >= 0) {
1099 if (m_env.exceptionalCircumstance()) {
1100 if ((m_env.subDisplayFile() ) &&
1101 (m_env.displayVerbosity() >= 0)) {
1102 P_V tmpVec(m_vectorSpace.zeroVector());
1103 for (
unsigned int i = 0; i < tmpLogLikelihoodValues.
subSequenceSize(); ++i) {
1105 *m_env.subDisplayFile() <<
"DEBUG finalChain[" << cumulativeNumPositions+i <<
"] "
1106 <<
"= tmpChain[" << i <<
"] = " << tmpVec
1107 <<
", tmpLogLikelihoodValues[" << i <<
"] = " << tmpLogLikelihoodValues[i]
1108 <<
", tmpLogTargetValues[" << i <<
"] = " << tmpLogTargetValues[i]
1114 cumulativeNumPositions += tmpChainSize;
1115 if (cumulativeNumPositions > 100) m_env.setExceptionalCircumstance(
false);
1117 if ((m_env.subDisplayFile() ) &&
1118 (m_env.displayVerbosity() >= 3)) {
1119 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1121 <<
", step " << m_currStep
1122 <<
", chainId = " << chainId
1123 <<
" < " << chainIdMax
1125 <<
" chain positions"
1131 if (currLogLikelihoodValues) {
1132 currLogLikelihoodValues->
append(tmpLogLikelihoodValues,1,tmpLogLikelihoodValues.
subSequenceSize()-1);
1133 if ((m_env.subDisplayFile() ) &&
1134 (m_env.displayVerbosity() >= 99) &&
1136 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1138 <<
", step " << m_currStep
1139 <<
", chainId = " << chainId
1140 <<
", tmpLogLikelihoodValues.subSequenceSize() = " << tmpLogLikelihoodValues.
subSequenceSize()
1141 <<
", tmpLogLikelihoodValues[0] = " << tmpLogLikelihoodValues[0]
1142 <<
", tmpLogLikelihoodValues[1] = " << tmpLogLikelihoodValues[1]
1143 <<
", currLogLikelihoodValues[0] = " << (*currLogLikelihoodValues)[0]
1147 if (currLogTargetValues) {
1153 struct timeval timevalBarrier;
1154 iRC = gettimeofday(&timevalBarrier, NULL);
1157 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1158 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1160 <<
", step " << m_currStep
1161 <<
": ended chain loop after " << loopTime <<
" seconds"
1162 <<
", calling fullComm().Barrier() at " << ctime(&timevalBarrier.tv_sec)
1166 m_env.fullComm().Barrier();
1168 struct timeval timevalLeaving;
1169 iRC = gettimeofday(&timevalLeaving, NULL);
1172 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1173 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1175 <<
", step " << m_currStep
1176 <<
": after " << barrierTime <<
" seconds in fullComm().Barrier()"
1177 <<
", at " << ctime(&timevalLeaving.tv_sec)
1184 #ifdef QUESO_HAS_GLPK
1185 template <
class P_V,
class P_M>
1188 std::vector<ExchangeInfoStruct>& exchangeStdVec)
1190 if (m_env.inter0Rank() != 0)
return;
1193 struct timeval timevalBIP;
1194 iRC = gettimeofday(&timevalBIP, NULL);
1197 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
1198 unsigned int Nc = exchangeStdVec.size();
1200 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1201 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::solveBIP_proc0()"
1203 <<
", step " << m_currStep
1213 lp = glp_create_prob();
1214 glp_set_prob_name(lp,
"sample");
1219 unsigned int m = Nc+Np-1;
1220 unsigned int n = Nc*Np;
1221 unsigned int ne = Nc*Np + 2*Nc*(Np -1);
1223 glp_add_rows(lp, m);
1224 for (
int i = 1; i <= (int) Nc; ++i) {
1225 glp_set_row_bnds(lp, i, GLP_FX, 1.0, 1.0);
1226 glp_set_row_name(lp, i,
"");
1228 for (
int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1229 glp_set_row_bnds(lp, i, GLP_UP, 0.0, 0.0);
1230 glp_set_row_name(lp, i,
"");
1233 glp_add_cols(lp, n);
1234 for (
int j = 1; j <= (int) n; ++j) {
1236 glp_set_col_kind(lp, j, GLP_IV);
1237 glp_set_col_bnds(lp, j, GLP_DB, 0.0, 1.0);
1238 glp_set_col_name(lp, j,
"");
1241 glp_set_obj_dir(lp, GLP_MIN);
1242 for (
int chainId = 0; chainId <= (int) (Nc-1); ++chainId) {
1243 glp_set_obj_coef(lp, (chainId*Np)+1, exchangeStdVec[chainId].numberOfPositions);
1246 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1247 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1249 <<
", step " << m_currStep
1250 <<
": finished setting BIP rows and cols"
1260 std::vector<int > iVec(ne+1,0);
1261 std::vector<int > jVec(ne+1,0);
1262 std::vector<double> aVec(ne+1,0.);
1264 for (
int i = 1; i <= (int) Nc; ++i) {
1265 for (
int j = 1; j <= (int) Np; ++j) {
1267 jVec[coefId] = (i-1)*Np + j;
1272 for (
int i = 1; i <= (int) (Np-1); ++i) {
1273 for (
int j = 1; j <= (int) Nc; ++j) {
1274 iVec[coefId] = Nc+i;
1275 jVec[coefId] = (j-1)*Np + i;
1276 aVec[coefId] = -((double) exchangeStdVec[j-1].numberOfPositions);
1279 iVec[coefId] = Nc+i;
1280 jVec[coefId] = (j-1)*Np + i + 1;
1281 aVec[coefId] = exchangeStdVec[j-1].numberOfPositions;
1285 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1286 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1288 <<
", step " << m_currStep
1289 <<
": finished setting BIP constraint matrix"
1291 <<
", coefId = " << coefId
1297 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1298 "invalid final coefId");
1300 glp_load_matrix(lp, ne, &iVec[0], &jVec[0], &aVec[0]);
1302 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1303 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1305 <<
", step " << m_currStep
1306 <<
": finished loading BIP constraint matrix"
1307 <<
", glp_get_num_rows(lp) = " << glp_get_num_rows(lp)
1308 <<
", glp_get_num_cols(lp) = " << glp_get_num_cols(lp)
1309 <<
", glp_get_num_nz(lp) = " << glp_get_num_nz(lp)
1310 <<
", glp_get_num_int(lp) = " << glp_get_num_int(lp)
1311 <<
", glp_get_num_bin(lp) = " << glp_get_num_bin(lp)
1320 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1321 "invalid number of rows");
1325 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1326 "invalid number of columnss");
1330 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1331 "invalid number of nonzero constraint coefficients");
1335 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1336 "invalid number of integer structural variables");
1340 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1341 "invalid number of binary structural variables");
1346 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1347 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1348 int j = chainId*Np + nodeId + 1;
1350 glp_set_col_stat(lp, j, GLP_BS);
1353 glp_set_col_stat(lp, j, GLP_BS);
1358 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1359 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1360 int j = chainId*Np + nodeId + 1;
1361 int initialState = glp_mip_col_val(lp, j);
1365 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1366 "for nodeId = 0, initial state should be '1'");
1371 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1372 "for nodeId > 0, initial state should be '0'");
1377 for (
int i = 1; i <= (int) Nc; ++i) {
1378 glp_set_row_stat(lp, i, GLP_NS);
1380 for (
int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1381 glp_set_row_stat(lp, i, GLP_BS);
1391 BIP_routine_struct BIP_routine_info;
1392 BIP_routine_info.env = &m_env;
1393 BIP_routine_info.currLevel = m_currLevel;
1395 glp_iocp BIP_params;
1396 glp_init_iocp(&BIP_params);
1397 BIP_params.presolve = GLP_ON;
1402 int BIP_rc = glp_intopt(lp, &BIP_params);
1404 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1405 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1407 <<
", step " << m_currStep
1408 <<
": finished solving BIP"
1409 <<
", BIP_rc = " << BIP_rc
1415 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1416 "BIP returned rc != 0");
1421 int BIP_Status = glp_mip_status(lp);
1422 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1423 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1425 <<
", step " << m_currStep
1426 <<
": BIP_Status = " << BIP_Status
1430 switch (BIP_Status) {
1433 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1434 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1436 <<
", step " << m_currStep
1437 <<
": BIP solution is optimal"
1444 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1445 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1447 <<
", step " << m_currStep
1448 <<
": BIP solution is guaranteed to be 'only' feasible"
1456 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1457 "BIP has an undefined solution or has no solution");
1461 for (
int i = 1; i <= (int) Nc; ++i) {
1464 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1465 "row should have value 1 at solution");
1467 for (
int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1470 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1471 "row should have value 0 or should be negative at solution");
1477 std::vector<unsigned int> finalNumChainsPerNode (Np,0);
1478 std::vector<unsigned int> finalNumPositionsPerNode(Np,0);
1479 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1480 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1481 int j = chainId*Np + nodeId + 1;
1482 if (glp_mip_col_val(lp, j) == 0) {
1485 else if (glp_mip_col_val(lp, j) == 1) {
1488 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1489 "chain has already been taken care of");
1490 exchangeStdVec[chainId].finalNodeOfInitialPosition = nodeId;
1491 finalNumChainsPerNode [nodeId] += 1;
1492 finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
1497 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1498 "control variable should be either '0' or '1'");
1503 unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1504 unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1506 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1507 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1509 <<
", step " << m_currStep
1510 <<
": finished preparing output information"
1517 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1518 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1520 <<
", step " << m_currStep
1521 <<
": solution gives the following redistribution"
1523 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1524 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1526 <<
", step " << m_currStep
1527 <<
", finalNumChainsPerNode[" << nodeId <<
"] = " << finalNumChainsPerNode[nodeId]
1528 <<
", finalNumPositionsPerNode[" << nodeId <<
"] = " << finalNumPositionsPerNode[nodeId]
1531 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1533 <<
", step " << m_currStep
1534 <<
", finalRatioOfPosPerNode = " << ((double) finalMaxPosPerNode) / ((double)finalMinPosPerNode)
1543 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1544 "Invalid objective value");
1546 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
1547 UQ_FATAL_TEST_MACRO(finalNumPositionsPerNode[nodeId-1] < finalNumPositionsPerNode[nodeId],
1549 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1550 "Next node should have a number of positions equal or less than the current node");
1553 for (
int i = (
int) (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1554 unsigned int nodeId = i - Nc;
1555 int diff = ((int) finalNumPositionsPerNode[nodeId]) - ((int) finalNumPositionsPerNode[nodeId-1]);
1558 "MLSampling<P_V,P_M>::solveBIP_proc0()",
1565 glp_delete_prob(lp);
1568 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1569 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::solveBIP_proc0()"
1571 <<
", step " << m_currStep
1572 <<
", after " << bipRunTime <<
" seconds"
1578 #endif // QUESO_HAS_GLPK
1580 template <
class P_V,
class P_M>
1584 std::vector<ExchangeInfoStruct>& exchangeStdVec)
1586 if (m_env.inter0Rank() != 0)
return;
1589 struct timeval timevalBal;
1590 iRC = gettimeofday(&timevalBal, NULL);
1593 unsigned int Np = m_env.numSubEnvironments();
1594 unsigned int Nc = exchangeStdVec.size();
1596 std::vector<ExchangeInfoStruct> currExchangeStdVec(Nc);
1597 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1598 currExchangeStdVec[chainId] = exchangeStdVec[chainId];
1599 currExchangeStdVec[chainId].finalNodeOfInitialPosition = currExchangeStdVec[chainId].originalNodeOfInitialPosition;
1605 unsigned int iterIdMax = 0;
1606 std::vector<unsigned int> currNumChainsPerNode (Np,0);
1607 std::vector<unsigned int> currNumPositionsPerNode(Np,0);
1608 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1609 unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1610 currNumChainsPerNode [nodeId] += 1;
1611 currNumPositionsPerNode[nodeId] += currExchangeStdVec[chainId].numberOfPositions;
1612 iterIdMax += currExchangeStdVec[chainId].numberOfPositions;
1614 unsigned int currMinPosPerNode = *std::min_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1615 unsigned int currMaxPosPerNode = *std::max_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1616 double currRatioOfPosPerNode = ((double) currMaxPosPerNode ) / ((double) currMinPosPerNode);
1623 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1624 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1626 <<
", step " << m_currStep
1627 <<
", iter " << iterId
1628 <<
", currRatioOfPosPerNode = " << currRatioOfPosPerNode
1630 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1631 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1633 <<
", step " << m_currStep
1634 <<
", iter " << iterId
1635 <<
", currNumChainsPerNode[" << nodeId <<
"] = " << currNumChainsPerNode[nodeId]
1636 <<
", currNumPositionsPerNode[" << nodeId <<
"] = " << currNumPositionsPerNode[nodeId]
1641 std::vector<std::vector<double> > vectorOfChainSizesPerNode(Np);
1642 while ((iterId < (
int) iterIdMax ) &&
1649 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1650 vectorOfChainSizesPerNode[nodeId].clear();
1652 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1653 unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1654 vectorOfChainSizesPerNode[nodeId].push_back(currExchangeStdVec[chainId].numberOfPositions);
1657 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1658 std::sort(vectorOfChainSizesPerNode[nodeId].begin(), vectorOfChainSizesPerNode[nodeId].end());
1659 UQ_FATAL_TEST_MACRO(vectorOfChainSizesPerNode[nodeId].size() != currNumChainsPerNode[nodeId],
1661 "MLSampling<P_V,P_M>::justBalance_proc0()",
1662 "inconsistent number of chains in node");
1668 unsigned int currBiggestAmountOfPositionsPerNode = currNumPositionsPerNode[0];
1669 unsigned int currSmallestAmountOfPositionsPerNode = currNumPositionsPerNode[0];
1670 unsigned int currNodeWithMostPositions = 0;
1671 unsigned int currNodeWithLeastPositions = 0;
1672 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1673 if (currNumPositionsPerNode[nodeId] > currBiggestAmountOfPositionsPerNode) {
1674 currBiggestAmountOfPositionsPerNode = currNumPositionsPerNode[nodeId];
1675 currNodeWithMostPositions = nodeId;
1677 if (currNumPositionsPerNode[nodeId] < currSmallestAmountOfPositionsPerNode) {
1678 currSmallestAmountOfPositionsPerNode = currNumPositionsPerNode[nodeId];
1679 currNodeWithLeastPositions = nodeId;
1683 UQ_FATAL_TEST_MACRO(currMinPosPerNode != currNumPositionsPerNode[currNodeWithLeastPositions],
1685 "MLSampling<P_V,P_M>::justBalance_proc0()",
1686 "inconsistent currMinPosPerNode");
1688 UQ_FATAL_TEST_MACRO(currMaxPosPerNode != currNumPositionsPerNode[currNodeWithMostPositions],
1690 "MLSampling<P_V,P_M>::justBalance_proc0()",
1691 "inconsistent currMaxPosPerNode");
1693 unsigned int numberOfPositionsToMove = vectorOfChainSizesPerNode[currNodeWithMostPositions][0];
1695 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1696 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1698 <<
", step " << m_currStep
1699 <<
", iter " << iterId
1700 <<
", before update"
1701 <<
", node w/ most pos is "
1702 << currNodeWithMostPositions <<
"(cs=" << currNumChainsPerNode[currNodeWithMostPositions ] <<
", ps=" << currNumPositionsPerNode[currNodeWithMostPositions ] <<
")"
1703 <<
", node w/ least pos is "
1704 << currNodeWithLeastPositions <<
"(cs=" << currNumChainsPerNode[currNodeWithLeastPositions] <<
", ps=" << currNumPositionsPerNode[currNodeWithLeastPositions] <<
")"
1705 <<
", number of pos to move = " << numberOfPositionsToMove
1712 std::vector<ExchangeInfoStruct> newExchangeStdVec(Nc);
1713 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1714 newExchangeStdVec[chainId] = currExchangeStdVec[chainId];
1717 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1718 if ((newExchangeStdVec[chainId].finalNodeOfInitialPosition == (
int) currNodeWithMostPositions) &&
1719 (newExchangeStdVec[chainId].numberOfPositions == numberOfPositionsToMove )) {
1720 newExchangeStdVec[chainId].finalNodeOfInitialPosition = currNodeWithLeastPositions;
1728 std::vector<unsigned int> newNumChainsPerNode (Np,0);
1729 std::vector<unsigned int> newNumPositionsPerNode(Np,0);
1730 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1731 unsigned int nodeId = newExchangeStdVec[chainId].finalNodeOfInitialPosition;
1732 newNumChainsPerNode [nodeId] += 1;
1733 newNumPositionsPerNode[nodeId] += newExchangeStdVec[chainId].numberOfPositions;
1736 unsigned int newBiggestAmountOfPositionsPerNode = newNumPositionsPerNode[0];
1737 unsigned int newSmallestAmountOfPositionsPerNode = newNumPositionsPerNode[0];
1738 unsigned int newNodeWithMostPositions = 0;
1739 unsigned int newNodeWithLeastPositions = 0;
1740 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1741 if (newNumPositionsPerNode[nodeId] > newBiggestAmountOfPositionsPerNode) {
1742 newBiggestAmountOfPositionsPerNode = newNumPositionsPerNode[nodeId];
1743 newNodeWithMostPositions = nodeId;
1745 if (newNumPositionsPerNode[nodeId] < newSmallestAmountOfPositionsPerNode) {
1746 newSmallestAmountOfPositionsPerNode = newNumPositionsPerNode[nodeId];
1747 newNodeWithLeastPositions = nodeId;
1751 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1752 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1754 <<
", step " << m_currStep
1755 <<
", iter " << iterId
1757 <<
", node w/ most pos is "
1758 << newNodeWithMostPositions <<
"(cs=" << newNumChainsPerNode[newNodeWithMostPositions ] <<
", ps=" << newNumPositionsPerNode[newNodeWithMostPositions ] <<
")"
1759 <<
", node w/ least pos is "
1760 << newNodeWithLeastPositions <<
"(cs=" << newNumChainsPerNode[newNodeWithLeastPositions] <<
", ps=" << newNumPositionsPerNode[newNodeWithLeastPositions] <<
")"
1764 unsigned int newMinPosPerNode = *std::min_element(newNumPositionsPerNode.begin(), newNumPositionsPerNode.end());
1765 unsigned int newMaxPosPerNode = *std::max_element(newNumPositionsPerNode.begin(), newNumPositionsPerNode.end());
1766 double newRatioOfPosPerNode = ((double) newMaxPosPerNode ) / ((double) newMinPosPerNode);
1768 UQ_FATAL_TEST_MACRO(newMinPosPerNode != newNumPositionsPerNode[newNodeWithLeastPositions],
1770 "MLSampling<P_V,P_M>::justBalance_proc0()",
1771 "inconsistent newMinPosPerNode");
1775 "MLSampling<P_V,P_M>::justBalance_proc0()",
1776 "inconsistent newMaxPosPerNode");
1778 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1779 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1781 <<
", step " << m_currStep
1782 <<
", iter " << iterId
1783 <<
", newMaxPosPerNode = " << newMaxPosPerNode
1784 <<
", newMinPosPerNode = " << newMinPosPerNode
1785 <<
", newRatioOfPosPerNode = " << newRatioOfPosPerNode
1787 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1788 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1790 <<
", step " << m_currStep
1791 <<
", iter " << iterId
1792 <<
", newNumChainsPerNode[" << nodeId <<
"] = " << newNumChainsPerNode [nodeId]
1793 <<
", newNumPositionsPerNode[" << nodeId <<
"] = " << newNumPositionsPerNode[nodeId]
1801 if (newRatioOfPosPerNode > currRatioOfPosPerNode) {
1808 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1809 currNumChainsPerNode [nodeId] = 0;
1810 currNumPositionsPerNode[nodeId] = 0;
1812 currRatioOfPosPerNode = newRatioOfPosPerNode;
1813 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1814 currExchangeStdVec[chainId] = newExchangeStdVec[chainId];
1815 unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1816 currNumChainsPerNode [nodeId] += 1;
1817 currNumPositionsPerNode[nodeId] += currExchangeStdVec[chainId].numberOfPositions;
1819 currMinPosPerNode = *std::min_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1820 currMaxPosPerNode = *std::max_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1821 currRatioOfPosPerNode = ((double) currMaxPosPerNode ) / ((double) currMinPosPerNode);
1827 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1828 exchangeStdVec[chainId].finalNodeOfInitialPosition = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1834 std::vector<unsigned int> finalNumChainsPerNode (Np,0);
1835 std::vector<unsigned int> finalNumPositionsPerNode(Np,0);
1836 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1837 unsigned int nodeId = exchangeStdVec[chainId].finalNodeOfInitialPosition;
1838 finalNumChainsPerNode [nodeId] += 1;
1839 finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
1841 unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1842 unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1843 double finalRatioOfPosPerNode = ((double) finalMaxPosPerNode ) / ((double) finalMinPosPerNode);
1845 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1846 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1848 <<
", step " << m_currStep
1849 <<
": solution gives the following redistribution"
1851 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1852 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1854 <<
", step " << m_currStep
1855 <<
", finalNumChainsPerNode[" << nodeId <<
"] = " << finalNumChainsPerNode[nodeId]
1856 <<
", finalNumPositionsPerNode[" << nodeId <<
"] = " << finalNumPositionsPerNode[nodeId]
1859 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1861 <<
", step " << m_currStep
1862 <<
", finalRatioOfPosPerNode = " << finalRatioOfPosPerNode
1870 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1871 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::justBalance_proc0()"
1873 <<
", step " << m_currStep
1874 <<
", iterId = " << iterId
1875 <<
", currRatioOfPosPerNode = " << currRatioOfPosPerNode
1876 <<
", after " << balRunTime <<
" seconds"
1883 template <
class P_V,
class P_M>
1887 double prevExponent,
1888 double currExponent,
1891 const std::vector<ExchangeInfoStruct>& exchangeStdVec,
1892 const std::vector<unsigned int>& finalNumChainsPerNode,
1893 const std::vector<unsigned int>& finalNumPositionsPerNode,
1896 if (m_env.inter0Rank() < 0) {
1900 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
1901 unsigned int Nc = exchangeStdVec.size();
1903 double expRatio = currExponent;
1904 if (prevExponent > 0.0) {
1905 expRatio /= prevExponent;
1908 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1909 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1911 <<
", step " << m_currStep
1922 for (
unsigned int r = 0; r < Np; ++r) {
1926 unsigned int numberOfInitialPositionsNodeRAlreadyHas = 0;
1927 std::vector<unsigned int> numberOfInitialPositionsNodeRHasToReceiveFromNode(Np,0);
1928 std::vector<unsigned int> indexesOfInitialPositionsNodeRHasToReceiveFromMe(0);
1930 unsigned int sumOfChainLenghtsNodeRAlreadyHas = 0;
1931 std::vector<unsigned int> chainLenghtsNodeRHasToInherit(0);
1933 for (
unsigned int i = 0; i < Nc; ++i) {
1934 if (exchangeStdVec[i].finalNodeOfInitialPosition == (
int) r) {
1935 if (exchangeStdVec[i].originalNodeOfInitialPosition == (
int) r) {
1936 numberOfInitialPositionsNodeRAlreadyHas++;
1937 sumOfChainLenghtsNodeRAlreadyHas += exchangeStdVec[i].numberOfPositions;
1940 numberOfInitialPositionsNodeRHasToReceiveFromNode[exchangeStdVec[i].originalNodeOfInitialPosition]++;
1941 chainLenghtsNodeRHasToInherit.push_back(exchangeStdVec[i].numberOfPositions);
1942 if (m_env.inter0Rank() == exchangeStdVec[i].originalNodeOfInitialPosition) {
1943 indexesOfInitialPositionsNodeRHasToReceiveFromMe.push_back(exchangeStdVec[i].originalIndexOfInitialPosition);
1949 unsigned int totalNumberOfInitialPositionsNodeRHasToReceive = 0;
1950 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1951 totalNumberOfInitialPositionsNodeRHasToReceive += numberOfInitialPositionsNodeRHasToReceiveFromNode[nodeId];
1954 unsigned int totalNumberOfChainLenghtsNodeRHasToInherit = chainLenghtsNodeRHasToInherit.size();
1955 unsigned int totalSumOfChainLenghtsNodeRHasToInherit = 0;
1956 for (
unsigned int i = 0; i < totalNumberOfChainLenghtsNodeRHasToInherit; ++i) {
1957 totalSumOfChainLenghtsNodeRHasToInherit += chainLenghtsNodeRHasToInherit[i];
1963 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1964 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1966 <<
", step " << m_currStep
1968 <<
", finalNumChainsPerNode[r] = " << finalNumChainsPerNode[r]
1969 <<
", totalNumberOfInitialPositionsNodeRHasToReceive = " << totalNumberOfInitialPositionsNodeRHasToReceive
1970 <<
", numberOfInitialPositionsNodeRAlreadyHas = " << numberOfInitialPositionsNodeRAlreadyHas
1972 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1974 <<
", step " << m_currStep
1976 <<
", finalNumPositionsPerNode[r] = " << finalNumPositionsPerNode[r]
1977 <<
", totalSumOfChainLenghtsNodeRHasToInherit = " << totalSumOfChainLenghtsNodeRHasToInherit
1978 <<
", sumOfChainLenghtsNodeRAlreadyHas = " << sumOfChainLenghtsNodeRAlreadyHas
1985 UQ_FATAL_TEST_MACRO(indexesOfInitialPositionsNodeRHasToReceiveFromMe.size() != numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()],
1987 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
1988 "inconsistent number of initial positions to send to node 'r'");
1990 UQ_FATAL_TEST_MACRO(finalNumChainsPerNode[r] != (totalNumberOfInitialPositionsNodeRHasToReceive + numberOfInitialPositionsNodeRAlreadyHas),
1992 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
1993 "inconsistent number of chains in node 'r'");
1995 UQ_FATAL_TEST_MACRO(finalNumPositionsPerNode[r] != (totalSumOfChainLenghtsNodeRHasToInherit + sumOfChainLenghtsNodeRAlreadyHas),
1997 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
1998 "inconsistent sum of chain lenghts in node 'r'");
2000 UQ_FATAL_TEST_MACRO(totalNumberOfInitialPositionsNodeRHasToReceive != totalNumberOfChainLenghtsNodeRHasToInherit,
2002 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
2003 "inconsistent on total number of initial positions to receive in node 'r'");
2006 indexesOfInitialPositionsNodeRHasToReceiveFromMe.resize(numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()]);
2007 chainLenghtsNodeRHasToInherit.resize (totalSumOfChainLenghtsNodeRHasToInherit);
2012 unsigned int dimSize = m_vectorSpace.dimLocal();
2013 unsigned int nValuesPerInitialPosition = dimSize + 2;
2014 P_V auxInitialPosition(m_vectorSpace.zeroVector());
2015 std::vector<double> sendbuf(0);
2016 unsigned int sendcnt = 0;
2017 if (m_env.inter0Rank() != (int) r) {
2018 sendcnt = numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()] * nValuesPerInitialPosition;
2019 sendbuf.resize(sendcnt);
2020 for (
unsigned int i = 0; i < numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()]; ++i) {
2021 unsigned int auxIndex = indexesOfInitialPositionsNodeRHasToReceiveFromMe[i];
2023 for (
unsigned int j = 0; j < dimSize; ++j) {
2024 sendbuf[i*nValuesPerInitialPosition + j] = auxInitialPosition[j];
2026 sendbuf[i*nValuesPerInitialPosition + dimSize] = prevLogLikelihoodValues[auxIndex];
2027 sendbuf[i*nValuesPerInitialPosition + dimSize + 1] = prevLogTargetValues[auxIndex];
2031 std::vector<double> recvbuf(0);
2032 std::vector<int> recvcnts(Np,0);
2033 if (m_env.inter0Rank() == (int) r) {
2034 recvbuf.resize(totalNumberOfInitialPositionsNodeRHasToReceive * nValuesPerInitialPosition);
2035 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
2036 recvcnts[nodeId] = numberOfInitialPositionsNodeRHasToReceiveFromNode[nodeId] * nValuesPerInitialPosition;
2040 std::vector<int> displs(Np,0);
2041 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
2042 displs[nodeId] = displs[nodeId-1] + recvcnts[nodeId-1];
2046 if (m_env.inter0Rank() == r) {
2048 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0(1)",
2049 "failed MPI.Gatherv()");
2053 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0(2)",
2054 "failed MPI.Gatherv()");
2058 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
2059 "failed MPI.Gatherv()");
2071 if (m_env.inter0Rank() == (int) r) {
2073 unsigned int auxIndex = 0;
2075 for (
unsigned int i = 0; i < Nc; ++i) {
2076 if ((exchangeStdVec[i].finalNodeOfInitialPosition == (
int) r) &&
2077 (exchangeStdVec[i].originalNodeOfInitialPosition == (
int) r)) {
2078 unsigned int originalIndex = exchangeStdVec[i].originalIndexOfInitialPosition;
2080 balancedLinkControl.
balLinkedChains[auxIndex].initialPosition =
new P_V(auxInitialPosition);
2081 balancedLinkControl.
balLinkedChains[auxIndex].initialLogPrior = prevLogTargetValues[originalIndex] - prevLogLikelihoodValues[originalIndex];
2082 balancedLinkControl.
balLinkedChains[auxIndex].initialLogLikelihood = expRatio*prevLogLikelihoodValues[originalIndex];
2083 balancedLinkControl.
balLinkedChains[auxIndex].numberOfPositions = exchangeStdVec[i].numberOfPositions;
2088 for (
unsigned int i = 0; i < totalNumberOfInitialPositionsNodeRHasToReceive; ++i) {
2089 for (
unsigned int j = 0; j < dimSize; ++j) {
2090 auxInitialPosition[j] = recvbuf[i*nValuesPerInitialPosition + j];
2092 balancedLinkControl.
balLinkedChains[auxIndex].initialPosition =
new P_V(auxInitialPosition);
2093 double prevLogLikelihood = recvbuf[i*nValuesPerInitialPosition + dimSize];
2094 double prevLogTarget = recvbuf[i*nValuesPerInitialPosition + dimSize + 1];
2095 balancedLinkControl.
balLinkedChains[auxIndex].initialLogPrior = prevLogTarget - prevLogLikelihood;
2096 balancedLinkControl.
balLinkedChains[auxIndex].initialLogLikelihood = expRatio*prevLogLikelihood;
2097 balancedLinkControl.
balLinkedChains[auxIndex].numberOfPositions = chainLenghtsNodeRHasToInherit[i];
2102 m_env.inter0Comm().Barrier();
2105 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2106 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
2108 <<
", step " << m_currStep
2116 template <
class P_V,
class P_M>
2119 double currExponent,
2125 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2126 *m_env.subDisplayFile() <<
"\n CHECKPOINTING initiating at level " << m_currLevel
2127 <<
"\n" << std::endl;
2134 unsigned int quantity2 = currLogLikelihoodValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2135 unsigned int quantity3 = currLogTargetValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2136 if (m_env.inter0Rank() >= 0) {
2139 "MLSampling<P_V,P_M>::checkpointML()",
2140 "number of evidence factors is not consistent");
2143 "MLSampling<P_V,P_M>::checkpointML()",
2144 "quantity2 is not consistent");
2147 "MLSampling<P_V,P_M>::checkpointML()",
2148 "quantity3 is not consistent");
2151 if (m_env.fullRank() == 0) {
2152 std::ofstream* ofsVar =
new std::ofstream((m_options.m_restartOutput_baseNameForFiles +
"Control.txt").c_str(),
2153 std::ofstream::out | std::ofstream::trunc);
2154 *ofsVar << m_currLevel << std::endl
2155 << m_vectorSpace.dimGlobal() << std::endl
2156 << currExponent << std::endl
2157 << currEta << std::endl
2158 << quantity1 << std::endl;
2159 unsigned int savedPrecision = ofsVar->precision();
2160 ofsVar->precision(16);
2161 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2162 *ofsVar << m_logEvidenceFactors[i] << std::endl;
2164 ofsVar->precision(savedPrecision);
2165 *ofsVar <<
"COMPLETE" << std::endl;
2169 m_env.fullComm().Barrier();
2174 char levelSufix[256];
2177 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2178 *m_env.subDisplayFile() <<
"\n CHECKPOINTING chain at level " << m_currLevel
2179 <<
"\n" << std::endl;
2181 currChain.
unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles +
"Chain_l" + levelSufix,
2182 m_options.m_restartOutput_fileType);
2183 m_env.fullComm().Barrier();
2185 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2186 *m_env.subDisplayFile() <<
"\n CHECKPOINTING like at level " << m_currLevel
2187 <<
"\n" << std::endl;
2189 currLogLikelihoodValues.
unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles +
"LogLike_l" + levelSufix,
2190 m_options.m_restartOutput_fileType);
2191 m_env.fullComm().Barrier();
2193 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2194 *m_env.subDisplayFile() <<
"\n CHECKPOINTING target at level " << m_currLevel
2195 <<
"\n" << std::endl;
2197 currLogTargetValues.
unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles +
"LogTarget_l" + levelSufix,
2198 m_options.m_restartOutput_fileType);
2199 m_env.fullComm().Barrier();
2204 if (m_env.fullRank() == 0) {
2205 std::ofstream* ofsVar =
new std::ofstream((m_options.m_restartOutput_baseNameForFiles +
"Control_l" + levelSufix +
".txt").c_str(),
2206 std::ofstream::out | std::ofstream::trunc);
2207 *ofsVar << m_currLevel << std::endl
2208 << m_vectorSpace.dimGlobal() << std::endl
2209 << currExponent << std::endl
2210 << currEta << std::endl
2211 << quantity1 << std::endl;
2212 unsigned int savedPrecision = ofsVar->precision();
2213 ofsVar->precision(16);
2214 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2215 *ofsVar << m_logEvidenceFactors[i] << std::endl;
2217 ofsVar->precision(savedPrecision);
2218 *ofsVar <<
"COMPLETE" << std::endl;
2222 m_env.fullComm().Barrier();
2224 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2225 *m_env.subDisplayFile() <<
"\n CHECKPOINTING done at level " << m_currLevel
2226 <<
"\n" << std::endl;
2232 template <
class P_V,
class P_M>
2235 double& currExponent,
2241 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2242 *m_env.subDisplayFile() <<
"\n RESTARTING initiating at level " << m_currLevel
2243 <<
"\n" << std::endl;
2249 unsigned int vectorSpaceDim = 0;
2250 unsigned int quantity1 = 0;
2251 std::string checkingString(
"");
2252 if (m_env.fullRank() == 0) {
2253 std::ifstream* ifsVar =
new std::ifstream((m_options.m_restartInput_baseNameForFiles +
"Control.txt").c_str(),
2259 unsigned int numLines = std::count(std::istreambuf_iterator<char>(*ifsVar),
2260 std::istreambuf_iterator<char>(),
2262 ifsVar->seekg(0,std::ios_base::beg);
2263 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2264 *m_env.subDisplayFile() <<
"Restart input file has " << numLines
2272 *ifsVar >> m_currLevel;
2275 "MLSampling<P_V,P_M>::restartML()",
2276 "number of lines read is different than pre-established number of lines in control file");
2278 m_logEvidenceFactors.clear();
2279 m_logEvidenceFactors.resize(m_currLevel,0.);
2280 *ifsVar >> vectorSpaceDim
2284 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2285 *ifsVar >> m_logEvidenceFactors[i];
2287 *ifsVar >> checkingString;
2290 "MLSampling<P_V,P_M>::restartML()",
2291 "control txt input file is not complete");
2293 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2294 *m_env.subDisplayFile() <<
"Restart input file has the following information:"
2295 <<
"\n m_currLevel = " << m_currLevel
2296 <<
"\n vectorSpaceDim = " << vectorSpaceDim
2297 <<
"\n currExponent = " << currExponent
2298 <<
"\n currEta = " << currEta
2299 <<
"\n quantity1 = " << quantity1;
2300 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2301 *m_env.subDisplayFile() <<
"\n [" << i <<
"] = " << m_logEvidenceFactors[i];
2303 *m_env.subDisplayFile() << std::endl;
2306 #if 0 // For debug only
2307 std::string tmpString;
2308 for (
unsigned int i = 0; i < 2; ++i) {
2309 *ifsVar >> tmpString;
2310 std::cout <<
"Just read '" << tmpString <<
"'" << std::endl;
2312 while ((lineId < numLines) && (ifsVar->eof() ==
false)) {
2314 ifsVar->ignore(maxCharsPerLine,
'\n');
2319 m_env.fullComm().Barrier();
2324 unsigned int tmpUint = (
unsigned int) m_currLevel;
2326 "MLSampling<P_V,P_M>::restartML()",
2327 "failed MPI.Bcast() for m_currLevel");
2328 if (m_env.fullRank() != 0) {
2329 m_currLevel = tmpUint;
2336 if (m_env.fullRank() == 0) {
2337 tmpData[0] = vectorSpaceDim;
2338 tmpData[1] = currExponent;
2339 tmpData[2] = currEta;
2340 tmpData[3] = quantity1;
2341 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2342 tmpData[4+i] = m_logEvidenceFactors[i];
2346 m_logEvidenceFactors.clear();
2347 m_logEvidenceFactors.resize(m_currLevel,0.);
2349 m_env.fullComm().Bcast((
void *) &tmpData[0], (
int) tmpData.size(),
RawValue_MPI_DOUBLE, 0,
2350 "MLSampling<P_V,P_M>::restartML()",
2351 "failed MPI.Bcast() for rest of information read from input file");
2352 if (m_env.fullRank() != 0) {
2353 vectorSpaceDim = tmpData[0];
2354 currExponent = tmpData[1];
2355 currEta = tmpData[2];
2356 quantity1 = tmpData[3];
2357 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2358 m_logEvidenceFactors[i] = tmpData[4+i];
2367 "MLSampling<P_V,P_M>::restartML()",
2368 "read vector space dimension is not consistent");
2371 "MLSampling<P_V,P_M>::restartML()",
2372 "read currExponent is not consistent");
2375 "MLSampling<P_V,P_M>::restartML()",
2376 "read size of chain should be a multiple of the number of subenvironments");
2377 unsigned int subSequenceSize = 0;
2378 subSequenceSize = ((double) quantity1) / ((double) m_env.numSubEnvironments());
2380 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2381 *m_env.subDisplayFile() <<
"Restart input file has the following information"
2382 <<
": subSequenceSize = " << subSequenceSize
2389 char levelSufix[256];
2392 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2393 *m_env.subDisplayFile() <<
"\n RESTARTING chain at level " << m_currLevel
2394 <<
"\n" << std::endl;
2396 currChain.
unifiedReadContents(m_options.m_restartInput_baseNameForFiles +
"Chain_l" + levelSufix,
2397 m_options.m_restartInput_fileType,
2399 m_env.fullComm().Barrier();
2401 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2402 *m_env.subDisplayFile() <<
"\n RESTARTING like at level " << m_currLevel
2403 <<
"\n" << std::endl;
2405 currLogLikelihoodValues.
unifiedReadContents(m_options.m_restartInput_baseNameForFiles +
"LogLike_l" + levelSufix,
2406 m_options.m_restartInput_fileType,
2408 m_env.fullComm().Barrier();
2410 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2411 *m_env.subDisplayFile() <<
"\n RESTARTING target at level " << m_currLevel
2412 <<
"\n" << std::endl;
2414 currLogTargetValues.
unifiedReadContents(m_options.m_restartInput_baseNameForFiles +
"LogTarget_l" + levelSufix,
2415 m_options.m_restartInput_fileType,
2417 m_env.fullComm().Barrier();
2419 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2420 *m_env.subDisplayFile() <<
"\n RESTARTING done at level " << m_currLevel
2421 <<
"\n" << std::endl;
2427 template <
class P_V,
class P_M>
2431 unsigned int& unifiedRequestedNumSamples,
2436 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2437 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateSequence()"
2439 <<
", currOptions.m_rawChainSize = " << currOptions.
m_rawChainSize
2444 struct timeval timevalLevel;
2445 iRC = gettimeofday(&timevalLevel, NULL);
2448 if (m_env.inter0Rank() >= 0) {
2451 "MLSampling<P_V,P_M>::generateSequence()",
2452 "failed MPI.Allreduce() for requested num samples in level 0");
2459 currLogLikelihoodValues.
setName(currOptions.
m_prefix +
"rawLogLikelihood");
2466 P_V auxVec(m_vectorSpace.zeroVector());
2470 bool outOfSupport =
true;
2473 m_priorRv.realizer().realization(auxVec);
2474 if (m_numDisabledParameters > 0) {
2475 unsigned int disabledCounter = 0;
2476 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2477 if (m_parameterEnabledStatus[paramId] ==
false) {
2483 auxVec.mpiBcast(0, m_env.subComm());
2485 outOfSupport = !(m_targetDomain->contains(auxVec));
2486 }
while (outOfSupport);
2490 #if 1 // prudencio 2010-08-01
2491 currLogLikelihoodValues[i] = likelihoodSynchronizer.callFunction(&auxVec,NULL,NULL,NULL,NULL,NULL,NULL);
2493 currLogLikelihoodValues[i] = m_likelihoodFunction.lnValue(auxVec,NULL,NULL,NULL,NULL);
2495 currLogTargetValues[i] = m_priorRv.pdf().lnValue(auxVec,NULL,NULL,NULL,NULL) + currLogLikelihoodValues[i];
2499 if (m_env.inter0Rank() >= 0) {
2500 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2501 if (currOptions.m_rawChainComputeStats) {
2510 currChain.computeStatistics(*currOptions.m_rawChainStatisticalOptionsObj,
2527 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2528 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2531 <<
" chain positions"
2547 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2548 if (currOptions.m_filteredChainComputeStats) {
2563 "MLSampling<P_V,P_M>::generateSequence()",
2564 "currChain (first one) has been generated with invalid size");
2567 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2568 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2570 <<
", total level time = " << levelRunTime <<
" seconds"
2577 template <
class P_V,
class P_M>
2581 unsigned int& unifiedRequestedNumSamples)
2584 struct timeval timevalStep;
2585 iRC = gettimeofday(&timevalStep, NULL);
2588 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2589 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2591 <<
", step " << m_currStep
2592 <<
": beginning step 1 of 11"
2600 "MLSampling<P_V,P_M>::generateSequence()",
2601 "failed MPI.Allreduce() for requested num samples in step 1");
2603 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2604 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateSequence()"
2606 <<
", step " << m_currStep
2607 <<
", currOptions->m_rawChainSize = " << currOptions->
m_rawChainSize
2612 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2613 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2615 <<
", step " << m_currStep
2616 <<
", after " << stepRunTime <<
" seconds"
2623 template <
class P_V,
class P_M>
2633 unsigned int& indexOfFirstWeight,
2634 unsigned int& indexOfLastWeight)
2637 struct timeval timevalStep;
2638 iRC = gettimeofday(&timevalStep, NULL);
2641 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2642 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2644 <<
", step " << m_currStep
2645 <<
": beginning step 2 of 11"
2649 prevChain = currChain;
2653 prevLogLikelihoodValues = currLogLikelihoodValues;
2654 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2655 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
2657 <<
", step " << m_currStep
2658 <<
", prevLogLikelihoodValues[0] = " << prevLogLikelihoodValues[0]
2661 prevLogTargetValues = currLogTargetValues;
2663 currLogLikelihoodValues.
clear();
2664 currLogLikelihoodValues.
setName(currOptions->
m_prefix +
"rawLogLikelihood");
2666 currLogTargetValues.
clear();
2669 #if 0 // For debug only
2670 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2671 P_V prevPosition(m_vectorSpace.zeroVector());
2672 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2674 <<
", step " << m_currStep
2679 *m_env.subDisplayFile() <<
" prevChain[" << i
2680 <<
"] = " << prevPosition
2681 <<
", prevLogLikelihoodValues[" << i
2682 <<
"] = " << prevLogLikelihoodValues[i]
2683 <<
", prevLogTargetValues[" << i
2684 <<
"] = " << prevLogTargetValues[i]
2692 unsigned int quantity3 = prevLogLikelihoodValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2693 unsigned int quantity4 = currLogLikelihoodValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2694 unsigned int quantity5 = prevLogTargetValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2695 unsigned int quantity6 = currLogTargetValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2696 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2697 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2699 <<
", step " << m_currStep
2700 <<
": prevChain.unifiedSequenceSize() = " << quantity1
2701 <<
", currChain.unifiedSequenceSize() = " << quantity2
2702 <<
", prevLogLikelihoodValues.unifiedSequenceSize() = " << quantity3
2703 <<
", currLogLikelihoodValues.unifiedSequenceSize() = " << quantity4
2704 <<
", prevLogTargetValues.unifiedSequenceSize() = " << quantity5
2705 <<
", currLogTargetValues.unifiedSequenceSize() = " << quantity6
2711 "MLSampling<P_V,P_M>::generateSequence()",
2712 "different sizes between previous chain and previous sequence of likelihood values");
2716 "MLSampling<P_V,P_M>::generateSequence()",
2717 "different sizes between previous chain and previous sequence of target values");
2720 indexOfFirstWeight = 0;
2721 indexOfLastWeight = indexOfFirstWeight + prevChain.
subSequenceSize()-1;
2724 int r = m_env.inter0Rank();
2726 m_env.inter0Comm().Barrier();
2727 unsigned int auxUint = 0;
2732 "MLSampling<P_V,P_M>::generateSequence()",
2733 "failed MPI.Recv()");
2735 indexOfFirstWeight = auxUint;
2736 indexOfLastWeight = indexOfFirstWeight + prevChain.
subSequenceSize()-1;
2738 if (r < (m_env.inter0Comm().NumProc()-1)) {
2739 auxUint = indexOfLastWeight + 1;
2742 "MLSampling<P_V,P_M>::generateSequence()",
2743 "failed MPI.Send()");
2746 m_env.inter0Comm().Barrier();
2750 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2751 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2753 <<
", step " << m_currStep
2754 <<
", after " << stepRunTime <<
" seconds"
2761 template <
class P_V,
class P_M>
2766 double prevExponent,
2767 double failedExponent,
2768 double& currExponent,
2772 struct timeval timevalStep;
2773 iRC = gettimeofday(&timevalStep, NULL);
2776 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2777 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2779 <<
", step " << m_currStep
2780 <<
", failedExponent = " << failedExponent
2781 <<
": beginning step 3 of 11"
2785 std::vector<double> exponents(2,0.);
2786 exponents[0] = prevExponent;
2789 double nowExponent = 1.;
2790 double nowEffectiveSizeRatio = 0.;
2792 unsigned int nowAttempt = 0;
2793 bool testResult =
false;
2797 double nowUnifiedEvidenceLnFactor = 0.;
2799 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2800 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2802 <<
", step " << m_currStep
2803 <<
", failedExponent = " << failedExponent
2804 <<
": entering loop for computing next exponent"
2805 <<
", with nowAttempt = " << nowAttempt
2809 if (failedExponent > 0.) {
2810 nowExponent = .5*(prevExponent+failedExponent);
2813 if (nowAttempt > 0) {
2814 if (nowEffectiveSizeRatio > meanEffectiveSizeRatio) {
2815 exponents[0] = nowExponent;
2818 exponents[1] = nowExponent;
2820 nowExponent = .5*(exponents[0] + exponents[1]);
2823 double auxExponent = nowExponent;
2824 if (prevExponent != 0.) {
2825 auxExponent /= prevExponent;
2828 double subWeightRatioSum = 0.;
2829 double unifiedWeightRatioSum = 0.;
2832 omegaLnDiffSequence[i] = prevLogLikelihoodValues[i]*auxExponent;
2835 #if 1 // prudenci-2012-07-06
2837 double unifiedOmegaLnMax = omegaLnDiffSequence.
unifiedMaxPlain(m_vectorSpace.numOfProcsForStorage() == 1);
2839 double unifiedOmegaLnMin = 0.;
2840 double unifiedOmegaLnMax = 0.;
2841 omegaLnDiffSequence.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
2843 omegaLnDiffSequence.subSequenceSize(),
2848 omegaLnDiffSequence[i] -= unifiedOmegaLnMax;
2849 weightSequence[i] = exp(omegaLnDiffSequence[i]);
2850 subWeightRatioSum += weightSequence[i];
2851 #if 0 // For debug only
2852 if ((m_currLevel == 1) && (nowAttempt == 6)) {
2853 if (m_env.subDisplayFile() && (m_env.displayVerbosity() >= 99)) {
2854 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2856 <<
", step " << m_currStep
2858 <<
", prevLogLikelihoodValues[i] = " << prevLogLikelihoodValues[i]
2859 <<
", omegaLnDiffSequence[i] = " << omegaLnDiffSequence[i]
2860 <<
", weightSequence[i] = " << weightSequence[i]
2868 "MLSampling<P_V,P_M>::generateSequence()",
2869 "failed MPI.Allreduce() for weight ratio sum");
2871 unsigned int auxQuantity = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2872 nowUnifiedEvidenceLnFactor = log(unifiedWeightRatioSum) + unifiedOmegaLnMax - log(auxQuantity);
2874 double effectiveSampleSize = 0.;
2876 weightSequence[i] /= unifiedWeightRatioSum;
2877 effectiveSampleSize += weightSequence[i]*weightSequence[i];
2888 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2889 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2891 <<
", step " << m_currStep
2892 <<
": nowAttempt = " << nowAttempt
2893 <<
", prevExponent = " << prevExponent
2894 <<
", exponents[0] = " << exponents[0]
2895 <<
", nowExponent = " << nowExponent
2896 <<
", exponents[1] = " << exponents[1]
2897 <<
", subWeightRatioSum = " << subWeightRatioSum
2898 <<
", unifiedWeightRatioSum = " << unifiedWeightRatioSum
2899 <<
", unifiedOmegaLnMax = " << unifiedOmegaLnMax
2900 <<
", weightSequence.unifiedSequenceSize() = " << auxQuantity
2901 <<
", nowUnifiedEvidenceLnFactor = " << nowUnifiedEvidenceLnFactor
2902 <<
", effectiveSampleSize = " << effectiveSampleSize
2906 #if 0 // For debug only
2907 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2908 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2910 <<
", step " << m_currStep
2915 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2916 *m_env.subDisplayFile() <<
" weightSequence[" << i
2917 <<
"] = " << weightSequence[i]
2923 double subQuantity = effectiveSampleSize;
2924 effectiveSampleSize = 0.;
2926 "MLSampling<P_V,P_M>::generateSequence()",
2927 "failed MPI.Allreduce() for effective sample size");
2929 effectiveSampleSize = 1./effectiveSampleSize;
2930 nowEffectiveSizeRatio = effectiveSampleSize/((double) weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1));
2933 "MLSampling<P_V,P_M>::generateSequence()",
2934 "effective sample size ratio cannot be > 1");
2940 if (failedExponent > 0.) {
2945 bool aux2 = (nowExponent == 1. )
2947 (nowEffectiveSizeRatio > meanEffectiveSizeRatio);
2950 (nowEffectiveSizeRatio <= currOptions->m_maxEffectiveSizeRatio);
2951 testResult = aux2 || aux3;
2954 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2955 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2957 <<
", step " << m_currStep
2958 <<
": nowAttempt = " << nowAttempt
2959 <<
", prevExponent = " << prevExponent
2960 <<
", failedExponent = " << failedExponent
2961 <<
", exponents[0] = " << exponents[0]
2962 <<
", nowExponent = " << nowExponent
2963 <<
", exponents[1] = " << exponents[1]
2964 <<
", effectiveSampleSize = " << effectiveSampleSize
2967 <<
", nowEffectiveSizeRatio = " << nowEffectiveSizeRatio
2971 <<
", testResult = " << testResult
2980 "MLSampling<P_V,P_M>::generateSequence(), step 3, nowExponent") ==
false) {
2981 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2982 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2984 <<
", step " << m_currStep
2985 <<
": nowAttempt = " << nowAttempt
2986 <<
", MiscCheck for 'nowExponent' detected a problem"
2995 "MLSampling<P_V,P_M>::generateSequence(), step 3, testResult") ==
false) {
2996 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2997 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2999 <<
", step " << m_currStep
3000 <<
": nowAttempt = " << nowAttempt
3001 <<
", MiscCheck for 'testResult' detected a problem"
3005 }
while (testResult ==
false);
3006 currExponent = nowExponent;
3007 if (failedExponent > 0.) {
3008 m_logEvidenceFactors[m_logEvidenceFactors.size()-1] = nowUnifiedEvidenceLnFactor;
3011 m_logEvidenceFactors.push_back(nowUnifiedEvidenceLnFactor);
3014 unsigned int quantity1 = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3015 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3016 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3018 <<
", step " << m_currStep
3019 <<
": weightSequence.subSequenceSize() = " << weightSequence.
subSequenceSize()
3020 <<
", weightSequence.unifiedSequenceSize() = " << quantity1
3021 <<
", failedExponent = " << failedExponent
3022 <<
", currExponent = " << currExponent
3023 <<
", effective ratio = " << nowEffectiveSizeRatio
3024 <<
", log(evidence factor) = " << m_logEvidenceFactors[m_logEvidenceFactors.size()-1]
3025 <<
", evidence factor = " << exp(m_logEvidenceFactors[m_logEvidenceFactors.size()-1])
3043 "MLSampling<P_V,P_M>::generateSequence(), step 3, logEvidenceFactor") ==
false) {
3044 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3045 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3047 <<
", step " << m_currStep
3048 <<
", failedExponent = " << failedExponent
3049 <<
": nowAttempt = " << nowAttempt
3050 <<
", MiscCheck for 'logEvidenceFactor' detected a problem"
3056 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3057 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3059 <<
", step " << m_currStep
3060 <<
", failedExponent = " << failedExponent
3061 <<
", after " << stepRunTime <<
" seconds"
3068 template <
class P_V,
class P_M>
3073 P_M& unifiedCovMatrix)
3076 struct timeval timevalStep;
3077 iRC = gettimeofday(&timevalStep, NULL);
3080 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3081 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3083 <<
", step " << m_currStep
3084 <<
": beginning step 4 of 11"
3088 P_V auxVec(m_vectorSpace.zeroVector());
3089 P_V subWeightedMeanVec(m_vectorSpace.zeroVector());
3092 subWeightedMeanVec += weightSequence[i]*auxVec;
3096 P_V unifiedWeightedMeanVec(m_vectorSpace.zeroVector());
3097 if (m_env.inter0Rank() >= 0) {
3098 subWeightedMeanVec.mpiAllReduce(
RawValue_MPI_SUM,m_env.inter0Comm(),unifiedWeightedMeanVec);
3101 unifiedWeightedMeanVec = subWeightedMeanVec;
3104 P_V diffVec(m_vectorSpace.zeroVector());
3105 P_M subCovMatrix(m_vectorSpace.zeroVector());
3108 diffVec = auxVec - unifiedWeightedMeanVec;
3109 subCovMatrix += weightSequence[i]*
matrixProduct(diffVec,diffVec);
3112 for (
unsigned int i = 0; i < unifiedCovMatrix.numRowsLocal(); ++i) {
3113 for (
unsigned int j = 0; j < unifiedCovMatrix.numCols(); ++j) {
3114 double localValue = subCovMatrix(i,j);
3115 double sumValue = 0.;
3116 if (m_env.inter0Rank() >= 0) {
3118 "MLSampling<P_V,P_M>::generateSequence()",
3119 "failed MPI.Allreduce() for cov matrix");
3122 sumValue = localValue;
3124 unifiedCovMatrix(i,j) = sumValue;
3128 if (m_numDisabledParameters > 0) {
3129 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
3130 if (m_parameterEnabledStatus[paramId] ==
false) {
3131 for (
unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
3132 unifiedCovMatrix(i,paramId) = 0.;
3134 for (
unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
3135 unifiedCovMatrix(paramId,j) = 0.;
3137 unifiedCovMatrix(paramId,paramId) = 1.;
3142 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3143 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3145 <<
", step " << m_currStep
3146 <<
": unifiedCovMatrix = " << unifiedCovMatrix
3151 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3152 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3154 <<
", step " << m_currStep
3155 <<
", after " << stepRunTime <<
" seconds"
3162 template <
class P_V,
class P_M>
3165 unsigned int unifiedRequestedNumSamples,
3167 std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
3168 std::vector<double>& unifiedWeightStdVectorAtProc0Only)
3171 struct timeval timevalStep;
3172 iRC = gettimeofday(&timevalStep, NULL);
3175 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3176 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3178 <<
", step " << m_currStep
3179 <<
": beginning step 5 of 11"
3183 #if 0 // For debug only
3184 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3185 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3187 <<
", step " << m_currStep
3188 <<
", before weightSequence.getUnifiedContentsAtProc0Only()"
3193 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3194 *m_env.subDisplayFile() <<
", weightSequence[" << i
3195 <<
"] = " << weightSequence[i]
3202 unifiedWeightStdVectorAtProc0Only);
3204 #if 0 // For debug only
3205 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3206 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3208 <<
", step " << m_currStep
3209 <<
", after weightSequence.getUnifiedContentsAtProc0Only()"
3213 for (
unsigned int i = 0; i < unifiedWeightStdVectorAtProc0Only.size(); ++i) {
3214 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3215 *m_env.subDisplayFile() <<
" unifiedWeightStdVectorAtProc0Only[" << i
3216 <<
"] = " << unifiedWeightStdVectorAtProc0Only[i]
3221 sampleIndexes_proc0(unifiedRequestedNumSamples,
3222 unifiedWeightStdVectorAtProc0Only,
3223 unifiedIndexCountersAtProc0Only);
3225 unsigned int auxUnifiedSize = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3226 if (m_env.inter0Rank() == 0) {
3229 "MLSampling<P_V,P_M>::generateSequence()",
3230 "wrong output from sampleIndexesAtProc0() in step 5");
3234 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3235 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3237 <<
", step " << m_currStep
3238 <<
", after " << stepRunTime <<
" seconds"
3245 template <
class P_V,
class P_M>
3249 unsigned int indexOfFirstWeight,
3250 unsigned int indexOfLastWeight,
3251 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
3252 bool& useBalancedChains,
3253 std::vector<ExchangeInfoStruct>& exchangeStdVec)
3256 struct timeval timevalStep;
3257 iRC = gettimeofday(&timevalStep, NULL);
3260 useBalancedChains = decideOnBalancedChains_all(currOptions,
3263 unifiedIndexCountersAtProc0Only,
3267 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3268 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3270 <<
", step " << m_currStep
3271 <<
", after " << stepRunTime <<
" seconds"
3278 template <
class P_V,
class P_M>
3281 bool useBalancedChains,
3282 unsigned int indexOfFirstWeight,
3283 unsigned int indexOfLastWeight,
3284 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
3288 double prevExponent,
3289 double currExponent,
3292 std::vector<ExchangeInfoStruct>& exchangeStdVec,
3296 struct timeval timevalStep;
3297 iRC = gettimeofday(&timevalStep, NULL);
3300 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3301 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3303 <<
", step " << m_currStep
3304 <<
": beginning step 7 of 11"
3308 if (useBalancedChains) {
3309 prepareBalLinkedChains_inter0(currOptions,
3313 prevLogLikelihoodValues,
3314 prevLogTargetValues,
3316 balancedLinkControl);
3319 prepareUnbLinkedChains_inter0(indexOfFirstWeight,
3321 unifiedIndexCountersAtProc0Only,
3322 unbalancedLinkControl);
3325 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3326 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3328 <<
", step " << m_currStep
3329 <<
": balancedLinkControl.balLinkedChains.size() = " << balancedLinkControl.
balLinkedChains.size()
3330 <<
", unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.
unbLinkedChains.size()
3335 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3336 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3338 <<
", step " << m_currStep
3339 <<
", after " << stepRunTime <<
" seconds"
3346 template <
class P_V,
class P_M>
3353 struct timeval timevalStep;
3354 iRC = gettimeofday(&timevalStep, NULL);
3357 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3358 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3360 <<
", step " << m_currStep
3361 <<
": beginning step 8 of 11"
3368 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3369 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3371 <<
", step " << m_currStep
3372 <<
", after " << stepRunTime <<
" seconds"
3379 template <
class P_V,
class P_M>
3383 double prevExponent,
3384 double currExponent,
3387 unsigned int indexOfFirstWeight,
3388 unsigned int indexOfLastWeight,
3389 const std::vector<double>& unifiedWeightStdVectorAtProc0Only,
3394 P_M& unifiedCovMatrix,
3398 struct timeval timevalStep;
3399 iRC = gettimeofday(&timevalStep, NULL);
3403 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3404 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3406 <<
", step " << m_currStep
3407 <<
": skipping step 9 of 11"
3412 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3413 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3415 <<
", step " << m_currStep
3416 <<
": beginning step 9 of 11"
3420 double beforeEta = prevEta;
3421 double beforeRejectionRate = 0.;
3422 bool beforeRejectionRateIsBelowRange =
true;
3424 double nowEta = prevEta;
3425 double nowRejectionRate = 0.;
3426 bool nowRejectionRateIsBelowRange =
true;
3428 std::vector<double> etas(2,0.);
3429 etas[0] = beforeEta;
3432 std::vector<double> rejs(2,0.);
3436 unsigned int nowAttempt = 0;
3437 bool testResult =
false;
3439 bool useMiddlePointLogicForEta =
false;
3440 P_M nowCovMatrix(unifiedCovMatrix);
3441 #if 0 // KAUST, to check
3442 std::vector<double> unifiedWeightStdVectorAtProc0Only(0);
3444 unifiedWeightStdVectorAtProc0Only);
3447 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3448 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3450 <<
", step " << m_currStep
3451 <<
": entering loop for assessing rejection rate"
3452 <<
", with nowAttempt = " << nowAttempt
3453 <<
", nowRejectionRate = " << nowRejectionRate
3456 nowCovMatrix = unifiedCovMatrix;
3458 if (nowRejectionRate < currOptions->m_minRejectionRate) {
3459 nowRejectionRateIsBelowRange =
true;
3462 nowRejectionRateIsBelowRange =
false;
3467 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3468 "nowRejectionRate should be out of the requested range at this point of the logic");
3471 if (m_env.inter0Rank() >= 0) {
3472 if (nowAttempt > 0) {
3473 if (useMiddlePointLogicForEta ==
false) {
3474 if (nowAttempt == 1) {
3477 else if ((beforeRejectionRateIsBelowRange ==
true) &&
3478 (nowRejectionRateIsBelowRange ==
true)) {
3481 else if ((beforeRejectionRateIsBelowRange ==
false) &&
3482 (nowRejectionRateIsBelowRange ==
false)) {
3485 else if ((beforeRejectionRateIsBelowRange ==
true ) &&
3486 (nowRejectionRateIsBelowRange ==
false)) {
3487 useMiddlePointLogicForEta =
true;
3490 etas[0] = std::min(beforeEta,nowEta);
3491 etas[1] = std::max(beforeEta,nowEta);
3493 if (etas[0] == beforeEta) {
3494 rejs[0] = beforeRejectionRate;
3495 rejs[1] = nowRejectionRate;
3498 rejs[0] = nowRejectionRate;
3499 rejs[1] = beforeRejectionRate;
3502 else if ((beforeRejectionRateIsBelowRange ==
false) &&
3503 (nowRejectionRateIsBelowRange ==
true )) {
3504 useMiddlePointLogicForEta =
true;
3507 etas[0] = std::min(beforeEta,nowEta);
3508 etas[1] = std::max(beforeEta,nowEta);
3513 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3514 "before and now range flags are inconsistent");
3519 beforeRejectionRate = nowRejectionRate;
3520 beforeRejectionRateIsBelowRange = nowRejectionRateIsBelowRange;
3521 if (useMiddlePointLogicForEta ==
false) {
3522 if (beforeRejectionRateIsBelowRange) nowEta *= 4.;
3524 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3525 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3527 <<
", step " << m_currStep
3528 <<
": in loop for assessing rejection rate"
3529 <<
", with nowAttempt = " << nowAttempt
3530 <<
", useMiddlePointLogicForEta = false"
3531 <<
", nowEta just updated to value (to be tested) " << nowEta
3536 if (nowRejectionRate > meanRejectionRate) {
3537 if (rejs[0] > meanRejectionRate) {
3547 if (rejs[0] < meanRejectionRate) {
3556 nowEta = .5*(etas[0] + etas[1]);
3557 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3558 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3560 <<
", step " << m_currStep
3561 <<
": in loop for assessing rejection rate"
3562 <<
", with nowAttempt = " << nowAttempt
3563 <<
", useMiddlePointLogicForEta = true"
3564 <<
", nowEta just updated to value (to be tested) " << nowEta
3565 <<
", etas[0] = " << etas[0]
3566 <<
", etas[1] = " << etas[1]
3573 nowCovMatrix *= nowEta;
3577 unsigned int originalSubNumSamples = 1 + (
unsigned int) (doubSubNumSamples);
3578 double auxDouble = (double) originalSubNumSamples;
3579 if ((auxDouble - doubSubNumSamples) < 1.e-8) {
3580 originalSubNumSamples += 1;
3583 if (m_env.inter0Rank() >= 0) {
3584 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3585 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3587 <<
", step " << m_currStep
3588 <<
": in loop for assessing rejection rate"
3589 <<
", about to sample " << originalSubNumSamples <<
" indexes"
3590 <<
", meanRejectionRate = " << meanRejectionRate
3596 std::vector<unsigned int> nowUnifiedIndexCountersAtProc0Only(0);
3597 if (m_env.inter0Rank() >= 0) {
3598 unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3599 sampleIndexes_proc0(tmpUnifiedNumSamples,
3600 unifiedWeightStdVectorAtProc0Only,
3601 nowUnifiedIndexCountersAtProc0Only);
3603 unsigned int auxUnifiedSize = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3604 if (m_env.inter0Rank() == 0) {
3607 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3608 "wrong output from sampleIndexesAtProc0() in step 9");
3611 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3612 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3614 <<
", step " << m_currStep
3615 <<
": in loop for assessing rejection rate"
3616 <<
", about to distribute sampled assessment indexes"
3621 std::vector<ExchangeInfoStruct> exchangeStdVec(0);
3626 bool useBalancedChains = decideOnBalancedChains_all(currOptions,
3629 nowUnifiedIndexCountersAtProc0Only,
3632 if (m_env.inter0Rank() >= 0) {
3633 if (useBalancedChains) {
3634 prepareBalLinkedChains_inter0(currOptions,
3638 prevLogLikelihoodValues,
3639 prevLogTargetValues,
3644 prepareUnbLinkedChains_inter0(indexOfFirstWeight,
3646 nowUnifiedIndexCountersAtProc0Only,
3651 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3652 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3654 <<
", step " << m_currStep
3655 <<
": in loop for assessing rejection rate"
3656 <<
", about to generate assessment chain"
3662 m_options.m_prefix+
"now_chain");
3663 double nowRunTime = 0.;
3664 unsigned int nowRejections = 0;
3669 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3670 bool savedRawChainComputeStats = currOptions->m_rawChainComputeStats;
3677 if (m_env.displayVerbosity() >= 999999) {
3681 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3682 currOptions->m_rawChainComputeStats =
false;
3689 if (useBalancedChains) {
3690 generateBalLinkedChains_all(*currOptions,
3701 generateUnbLinkedChains_all(*currOptions,
3709 prevLogLikelihoodValues,
3710 prevLogTargetValues,
3721 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3722 currOptions->m_rawChainComputeStats = savedRawChainComputeStats;
3728 for (
unsigned int i = 0; i < nowBalLinkControl.
balLinkedChains.size(); ++i) {
3731 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3732 "Initial position pointer in step 9 should not be NULL");
3738 if (m_env.inter0Rank() >= 0) {
3740 unsigned int nowUnifiedRejections = 0;
3742 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3743 "failed MPI.Allreduce() for now rejections");
3745 #if 0 // Round Rock 2009 12 29
3746 unsigned int tmpUnifiedNumSamples = 0;
3748 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3749 "failed MPI.Allreduce() for num samples in step 9");
3752 unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3753 nowRejectionRate = ((double) nowUnifiedRejections) / ((double) tmpUnifiedNumSamples);
3758 (nowRejectionRate <= currOptions->m_maxRejectionRate);
3765 "MLSampling<P_V,P_M>::generateSequence_Step09_all(), step 9, testResult") ==
false) {
3766 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3767 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3769 <<
", step " << m_currStep
3770 <<
": nowAttempt = " << nowAttempt
3771 <<
", MiscCheck for 'testResult' detected a problem"
3778 unsigned int tmpUint = (
unsigned int) testResult;
3780 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3781 "failed MPI.Bcast() for testResult");
3782 testResult = (bool) tmpUint;
3784 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3785 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3787 <<
", step " << m_currStep
3788 <<
": in loop for assessing rejection rate"
3789 <<
", nowAttempt = " << nowAttempt
3790 <<
", beforeEta = " << beforeEta
3791 <<
", etas[0] = " << etas[0]
3792 <<
", nowEta = " << nowEta
3793 <<
", etas[1] = " << etas[1]
3795 <<
", nowRejectionRate = " << nowRejectionRate
3801 if (m_env.inter0Rank() >= 0) {
3806 "MLSampling<P_V,P_M>::generateSequence_Step09_all(), step 9, nowEta") ==
false) {
3807 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3808 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3810 <<
", step " << m_currStep
3811 <<
": nowAttempt = " << nowAttempt
3812 <<
", MiscCheck for 'nowEta' detected a problem"
3817 }
while (testResult ==
false);
3819 if (currEta != 1.) {
3820 unifiedCovMatrix *= currEta;
3821 if (m_numDisabledParameters > 0) {
3822 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
3823 if (m_parameterEnabledStatus[paramId] ==
false) {
3824 for (
unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
3825 unifiedCovMatrix(i,paramId) = 0.;
3827 for (
unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
3828 unifiedCovMatrix(paramId,j) = 0.;
3830 unifiedCovMatrix(paramId,paramId) = 1.;
3836 unsigned int quantity1 = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3837 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3838 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3840 <<
", step " << m_currStep
3841 <<
": weightSequence.subSequenceSize() = " << weightSequence.
subSequenceSize()
3842 <<
", weightSequence.unifiedSequenceSize() = " << quantity1
3843 <<
", currEta = " << currEta
3844 <<
", assessed rejection rate = " << nowRejectionRate
3850 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3851 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3853 <<
", step " << m_currStep
3854 <<
", after " << stepRunTime <<
" seconds"
3861 template <
class P_V,
class P_M>
3865 const P_M& unifiedCovMatrix,
3867 bool useBalancedChains,
3869 unsigned int indexOfFirstWeight,
3871 double prevExponent,
3872 double currExponent,
3877 double& cumulativeRawChainRunTime,
3878 unsigned int& cumulativeRawChainRejections,
3883 struct timeval timevalStep;
3884 iRC = gettimeofday(&timevalStep, NULL);
3887 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3888 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3890 <<
", step " << m_currStep
3891 <<
": beginning step 10 of 11"
3892 <<
", currLogLikelihoodValues = " << currLogLikelihoodValues
3899 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3900 bool savedRawChainComputeStats = currOptions.m_rawChainComputeStats;
3905 if (m_env.displayVerbosity() >= 999999) {
3909 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3910 currOptions.m_rawChainComputeStats =
false;
3915 if (useBalancedChains) {
3916 generateBalLinkedChains_all(currOptions,
3919 balancedLinkControl,
3921 cumulativeRawChainRunTime,
3922 cumulativeRawChainRejections,
3923 currLogLikelihoodValues,
3924 currLogTargetValues);
3927 generateUnbLinkedChains_all(currOptions,
3930 unbalancedLinkControl,
3935 prevLogLikelihoodValues,
3936 prevLogTargetValues,
3938 cumulativeRawChainRunTime,
3939 cumulativeRawChainRejections,
3940 currLogLikelihoodValues,
3941 currLogTargetValues);
3944 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3945 double tmpValue = INFINITY;
3946 if (currLogLikelihoodValues) tmpValue = (*currLogLikelihoodValues)[0];
3947 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
3949 <<
", step " << m_currStep
3950 <<
", after chain generatrion"
3951 <<
", currLogLikelihoodValues[0] = " << tmpValue
3958 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3959 currOptions.m_rawChainComputeStats = savedRawChainComputeStats;
3964 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3965 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3967 <<
", step " << m_currStep
3968 <<
", after " << stepRunTime <<
" seconds"
3975 template <
class P_V,
class P_M>
3979 unsigned int unifiedRequestedNumSamples,
3980 unsigned int cumulativeRawChainRejections,
3984 unsigned int& unifiedNumberOfRejections)
3987 struct timeval timevalStep;
3988 iRC = gettimeofday(&timevalStep, NULL);
3991 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3992 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3994 <<
", step " << m_currStep
3995 <<
": beginning step 11 of 11"
4001 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4002 if (currOptions->m_rawChainComputeStats) {
4010 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
4011 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
4013 <<
", step " << m_currStep
4014 <<
", calling computeStatistics for raw chain"
4015 <<
". Ofstream pointer value = " << filePtrSet.
ofsVar
4016 <<
", statistical options are"
4017 <<
"\n" << *currOptions->m_rawChainStatisticalOptionsObj
4021 currChain.computeStatistics(*currOptions->m_rawChainStatisticalOptionsObj,
4032 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4033 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
4035 <<
", step " << m_currStep
4036 <<
", before calling currLogLikelihoodValues.unifiedWriteContents()"
4037 <<
", currLogLikelihoodValues[0] = " << currLogLikelihoodValues[0]
4056 if (filterSpacing == 0) {
4063 currChain.
filter(filterInitialPos,
4067 currLogLikelihoodValues.
filter(filterInitialPos,
4069 currLogLikelihoodValues.
setName(currOptions->
m_prefix +
"filtLogLikelihood");
4071 currLogTargetValues.
filter(filterInitialPos,
4075 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4076 if (currOptions->m_filteredChainComputeStats) {
4077 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
4078 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
4080 <<
", step " << m_currStep
4081 <<
", calling computeStatistics for filtered chain"
4082 <<
". Ofstream pointer value = " << filePtrSet.
ofsVar
4083 <<
", statistical options are"
4084 <<
"\n" << *currOptions->m_rawChainStatisticalOptionsObj
4089 currChain.computeStatistics(*currOptions->m_filteredChainStatisticalOptionsObj,
4114 unsigned int unifiedGeneratedNumSamples = 0;
4116 "MLSampling<P_V,P_M>::generateSequence()",
4117 "failed MPI.Allreduce() for generated num samples in step 11");
4123 "MLSampling<P_V,P_M>::generateSequence()",
4124 "currChain (linked one) has been generated with invalid size");
4129 "MLSampling<P_V,P_M>::generateSequence()",
4130 "failed MPI.Allreduce() for number of rejections");
4133 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4134 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
4136 <<
", step " << m_currStep
4137 <<
", after " << stepRunTime <<
" seconds"
4145 template<
class P_V,
class P_M>
4151 m_env (priorRv.env()),
4152 m_priorRv (priorRv),
4153 m_likelihoodFunction(likelihoodFunction),
4154 m_vectorSpace (m_priorRv.imageSet().vectorSpace()),
4156 m_numDisabledParameters (0),
4157 m_parameterEnabledStatus(m_vectorSpace.dimLocal(),true),
4158 m_options (m_env,prefix),
4161 m_debugExponent (0.),
4162 m_logEvidenceFactors(0),
4164 m_meanLogLikelihood (0.),
4180 template<
class P_V,
class P_M>
4183 m_numDisabledParameters = 0;
4184 m_parameterEnabledStatus.clear();
4185 if (m_targetDomain)
delete m_targetDomain;
4191 template <
class P_V,
class P_M>
4198 struct timeval timevalRoutineBegin;
4200 iRC = gettimeofday(&timevalRoutineBegin, NULL);
4203 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4204 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::generateSequence()"
4205 <<
", at " << ctime(&timevalRoutineBegin.tv_sec)
4206 <<
", after " << timevalRoutineBegin.tv_sec - m_env.timevalBegin().tv_sec
4207 <<
" seconds from queso environment instatiation..."
4214 double currExponent = 0.;
4215 double currEta = 1.;
4216 unsigned int currUnifiedRequestedNumSamples = 0;
4219 m_options.m_prefix+
"curr_chain");
4227 bool stopAtEndOfLevel =
false;
4228 char levelPrefix[256];
4239 if (m_options.m_restartInput_baseNameForFiles !=
".") {
4240 restartML(currExponent,
4243 currLogLikelihoodValues,
4244 currLogTargetValues);
4246 if (currExponent == 1.) {
4247 if (lastLevelOptions.m_parameterDisabledSet.size() > 0) {
4248 for (std::set<unsigned int>::iterator setIt = lastLevelOptions.m_parameterDisabledSet.begin(); setIt != lastLevelOptions.m_parameterDisabledSet.end(); ++setIt) {
4249 unsigned int paramId = *setIt;
4250 if (paramId < m_vectorSpace.dimLocal()) {
4251 m_numDisabledParameters++;
4252 m_parameterEnabledStatus[paramId] =
false;
4257 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4258 if (lastLevelOptions.m_rawChainComputeStats) {
4260 m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4262 lastLevelOptions.m_dataOutputAllowedSet,
4266 currChain.computeStatistics(*lastLevelOptions.m_rawChainStatisticalOptionsObj,
4276 currChain.
unifiedWriteContents (lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4277 currLogLikelihoodValues.
unifiedWriteContents(lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4278 currLogTargetValues.
unifiedWriteContents (lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4281 if (lastLevelOptions.m_filteredChainGenerate) {
4283 m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4285 lastLevelOptions.m_dataOutputAllowedSet,
4289 unsigned int filterInitialPos = (
unsigned int) (lastLevelOptions.m_filteredChainDiscardedPortion * (
double) currChain.
subSequenceSize());
4290 unsigned int filterSpacing = lastLevelOptions.m_filteredChainLag;
4291 if (filterSpacing == 0) {
4298 currChain.
filter(filterInitialPos,
4300 currChain.
setName(lastLevelOptions.m_prefix +
"filtChain");
4302 currLogLikelihoodValues.
filter(filterInitialPos,
4304 currLogLikelihoodValues.
setName(lastLevelOptions.m_prefix +
"filtLogLikelihood");
4306 currLogTargetValues.
filter(filterInitialPos,
4308 currLogTargetValues.
setName(lastLevelOptions.m_prefix +
"filtLogTarget");
4310 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4311 if (lastLevelOptions.m_filteredChainComputeStats) {
4312 currChain.computeStatistics(*lastLevelOptions.m_filteredChainStatisticalOptionsObj,
4321 currChain.
unifiedWriteContents (lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4322 currLogLikelihoodValues.
unifiedWriteContents(lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4323 currLogTargetValues.
unifiedWriteContents (lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4333 if (currOptions.m_parameterDisabledSet.size() > 0) {
4334 for (std::set<unsigned int>::iterator setIt = currOptions.m_parameterDisabledSet.begin(); setIt != currOptions.m_parameterDisabledSet.end(); ++setIt) {
4335 unsigned int paramId = *setIt;
4336 if (paramId < m_vectorSpace.dimLocal()) {
4337 m_numDisabledParameters++;
4338 m_parameterEnabledStatus[paramId] =
false;
4343 generateSequence_Level0_all(currOptions,
4344 currUnifiedRequestedNumSamples,
4346 currLogLikelihoodValues,
4347 currLogTargetValues);
4349 stopAtEndOfLevel = currOptions.m_stopAtEnd;
4350 bool performCheckpoint = stopAtEndOfLevel;
4351 if (m_options.m_restartOutput_levelPeriod > 0) {
4352 performCheckpoint = performCheckpoint || ( ((m_currLevel + 1) % m_options.m_restartOutput_levelPeriod) == 0 );
4354 if (performCheckpoint) {
4355 checkpointML(currExponent,
4358 currLogLikelihoodValues,
4359 currLogTargetValues);
4365 double minLogLike = -INFINITY;
4366 double maxLogLike = INFINITY;
4371 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4372 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4374 <<
", sub minLogLike = " << minLogLike
4375 <<
", sub maxLogLike = " << maxLogLike
4379 m_env.fullComm().Barrier();
4381 minLogLike = -INFINITY;
4382 maxLogLike = INFINITY;
4383 currLogLikelihoodValues.
unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4388 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4389 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4391 <<
", unified minLogLike = " << minLogLike
4392 <<
", unified maxLogLike = " << maxLogLike
4399 while ((currExponent < 1. ) &&
4400 (stopAtEndOfLevel ==
false)) {
4403 struct timeval timevalLevelBegin;
4405 iRC = gettimeofday(&timevalLevelBegin, NULL);
4407 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4408 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4410 <<
", at " << ctime(&timevalLevelBegin.tv_sec)
4411 <<
", after " << timevalLevelBegin.tv_sec - timevalRoutineBegin.tv_sec
4412 <<
" seconds from entering the routine"
4413 <<
", after " << timevalLevelBegin.tv_sec - m_env.timevalBegin().tv_sec
4414 <<
" seconds from queso environment instatiation"
4419 struct timeval timevalLevel;
4420 iRC = gettimeofday(&timevalLevel, NULL);
4421 double cumulativeRawChainRunTime = 0.;
4422 unsigned int cumulativeRawChainRejections = 0;
4424 bool tryExponentEta =
true;
4425 double failedExponent = 0.;
4426 double failedEta = 0.;
4432 double prevExponent = 0.;
4433 unsigned int indexOfFirstWeight = 0;
4434 unsigned int indexOfLastWeight = 0;
4435 P_M* unifiedCovMatrix = NULL;
4436 bool useBalancedChains =
false;
4442 unsigned int exponentEtaTriedAmount = 0;
4443 while (tryExponentEta) {
4444 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4445 *m_env.subDisplayFile() <<
"In IMLSampling<P_V,P_M>::generateSequence()"
4447 <<
", beginning 'do-while(tryExponentEta)"
4448 <<
": failedExponent = " << failedExponent
4449 <<
", failedEta = " << failedEta
4463 unsigned int paramId = *setIt;
4464 if (paramId < m_vectorSpace.dimLocal()) {
4465 m_numDisabledParameters++;
4466 m_parameterEnabledStatus[paramId] =
false;
4471 if (m_env.inter0Rank() >= 0) {
4472 generateSequence_Step01_inter0(currOptions,
4473 currUnifiedRequestedNumSamples);
4480 prevExponent = currExponent;
4481 double prevEta = currEta;
4482 unsigned int prevUnifiedRequestedNumSamples = currUnifiedRequestedNumSamples;
4485 m_options.m_prefix+
"prev_chain");
4487 indexOfFirstWeight = 0;
4488 indexOfLastWeight = 0;
4491 if (m_env.inter0Rank() >= 0) {
4492 generateSequence_Step02_inter0(currOptions,
4494 currLogLikelihoodValues,
4496 currLogTargetValues,
4498 prevLogLikelihoodValues,
4499 prevLogTargetValues,
4510 if (m_env.inter0Rank() >= 0) {
4511 generateSequence_Step03_inter0(currOptions,
4512 prevLogLikelihoodValues,
4521 "MLSampling<P_V,P_M>::generateSequence()",
4522 "failed MPI.Bcast() for currExponent");
4523 m_debugExponent = currExponent;
4525 if (currExponent == 1.) {
4526 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4527 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4529 <<
", step " << m_currStep
4530 <<
": copying 'last' level options to current options"
4534 currOptions = &lastLevelOptions;
4536 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4537 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4539 <<
", step " << m_currStep
4540 <<
": after copying 'last' level options to current options, the current options are"
4541 <<
"\n" << *currOptions
4545 if (m_env.inter0Rank() >= 0) {
4550 "MLSampling<P_V,P_M>::generateSequence()",
4551 "failed MPI.Allreduce() for requested num samples in step 3");
4559 P_V oneVec(m_vectorSpace.zeroVector());
4562 unifiedCovMatrix = NULL;
4563 if (m_env.inter0Rank() >= 0) {
4564 unifiedCovMatrix = m_vectorSpace.newMatrix();
4567 unifiedCovMatrix =
new P_M(oneVec);
4570 if (m_env.inter0Rank() >= 0) {
4571 generateSequence_Step04_inter0(*prevChain,
4580 std::vector<unsigned int> unifiedIndexCountersAtProc0Only(0);
4581 std::vector<double> unifiedWeightStdVectorAtProc0Only(0);
4582 if (m_env.inter0Rank() >= 0) {
4583 generateSequence_Step05_inter0(currUnifiedRequestedNumSamples,
4585 unifiedIndexCountersAtProc0Only,
4586 unifiedWeightStdVectorAtProc0Only);
4593 useBalancedChains =
false;
4594 std::vector<ExchangeInfoStruct> exchangeStdVec(0);
4596 generateSequence_Step06_all(currOptions,
4599 unifiedIndexCountersAtProc0Only,
4610 if (m_env.inter0Rank() >= 0) {
4611 generateSequence_Step07_inter0(useBalancedChains,
4614 unifiedIndexCountersAtProc0Only,
4615 *unbalancedLinkControl,
4620 prevLogLikelihoodValues,
4621 prevLogTargetValues,
4623 *balancedLinkControl);
4634 m_likelihoodFunction,
4642 generateSequence_Step08_all(*currPdf,
4649 generateSequence_Step09_all(*prevChain,
4652 prevLogLikelihoodValues,
4653 prevLogTargetValues,
4656 unifiedWeightStdVectorAtProc0Only,
4664 tryExponentEta =
false;
4666 (currEta < currOptions->m_minAcceptableEta)) {
4667 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4668 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4670 <<
", preparing to retry ExponentEta"
4671 <<
": currExponent = " << currExponent
4672 <<
", currEta = " << currEta
4675 exponentEtaTriedAmount++;
4676 tryExponentEta =
true;
4677 failedExponent = currExponent;
4678 failedEta = currEta;
4686 delete balancedLinkControl;
4687 balancedLinkControl = NULL;
4688 delete unbalancedLinkControl;
4689 unbalancedLinkControl = NULL;
4691 delete unifiedCovMatrix;
4692 unifiedCovMatrix = NULL;
4694 currExponent = prevExponent;
4696 currUnifiedRequestedNumSamples = prevUnifiedRequestedNumSamples;
4699 currChain = (*prevChain);
4703 currLogLikelihoodValues = prevLogLikelihoodValues;
4704 currLogTargetValues = prevLogTargetValues;
4708 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4709 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4711 <<
", exited 'do-while(tryExponentEta)"
4712 <<
", failedExponent = " << failedExponent
4713 <<
", failedEta = " << failedEta
4722 generateSequence_Step10_all(*currOptions,
4726 *unbalancedLinkControl,
4731 prevLogLikelihoodValues,
4732 prevLogTargetValues,
4733 *balancedLinkControl,
4735 cumulativeRawChainRunTime,
4736 cumulativeRawChainRejections,
4737 &currLogLikelihoodValues,
4738 &currLogTargetValues);
4744 bool performCheckpoint = stopAtEndOfLevel;
4745 if (m_options.m_restartOutput_levelPeriod > 0) {
4746 performCheckpoint = performCheckpoint || ( ((m_currLevel + 1) % m_options.m_restartOutput_levelPeriod) == 0 );
4747 if (currExponent == 1.) {
4748 performCheckpoint =
true;
4751 if (performCheckpoint) {
4752 checkpointML(currExponent,
4755 currLogLikelihoodValues,
4756 currLogTargetValues);
4763 delete unifiedCovMatrix;
4765 for (
unsigned int i = 0; i < balancedLinkControl->
balLinkedChains.size(); ++i) {
4768 "MLSampling<P_V,P_M>::generateSequence()",
4769 "Initial position pointer in step 9 should not be NULL");
4780 unsigned int unifiedNumberOfRejections = 0;
4781 if (m_env.inter0Rank() >= 0) {
4782 generateSequence_Step11_inter0(currOptions,
4783 currUnifiedRequestedNumSamples,
4784 cumulativeRawChainRejections,
4786 currLogLikelihoodValues,
4787 currLogTargetValues,
4788 unifiedNumberOfRejections);
4791 minLogLike = -INFINITY;
4792 maxLogLike = INFINITY;
4797 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4798 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4800 <<
", sub minLogLike = " << minLogLike
4801 <<
", sub maxLogLike = " << maxLogLike
4805 m_env.fullComm().Barrier();
4807 minLogLike = -INFINITY;
4808 maxLogLike = INFINITY;
4809 currLogLikelihoodValues.
unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4814 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4815 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4817 <<
", unified minLogLike = " << minLogLike
4818 <<
", unified maxLogLike = " << maxLogLike
4826 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4827 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4830 <<
" chain positions"
4831 <<
", cumulativeRawChainRunTime = " << cumulativeRawChainRunTime <<
" seconds"
4832 <<
", total level time = " << levelRunTime <<
" seconds"
4833 <<
", cumulativeRawChainRejections = " << cumulativeRawChainRejections
4834 <<
" (" << 100.*((double) cumulativeRawChainRejections)/((double) currOptions->
m_rawChainSize)
4835 <<
"% at this processor)"
4836 <<
" (" << 100.*((double) unifiedNumberOfRejections)/((double) currUnifiedRequestedNumSamples)
4837 <<
"% over all processors)"
4838 <<
", stopAtEndOfLevel = " << stopAtEndOfLevel
4842 if (m_env.inter0Rank() >= 0) {
4843 double minCumulativeRawChainRunTime = 0.;
4845 "MLSampling<P_V,P_M>::generateSequence()",
4846 "failed MPI.Allreduce() for min cumulative raw chain run time");
4848 double maxCumulativeRawChainRunTime = 0.;
4850 "MLSampling<P_V,P_M>::generateSequence()",
4851 "failed MPI.Allreduce() for max cumulative raw chain run time");
4853 double avgCumulativeRawChainRunTime = 0.;
4855 "MLSampling<P_V,P_M>::generateSequence()",
4856 "failed MPI.Allreduce() for sum cumulative raw chain run time");
4857 avgCumulativeRawChainRunTime /= ((double) m_env.inter0Comm().NumProc());
4859 double minLevelRunTime = 0.;
4861 "MLSampling<P_V,P_M>::generateSequence()",
4862 "failed MPI.Allreduce() for min level run time");
4864 double maxLevelRunTime = 0.;
4866 "MLSampling<P_V,P_M>::generateSequence()",
4867 "failed MPI.Allreduce() for max level run time");
4869 double avgLevelRunTime = 0.;
4871 "MLSampling<P_V,P_M>::generateSequence()",
4872 "failed MPI.Allreduce() for sum level run time");
4873 avgLevelRunTime /= ((double) m_env.inter0Comm().NumProc());
4875 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4876 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4878 <<
": min cumul seconds = " << minCumulativeRawChainRunTime
4879 <<
", avg cumul seconds = " << avgCumulativeRawChainRunTime
4880 <<
", max cumul seconds = " << maxCumulativeRawChainRunTime
4881 <<
", min level seconds = " << minLevelRunTime
4882 <<
", avg level seconds = " << avgLevelRunTime
4883 <<
", max level seconds = " << maxLevelRunTime
4888 if (currExponent != 1.)
delete currOptions;
4890 struct timeval timevalLevelEnd;
4892 iRC = gettimeofday(&timevalLevelEnd, NULL);
4894 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4895 *m_env.subDisplayFile() <<
"Getting at the end of level " << m_currLevel+
LEVEL_REF_ID
4896 <<
", as part of a 'while' on levels"
4897 <<
", at " << ctime(&timevalLevelEnd.tv_sec)
4898 <<
", after " << timevalLevelEnd.tv_sec - timevalRoutineBegin.tv_sec
4899 <<
" seconds from entering the routine"
4900 <<
", after " << timevalLevelEnd.tv_sec - m_env.timevalBegin().tv_sec
4901 <<
" seconds from queso environment instatiation"
4915 if (m_env.inter0Rank() >= 0) {
4918 "MLSampling<P_V,P_M>::generateSequence()",
4919 "invalid m_currLevel at the exit of the level loop");
4921 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
4922 m_logEvidence += m_logEvidenceFactors[i];
4925 #if 1 // prudenci-2012-07-06
4926 m_meanLogLikelihood = currLogLikelihoodValues.
unifiedMeanPlain(m_vectorSpace.numOfProcsForStorage() == 1);
4928 m_meanLogLikelihood = currLogLikelihoodValues.
unifiedMeanExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4933 m_eig = m_meanLogLikelihood - m_logEvidence;
4935 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4936 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4937 <<
", log(evidence) = " << m_logEvidence
4938 <<
", evidence = " << exp(m_logEvidence)
4939 <<
", meanLogLikelihood = " << m_meanLogLikelihood
4940 <<
", eig = " << m_eig
4946 "MLSampling<P_V,P_M>::generateSequence()",
4947 "failed MPI.Bcast() for m_logEvidence");
4950 "MLSampling<P_V,P_M>::generateSequence()",
4951 "failed MPI.Bcast() for m_meanLogLikelihood");
4954 "MLSampling<P_V,P_M>::generateSequence()",
4955 "failed MPI.Bcast() for m_eig");
4960 workingChain.
clear();
4962 P_V auxVec(m_vectorSpace.zeroVector());
4964 if (m_env.inter0Rank() >= 0) {
4970 if (workingLogLikelihoodValues) *workingLogLikelihoodValues = currLogLikelihoodValues;
4971 if (workingLogTargetValues ) *workingLogTargetValues = currLogTargetValues;
4973 struct timeval timevalRoutineEnd;
4975 iRC = gettimeofday(&timevalRoutineEnd, NULL);
4977 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4978 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence()"
4979 <<
", at " << ctime(&timevalRoutineEnd.tv_sec)
4980 <<
", after " << timevalRoutineEnd.tv_sec - timevalRoutineBegin.tv_sec
4981 <<
" seconds from entering the routine"
4982 <<
", after " << timevalRoutineEnd.tv_sec - m_env.timevalBegin().tv_sec
4983 <<
" seconds from queso environment instatiation"
4990 template<
class P_V,
class P_M>
4991 std::ostream& operator<<(std::ostream& os, const MLSampling<P_V,P_M>& obj)
4998 template <
class P_V,
class P_M>
5001 return m_logEvidence;
5004 template <
class P_V,
class P_M>
5007 return m_meanLogLikelihood;
5010 template <
class P_V,
class P_M>
bool m_totallyMute
Whether or not to be totally mute (no printout message).
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
#define ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
void clear()
Reset the values and the size of the sequence of vectors.
double m_minRejectionRate
minimum allowed attempted rejection rate at current level
void append(const BaseVectorSequence< V, M > &src, unsigned int initialPos, unsigned int numPos)
Appends the vector src to this vector.
void generateSequence_Step07_inter0(bool useBalancedChains, unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, UnbalancedLinkedChainsPerNodeStruct &unbalancedLinkControl, const MLSamplingLevelOptions *currOptions, const SequenceOfVectors< P_V, P_M > &prevChain, double prevExponent, double currExponent, const ScalarSequence< double > &prevLogLikelihoodValues, const ScalarSequence< double > &prevLogTargetValues, std::vector< ExchangeInfoStruct > &exchangeStdVec, BalancedLinkedChainsPerNodeStruct< P_V > &balancedLinkControl)
Plans for number of linked chains for each node so that all nodes generate the closest possible to th...
void generateSequence_Step08_all(BayesianJointPdf< P_V, P_M > &currPdf, GenericVectorRV< P_V, P_M > &currRv)
Creates a vector RV for current level (Step 08 from ML algorithm).
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
void generateSequence_Step02_inter0(const MLSamplingLevelOptions *currOptions, SequenceOfVectors< P_V, P_M > &currChain, ScalarSequence< double > &currLogLikelihoodValues, ScalarSequence< double > &currLogTargetValues, SequenceOfVectors< P_V, P_M > &prevChain, ScalarSequence< double > &prevLogLikelihoodValues, ScalarSequence< double > &prevLogTargetValues, unsigned int &indexOfFirstWeight, unsigned int &indexOfLastWeight)
Saves chain and corresponding target pdf values from previous level (Step 02 from ML algorithm)...
void prepareUnbLinkedChains_inter0(unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, UnbalancedLinkedChainsPerNodeStruct &unbalancedLinkControl)
std::ofstream * ofsVar
Provides a stream interface to write data to files.
void generateSequence_Step03_inter0(const MLSamplingLevelOptions *currOptions, const ScalarSequence< double > &prevLogLikelihoodValues, double prevExponent, double failedExponent, double &currExponent, ScalarSequence< double > &weightSequence)
Computes currExponent and sequence of weights for current level and update 'm_logEvidenceFactors' (St...
double logEvidence() const
Method to calculate the logarithm of the evidence.
A templated class for a finite distribution.
void getUnifiedContentsAtProc0Only(bool useOnlyInter0Comm, std::vector< T > &outputVec) const
Gets the unified contents of processor of rank equals to 0.
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define RawValue_MPI_DOUBLE
unsigned int sample() const
Samples.
void unifiedMinMaxExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, T &unifiedMinValue, T &unifiedMaxValue) const
Finds the minimum and the maximum values of the unified sequence, considering numPos positions starti...
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
void prepareBalLinkedChains_inter0(const MLSamplingLevelOptions *currOptions, const SequenceOfVectors< P_V, P_M > &prevChain, double prevExponent, double currExponent, const ScalarSequence< double > &prevLogLikelihoodValues, const ScalarSequence< double > &prevLogTargetValues, std::vector< ExchangeInfoStruct > &exchangeStdVec, BalancedLinkedChainsPerNodeStruct< P_V > &balancedLinkControl)
void justBalance_proc0(const MLSamplingLevelOptions *currOptions, std::vector< ExchangeInfoStruct > &exchangeStdVec)
void append(const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
Appends the scalar sequence src to this sequence.
MPI_Status RawType_MPI_Status
std::string m_filteredChainDataOutputFileName
Name of output file for filtered chain.
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
T unifiedMeanExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos) const
Finds the mean value of the unified sequence of numPos positions starting at position initialPos...
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
std::string m_rawChainDataOutputFileType
Type of output file for raw chain.
std::string m_filteredChainDataOutputFileType
Type of output file for filtered chain.
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
#define RawValue_MPI_UNSIGNED
bool m_filteredChainGenerate
Whether or not to generate filtered chain.
#define RawValue_MPI_IN_PLACE
void generateSequence_Step06_all(const MLSamplingLevelOptions *currOptions, unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, bool &useBalancedChains, std::vector< ExchangeInfoStruct > &exchangeStdVec)
Decides on wheter or not to use balanced chains (Step 06 from ML algorithm).
A templated class that represents a Multilevel generator of samples.
unsigned int numberOfPositions
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
std::vector< double > m_initialValuesOfDisabledParameters
double m_maxRejectionRate
maximum allowed attempted rejection rate at current level.
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
double meanLogLikelihood() const
Method to calculate the mean of the logarithm of the likelihood.
unsigned int m_rawChainSize
Size of raw chain.
double m_maxEffectiveSizeRatio
Maximum allowed effective size ratio wrt previous level.
unsigned int numberOfPositions
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
double m_loadBalanceTreshold
Perform load balancing if load unbalancing ratio > threshold.
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
A templated class that represents a Metropolis-Hastings generator of samples.
void mpiExchangePositions_inter0(const SequenceOfVectors< P_V, P_M > &prevChain, double prevExponent, double currExponent, const ScalarSequence< double > &prevLogLikelihoodValues, const ScalarSequence< double > &prevLogTargetValues, const std::vector< ExchangeInfoStruct > &exchangeStdVec, const std::vector< unsigned int > &finalNumChainsPerNode, const std::vector< unsigned int > &finalNumPositionsPerNode, BalancedLinkedChainsPerNodeStruct< P_V > &balancedLinkControl)
bool MiscCheckForSameValueInAllNodes(T &inputValue, double acceptableTreshold, const MpiComm &comm, const char *whereString)
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Struct for handling data input and output from files.
void generateUnbLinkedChains_all(MLSamplingLevelOptions &inputOptions, const P_M &unifiedCovMatrix, const GenericVectorRV< P_V, P_M > &rv, const UnbalancedLinkedChainsPerNodeStruct &unbalancedLinkControl, unsigned int indexOfFirstWeight, const SequenceOfVectors< P_V, P_M > &prevChain, double prevExponent, double currExponent, const ScalarSequence< double > &prevLogLikelihoodValues, const ScalarSequence< double > &prevLogTargetValues, SequenceOfVectors< P_V, P_M > &workingChain, double &cumulativeRunTime, unsigned int &cumulativeRejections, ScalarSequence< double > *currLogLikelihoodValues, ScalarSequence< double > *currLogTargetValues)
unsigned int m_filteredChainLag
Spacing for chain filtering.
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
void generateSequence_Step09_all(const SequenceOfVectors< P_V, P_M > &prevChain, double prevExponent, double currExponent, const ScalarSequence< double > &prevLogLikelihoodValues, const ScalarSequence< double > &prevLogTargetValues, unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< double > &unifiedWeightStdVectorAtProc0Only, const ScalarSequence< double > &weightSequence, double prevEta, const GenericVectorRV< P_V, P_M > &currRv, MLSamplingLevelOptions *currOptions, P_M &unifiedCovMatrix, double &currEta)
Scales the unified covariance matrix until min <= rejection rate <= max (Step 09 from ML algorithm)...
A struct that represents a Metropolis-Hastings sample.
unsigned int numRejections
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
bool m_initialPositionUsePreviousLevelLikelihood
Use previous level likelihood for initial chain position instead of re-computing it from target pdf...
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
std::string m_rawChainDataOutputFileName
Name of output file for raw chain.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
void sampleIndexes_proc0(unsigned int unifiedRequestedNumSamples, const std::vector< double > &unifiedWeightStdVectorAtProc0Only, std::vector< unsigned int > &unifiedIndexCountersAtProc0Only)
void scanOptionsValues()
It scans the option values from the options input file.
void generateSequence_Level0_all(const MLSamplingLevelOptions &currOptions, unsigned int &unifiedRequestedNumSamples, SequenceOfVectors< P_V, P_M > &currChain, ScalarSequence< double > &currLogLikelihoodValues, ScalarSequence< double > &currLogTargetValues)
Generates the sequence at the level 0.
void generateSequence_Step01_inter0(const MLSamplingLevelOptions *currOptions, unsigned int &unifiedRequestedNumSamples)
Reads options for the ML algorithm (Step 01 from ML algorithm).
unsigned int originalIndexOfInitialPosition
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
double m_minEffectiveSizeRatio
Minimum allowed effective size ratio wrt previous level.
double eig() const
Calculates the expected information gain value, EIG.
unsigned int m_drMaxNumExtraStages
'dr' maximum number of extra stages.
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
std::string m_dataOutputFileName
Name of generic output file.
std::set< unsigned int > m_dataOutputAllowedSet
subEnvs that will write to generic output file.
MLSampling(const char *prefix, const BaseVectorRV< P_V, P_M > &priorRv, const BaseScalarFunction< P_V, P_M > &likelihoodFunction)
Constructor.
#define RawValue_MPI_CHAR
bool m_scaleCovMatrix
Whether or not scale proposal covariance matrix.
void restartML(double &currExponent, double &currEta, SequenceOfVectors< P_V, P_M > &currChain, ScalarSequence< double > &currLogLikelihoodValues, ScalarSequence< double > &currLogTargetValues)
Restarts ML algorithm.
Class for handling vector samples (sequence of vectors).
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
void scanOptionsValues(const MLSamplingLevelOptions *defaultOptions)
It scans the option values from the options input file.
std::set< unsigned int > m_parameterDisabledSet
int originalNodeOfInitialPosition
bool m_stopAtEnd
Stop at end of such level.
void generateSequence_Step10_all(MLSamplingLevelOptions &currOptions, const P_M &unifiedCovMatrix, const GenericVectorRV< P_V, P_M > &currRv, bool useBalancedChains, const UnbalancedLinkedChainsPerNodeStruct &unbalancedLinkControl, unsigned int indexOfFirstWeight, const SequenceOfVectors< P_V, P_M > &prevChain, double prevExponent, double currExponent, const ScalarSequence< double > &prevLogLikelihoodValues, const ScalarSequence< double > &prevLogTargetValues, const BalancedLinkedChainsPerNodeStruct< P_V > &balancedLinkControl, SequenceOfVectors< P_V, P_M > &currChain, double &cumulativeRawChainRunTime, unsigned int &cumulativeRawChainRejections, ScalarSequence< double > *currLogLikelihoodValues, ScalarSequence< double > *currLogTargetValues)
Samples the vector RV of current level (Step 10 from ML algorithm).
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
void generateBalLinkedChains_all(MLSamplingLevelOptions &inputOptions, const P_M &unifiedCovMatrix, const GenericVectorRV< P_V, P_M > &rv, const BalancedLinkedChainsPerNodeStruct< P_V > &balancedLinkControl, SequenceOfVectors< P_V, P_M > &workingChain, double &cumulativeRunTime, unsigned int &cumulativeRejections, ScalarSequence< double > *currLogLikelihoodValues, ScalarSequence< double > *currLogTargetValues)
void clear()
Clears the sequence of scalars.
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
void checkpointML(double currExponent, double currEta, const SequenceOfVectors< P_V, P_M > &currChain, const ScalarSequence< double > &currLogLikelihoodValues, const ScalarSequence< double > &currLogTargetValues)
Writes checkpoint data for the ML method.
unsigned int m_amAdaptInterval
'am' adaptation interval.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
void generateSequence_Step05_inter0(unsigned int unifiedRequestedNumSamples, const ScalarSequence< double > &weightSequence, std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, std::vector< double > &unifiedWeightStdVectorAtProc0Only)
Creates unified finite distribution for current level (Step 05 from ML algorithm).
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
bool decideOnBalancedChains_all(const MLSamplingLevelOptions *currOptions, unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, std::vector< ExchangeInfoStruct > &exchangeStdVec)
double m_minAcceptableEta
minimum acceptable eta,
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
double m_filteredChainDiscardedPortion
Initial discarded portion for chain filtering.
std::vector< BalancedLinkedChainControlStruct< P_V > > balLinkedChains
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
unsigned int initialPositionIndexInPreviousChain
double MiscGetEllapsedSeconds(struct timeval *timeval0)
double m_covRejectionRate
c.o.v. for judging attempted rejection rate at current level.
unsigned int m_loadBalanceAlgorithmId
Perform load balancing with chosen algorithm (0 = no balancing).
const BaseEnvironment & m_env
Queso enviroment.
unsigned int unifiedSequenceSize() const
Calculates the size of the unified sequence of vectors.
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
std::vector< UnbalancedLinkedChainControlStruct > unbLinkedChains
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
void generateSequence_Step04_inter0(const SequenceOfVectors< P_V, P_M > &prevChain, const ScalarSequence< double > &weightSequence, P_M &unifiedCovMatrix)
Creates covariance matrix for current level (Step 04 from ML algorithm).
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
MLSamplingOptions m_options
Options for the ML algorithm.
This class provides options for each level of the Multilevel sequence generator if no input file is a...
void generateSequence_Step11_inter0(const MLSamplingLevelOptions *currOptions, unsigned int unifiedRequestedNumSamples, unsigned int cumulativeRawChainRejections, SequenceOfVectors< P_V, P_M > &currChain, ScalarSequence< double > &currLogLikelihoodValues, ScalarSequence< double > &currLogTargetValues, unsigned int &unifiedNumberOfRejections)
Filters chain (Step 11 from ML algorithm).
int finalNodeOfInitialPosition