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;
86 #endif // QUESO_HAS_GLPK
88 template <
class P_V,
class P_M>
91 unsigned int unifiedRequestedNumSamples,
92 const std::vector<double>& unifiedWeightStdVectorAtProc0Only,
93 std::vector<unsigned int>& unifiedIndexCountersAtProc0Only)
95 if (m_env.inter0Rank() != 0)
return;
97 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
98 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::sampleIndexes_proc0()"
100 <<
", step " << m_currStep
101 <<
": unifiedRequestedNumSamples = " << unifiedRequestedNumSamples
102 <<
", unifiedWeightStdVectorAtProc0Only.size() = " << unifiedWeightStdVectorAtProc0Only.size()
106 #if 0 // For debug only
107 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
108 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::sampleIndexes_proc0()"
110 <<
", step " << m_currStep
113 unsigned int numZeros = 0;
114 for (
unsigned int i = 0; i < unifiedWeightStdVectorAtProc0Only.size(); ++i) {
115 *m_env.subDisplayFile() <<
" unifiedWeightStdVectorAtProc0Only[" << i
116 <<
"] = " << unifiedWeightStdVectorAtProc0Only[i]
118 if (unifiedWeightStdVectorAtProc0Only[i] == 0.) numZeros++;
120 *m_env.subDisplayFile() <<
"Number of zeros in unifiedWeightStdVectorAtProc0Only = " << numZeros
125 if (m_env.inter0Rank() == 0) {
126 unsigned int resizeSize = unifiedWeightStdVectorAtProc0Only.size();
127 unifiedIndexCountersAtProc0Only.resize(resizeSize,0);
132 unifiedWeightStdVectorAtProc0Only);
133 for (
unsigned int i = 0; i < unifiedRequestedNumSamples; ++i) {
134 unsigned int index = tmpFd.
sample();
135 unifiedIndexCountersAtProc0Only[index] += 1;
142 template <
class P_V,
class P_M>
146 unsigned int indexOfFirstWeight,
147 unsigned int indexOfLastWeight,
148 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
149 std::vector<ExchangeInfoStruct>& exchangeStdVec)
153 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
154 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
156 <<
", step " << m_currStep
157 <<
": indexOfFirstWeight = " << indexOfFirstWeight
158 <<
", indexOfLastWeight = " << indexOfLastWeight
164 if (m_env.inter0Rank() >= 0) {
165 Np = (
unsigned int) m_env.inter0Comm().NumProc();
167 std::vector<unsigned int> allFirstIndexes(Np,0);
168 std::vector<unsigned int> allLastIndexes(Np,0);
170 if (m_env.inter0Rank() >= 0) {
174 unsigned int auxUInt = indexOfFirstWeight;
175 m_env.inter0Comm().template Gather<unsigned int>(&auxUInt, 1, &allFirstIndexes[0], (int) 1, 0,
176 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
177 "failed MPI.Gather() for first indexes");
179 if (m_env.inter0Rank() == 0) {
180 queso_require_equal_to_msg(allFirstIndexes[0], indexOfFirstWeight,
"failed MPI.Gather() result for first indexes, at proc 0");
183 auxUInt = indexOfLastWeight;
184 m_env.inter0Comm().template Gather<unsigned int>(&auxUInt, 1, &allLastIndexes[0], (int) 1, 0,
185 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
186 "failed MPI.Gather() for last indexes");
188 if (m_env.inter0Rank() == 0) {
197 if (m_env.inter0Rank() == 0) {
198 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
199 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
201 <<
", step " << m_currStep
202 <<
": original distribution of unified indexes in 'inter0Comm' is as follows"
204 for (
unsigned int r = 0; r < Np; ++r) {
205 *m_env.subDisplayFile() <<
" allFirstIndexes[" << r <<
"] = " << allFirstIndexes[r]
206 <<
" allLastIndexes[" << r <<
"] = " << allLastIndexes[r]
210 for (
unsigned int r = 0; r < (Np-1); ++r) {
214 for (
unsigned int r = 0; r < (Np-1); ++r) {
218 std::vector<unsigned int> origNumChainsPerNode (Np,0);
219 std::vector<unsigned int> origNumPositionsPerNode(Np,0);
221 for (
unsigned int i = 0; i < unifiedIndexCountersAtProc0Only.size(); ++i) {
222 if ((allFirstIndexes[r] <= i) &&
223 (i <= allLastIndexes[r] )) {
228 if ((r < (
int) Np ) &&
229 (allFirstIndexes[r] <= i) &&
230 (i <= allLastIndexes[r] )) {
234 std::cerr <<
"In MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
237 <<
", allFirstIndexes[r] = " << allFirstIndexes[r]
238 <<
", allLastIndexes[r] = " << allLastIndexes[r]
243 if (unifiedIndexCountersAtProc0Only[i] != 0) {
244 origNumChainsPerNode [r] += 1;
245 origNumPositionsPerNode[r] += unifiedIndexCountersAtProc0Only[i];
252 exchangeStdVec.push_back(auxInfo);
258 unsigned int totalNumberOfChains = 0;
259 for (
unsigned int r = 0; r < Np; ++r) {
260 totalNumberOfChains += origNumChainsPerNode[r];
262 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
263 *m_env.subDisplayFile() <<
" KEY"
265 <<
", step " << m_currStep
267 <<
", totalNumberOfChains = " << totalNumberOfChains
272 unsigned int origMinPosPerNode = *std::min_element(origNumPositionsPerNode.begin(), origNumPositionsPerNode.end());
273 unsigned int origMaxPosPerNode = *std::max_element(origNumPositionsPerNode.begin(), origNumPositionsPerNode.end());
274 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
275 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
276 *m_env.subDisplayFile() <<
" KEY"
278 <<
", step " << m_currStep
279 <<
", origNumChainsPerNode[" << nodeId <<
"] = " << origNumChainsPerNode[nodeId]
280 <<
", origNumPositionsPerNode[" << nodeId <<
"] = " << origNumPositionsPerNode[nodeId]
284 double origRatioOfPosPerNode = ((double) origMaxPosPerNode ) / ((double) origMinPosPerNode);
285 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
286 *m_env.subDisplayFile() <<
" KEY"
288 <<
", step " << m_currStep
289 <<
", origRatioOfPosPerNode = " << origRatioOfPosPerNode
297 (m_env.numSubEnvironments() > 1 ) &&
298 (Np < totalNumberOfChains ) &&
305 m_env.fullComm().Barrier();
306 unsigned int tmpValue = result;
308 "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
309 "failed MPI.Bcast() for 'result'");
310 if (m_env.inter0Rank() != 0) result = tmpValue;
312 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
313 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
315 <<
", step " << m_currStep
316 <<
": result = " << result
323 template <
class P_V,
class P_M>
332 std::vector<ExchangeInfoStruct>& exchangeStdVec,
335 if (m_env.inter0Rank() < 0)
return;
337 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
338 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
340 <<
", step " << m_currStep
344 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
345 if (m_env.inter0Rank() == 0) {
348 justBalance_proc0(currOptions,
354 #ifdef QUESO_HAS_GLPK
356 solveBIP_proc0(exchangeStdVec);
358 if (m_env.subDisplayFile()) {
359 *m_env.subDisplayFile() <<
"WARNING in MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
361 <<
"' has been requested, but this QUESO library has not been built with 'hdf5'"
362 <<
". Code will therefore process the algorithm id '" << 2
366 if (m_env.subRank() == 0) {
367 std::cerr <<
"WARNING in MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
369 <<
"' has been requested, but this QUESO library has not been built with 'hdf5'"
370 <<
". Code will therefore process the algorithm id '" << 2
374 justBalance_proc0(currOptions,
381 m_env.inter0Comm().Barrier();
386 unsigned int exchangeStdVecSize = exchangeStdVec.size();
388 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
389 "failed MPI.Bcast() for exchangeStdVec size");
390 if (m_env.inter0Rank() > 0) exchangeStdVec.resize(exchangeStdVecSize);
393 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
394 "failed MPI.Bcast() for exchangeStdVec data");
399 std::vector<unsigned int> finalNumChainsPerNode (Np,0);
400 std::vector<unsigned int> finalNumPositionsPerNode(Np,0);
401 unsigned int Nc = exchangeStdVec.size();
402 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
403 unsigned int nodeId = exchangeStdVec[chainId].finalNodeOfInitialPosition;
404 finalNumChainsPerNode [nodeId] += 1;
405 finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
411 unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
412 unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
413 double finalRatioOfPosPerNode = ((double) finalMaxPosPerNode) / ((double)finalMinPosPerNode);
416 std::vector<double> auxBuf(1,0.);
417 double minRatio = 0.;
418 auxBuf[0] = finalRatioOfPosPerNode;
419 m_env.inter0Comm().template Allreduce<double>(&auxBuf[0], &minRatio, (int) auxBuf.size(),
RawValue_MPI_MIN,
420 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
421 "failed MPI.Allreduce() for min");
425 double maxRatio = 0.;
426 auxBuf[0] = finalRatioOfPosPerNode;
427 m_env.inter0Comm().template Allreduce<double>(&auxBuf[0], &maxRatio, (int) auxBuf.size(),
RawValue_MPI_MAX,
428 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
429 "failed MPI.Allreduce() for max");
436 unsigned int finalNumChainsPerNodeSize = finalNumChainsPerNode.size();
438 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
439 "failed MPI.Bcast() for finalNumChainsPerNode size");
440 if (m_env.inter0Rank() > 0) finalNumChainsPerNode.resize(finalNumChainsPerNodeSize);
442 m_env.inter0Comm().Bcast((
void *) &finalNumChainsPerNode[0], (
int) finalNumChainsPerNodeSize,
RawValue_MPI_UNSIGNED, 0,
443 "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
444 "failed MPI.Bcast() for finalNumChainsPerNode data");
450 mpiExchangePositions_inter0(prevChain,
453 prevLogLikelihoodValues,
456 finalNumChainsPerNode,
457 finalNumPositionsPerNode,
458 balancedLinkControl);
463 template <
class P_V,
class P_M>
466 unsigned int indexOfFirstWeight,
467 unsigned int indexOfLastWeight,
468 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
471 if (m_env.inter0Rank() < 0)
return;
473 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
474 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
476 <<
", step " << m_currStep
477 <<
": indexOfFirstWeight = " << indexOfFirstWeight
478 <<
", indexOfLastWeight = " << indexOfLastWeight
482 unsigned int subNumSamples = 0;
483 std::vector<unsigned int> unifiedIndexCountersAtAllProcs(0);
486 unsigned int resizeSize = unifiedIndexCountersAtProc0Only.size();
488 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
489 "failed MPI.Bcast() for resizeSize");
490 unifiedIndexCountersAtAllProcs.resize(resizeSize,0);
492 if (m_env.inter0Rank() == 0) unifiedIndexCountersAtAllProcs = unifiedIndexCountersAtProc0Only;
495 m_env.inter0Comm().Bcast((
void *) &unifiedIndexCountersAtAllProcs[0], (
int) unifiedIndexCountersAtAllProcs.size(),
RawValue_MPI_UNSIGNED, 0,
496 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
497 "failed MPI.Bcast() for unified index counters");
498 #if 0 // Use allgatherv ??? for subNumSamples instead
499 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
500 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
502 <<
", step " << m_currStep
505 for (
int r = 0; r < m_env.inter0Comm().NumProc(); ++r) {
506 *m_env.subDisplayFile() <<
" unifiedIndexCountersAtAllProcs[" << r <<
"] = " << unifiedIndexCountersAtAllProcs[r]
518 queso_require_less_msg(indexOfFirstWeight, unifiedIndexCountersAtAllProcs.size(),
"invalid indexOfFirstWeight");
519 queso_require_less_msg(indexOfLastWeight, unifiedIndexCountersAtAllProcs.size(),
"invalid indexOfLastWeight");
521 for (
unsigned int i = indexOfFirstWeight; i <= indexOfLastWeight; ++i) {
522 subNumSamples += unifiedIndexCountersAtAllProcs[i];
525 std::vector<unsigned int> auxBuf(1,0);
527 unsigned int minModifiedSubNumSamples = 0;
528 auxBuf[0] = subNumSamples;
529 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &minModifiedSubNumSamples, (int) auxBuf.size(),
RawValue_MPI_MIN,
530 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
531 "failed MPI.Allreduce() for min");
533 unsigned int maxModifiedSubNumSamples = 0;
534 auxBuf[0] = subNumSamples;
535 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &maxModifiedSubNumSamples, (int) auxBuf.size(),
RawValue_MPI_MAX,
536 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
537 "failed MPI.Allreduce() for max");
539 unsigned int sumModifiedSubNumSamples = 0;
540 auxBuf[0] = subNumSamples;
541 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &sumModifiedSubNumSamples, (int) auxBuf.size(),
RawValue_MPI_SUM,
542 "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
543 "failed MPI.Allreduce() for sum");
550 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
551 *m_env.subDisplayFile() <<
"KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
553 <<
", step " << m_currStep
554 <<
": subNumSamples = " << subNumSamples
555 <<
", unifiedIndexCountersAtAllProcs.size() = " << unifiedIndexCountersAtAllProcs.size()
557 *m_env.subDisplayFile() <<
"KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
559 <<
", step " << m_currStep
560 <<
": minModifiedSubNumSamples = " << minModifiedSubNumSamples
561 <<
", avgModifiedSubNumSamples = " << ((double) sumModifiedSubNumSamples)/((double) m_env.inter0Comm().NumProc())
562 <<
", maxModifiedSubNumSamples = " << maxModifiedSubNumSamples
566 unsigned int numberOfPositionsToGuaranteeForNode = subNumSamples;
567 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
568 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
570 <<
", step " << m_currStep
571 <<
": numberOfPositionsToGuaranteeForNode = " << numberOfPositionsToGuaranteeForNode
574 for (
unsigned int i = indexOfFirstWeight; i <= indexOfLastWeight; ++i) {
576 while (unifiedIndexCountersAtAllProcs[i] != 0) {
577 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 30)) {
578 *m_env.subDisplayFile() <<
", numberOfPositionsToGuaranteeForNode = " << numberOfPositionsToGuaranteeForNode
579 <<
", unifiedIndexCountersAtAllProcs[" << i
580 <<
"] = " << unifiedIndexCountersAtAllProcs[i]
583 if (unifiedIndexCountersAtAllProcs[i] < numberOfPositionsToGuaranteeForNode) {
589 numberOfPositionsToGuaranteeForNode -= unifiedIndexCountersAtAllProcs[i];
590 unifiedIndexCountersAtAllProcs[i] = 0;
592 else if ((unifiedIndexCountersAtAllProcs[i] == numberOfPositionsToGuaranteeForNode) &&
593 (unifiedIndexCountersAtAllProcs[i] > 0 )) {
600 unifiedIndexCountersAtAllProcs[i] -= numberOfPositionsToGuaranteeForNode;
601 numberOfPositionsToGuaranteeForNode = 0;
603 else if ((unifiedIndexCountersAtAllProcs[i] == numberOfPositionsToGuaranteeForNode) &&
604 (unifiedIndexCountersAtAllProcs[i] == 0 )) {
612 queso_require_equal_to_msg(numberOfPositionsToGuaranteeForNode, 0,
"numberOfPositionsToGuaranteeForNode exited loop with wrong value");
615 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
616 *m_env.subDisplayFile() <<
"KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
618 <<
", step " << m_currStep
619 <<
": unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.
unbLinkedChains.size()
626 template <
class P_V,
class P_M>
630 const P_M& unifiedCovMatrix,
634 double& cumulativeRunTime,
635 unsigned int& cumulativeRejections,
639 m_env.fullComm().Barrier();
641 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
642 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
643 <<
": balancedLinkControl.balLinkedChains.size() = " << balancedLinkControl.
balLinkedChains.size()
647 P_V auxInitialPosition(m_vectorSpace.zeroVector());
648 double auxInitialLogPrior;
649 double auxInitialLogLikelihood;
651 unsigned int chainIdMax = 0;
652 if (m_env.inter0Rank() >= 0) {
657 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
658 "failed MPI.Bcast() for chainIdMax");
660 struct timeval timevalEntering;
662 iRC = gettimeofday(&timevalEntering, NULL);
665 if (m_env.inter0Rank() >= 0) {
666 unsigned int numberOfPositions = 0;
667 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
668 numberOfPositions += balancedLinkControl.
balLinkedChains[chainId].numberOfPositions;
671 std::vector<unsigned int> auxBuf(1,0);
673 unsigned int minNumberOfPositions = 0;
674 auxBuf[0] = numberOfPositions;
675 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &minNumberOfPositions, (int) auxBuf.size(),
RawValue_MPI_MIN,
676 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
677 "failed MPI.Allreduce() for min");
679 unsigned int maxNumberOfPositions = 0;
680 auxBuf[0] = numberOfPositions;
681 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &maxNumberOfPositions, (int) auxBuf.size(),
RawValue_MPI_MAX,
682 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
683 "failed MPI.Allreduce() for max");
685 unsigned int sumNumberOfPositions = 0;
686 auxBuf[0] = numberOfPositions;
687 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &sumNumberOfPositions, (int) auxBuf.size(),
RawValue_MPI_SUM,
688 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
689 "failed MPI.Allreduce() for sum");
691 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
692 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
694 <<
", step " << m_currStep
695 <<
": chainIdMax = " << chainIdMax
696 <<
", numberOfPositions = " << numberOfPositions
697 <<
", at " << ctime(&timevalEntering.tv_sec)
699 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
701 <<
", step " << m_currStep
702 <<
": minNumberOfPositions = " << minNumberOfPositions
703 <<
", avgNumberOfPositions = " << ((double) sumNumberOfPositions)/((double) m_env.inter0Comm().NumProc())
704 <<
", maxNumberOfPositions = " << maxNumberOfPositions
710 if ((m_debugExponent == 1.) &&
711 (m_currStep == 10)) {
714 unsigned int cumulativeNumPositions = 0;
715 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
716 unsigned int tmpChainSize = 0;
717 if (m_env.inter0Rank() >= 0) {
719 auxInitialPosition = *(balancedLinkControl.
balLinkedChains[chainId].initialPosition);
720 auxInitialLogPrior = balancedLinkControl.
balLinkedChains[chainId].initialLogPrior;
721 auxInitialLogLikelihood = balancedLinkControl.
balLinkedChains[chainId].initialLogLikelihood;
722 tmpChainSize = balancedLinkControl.
balLinkedChains[chainId].numberOfPositions+1;
723 if ((m_env.subDisplayFile() ) &&
724 (m_env.displayVerbosity() >= 3)) {
725 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
727 <<
", step " << m_currStep
728 <<
", chainId = " << chainId
729 <<
" < " << chainIdMax
730 <<
": begin generating " << tmpChainSize
731 <<
" chain positions"
735 auxInitialPosition.mpiBcast(0, m_env.subComm());
737 #if 0 // For debug only
738 for (
int r = 0; r < m_env.subComm().NumProc(); ++r) {
739 if (r == m_env.subComm().MyPID()) {
740 std::cout <<
"Vector 'auxInitialPosition at rank " << r
741 <<
" has contents " << auxInitialPosition
744 m_env.subComm().Barrier();
751 "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
752 "failed MPI.Bcast() for tmpChainSize");
757 m_options.m_prefix+
"tmp_chain");
764 "MLSamplingClass<P_V,P_M>::generateBalLinkedChains_all()",
765 "failed MPI.Bcast() for auxInitialLogPrior");
767 "MLSamplingClass<P_V,P_M>::generateBalLinkedChains_all()",
768 "failed MPI.Bcast() for auxInitialLogLikelihood");
774 auxInitialLogLikelihood,
777 &tmpLogLikelihoodValues,
778 &tmpLogTargetValues);
787 &tmpLogLikelihoodValues,
788 &tmpLogTargetValues);
792 cumulativeRunTime += mcRawInfo.
runTime;
795 if (m_env.inter0Rank() >= 0) {
796 if (m_env.exceptionalCircumstance()) {
797 if ((m_env.subDisplayFile() ) &&
798 (m_env.displayVerbosity() >= 0)) {
799 P_V tmpVec(m_vectorSpace.zeroVector());
800 for (
unsigned int i = 0; i < tmpLogLikelihoodValues.
subSequenceSize(); ++i) {
802 *m_env.subDisplayFile() <<
"DEBUG finalChain[" << cumulativeNumPositions+i <<
"] "
803 <<
"= tmpChain[" << i <<
"] = " << tmpVec
804 <<
", tmpLogLikelihoodValues[" << i <<
"] = " << tmpLogLikelihoodValues[i]
805 <<
", tmpLogTargetValues[" << i <<
"] = " << tmpLogTargetValues[i]
811 cumulativeNumPositions += tmpChainSize;
812 if (cumulativeNumPositions > 100) m_env.setExceptionalCircumstance(
false);
814 if ((m_env.subDisplayFile() ) &&
815 (m_env.displayVerbosity() >= 3)) {
816 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
818 <<
", step " << m_currStep
819 <<
", chainId = " << chainId
820 <<
" < " << chainIdMax
822 <<
" chain positions"
828 if (currLogLikelihoodValues) {
829 currLogLikelihoodValues->
append(tmpLogLikelihoodValues,1,tmpLogLikelihoodValues.
subSequenceSize()-1);
830 if ((m_env.subDisplayFile() ) &&
831 (m_env.displayVerbosity() >= 99) &&
833 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
835 <<
", step " << m_currStep
836 <<
", chainId = " << chainId
837 <<
", tmpLogLikelihoodValues.subSequenceSize() = " << tmpLogLikelihoodValues.
subSequenceSize()
838 <<
", tmpLogLikelihoodValues[0] = " << tmpLogLikelihoodValues[0]
839 <<
", tmpLogLikelihoodValues[1] = " << tmpLogLikelihoodValues[1]
840 <<
", currLogLikelihoodValues[0] = " << (*currLogLikelihoodValues)[0]
844 if (currLogTargetValues) {
853 struct timeval timevalBarrier;
854 iRC = gettimeofday(&timevalBarrier, NULL);
857 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
858 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
860 <<
", step " << m_currStep
861 <<
": ended chain loop after " << loopTime <<
" seconds"
862 <<
", calling fullComm().Barrier() at " << ctime(&timevalBarrier.tv_sec)
866 m_env.fullComm().Barrier();
868 struct timeval timevalLeaving;
869 iRC = gettimeofday(&timevalLeaving, NULL);
872 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
873 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
875 <<
", step " << m_currStep
876 <<
": after " << barrierTime <<
" seconds in fullComm().Barrier()"
877 <<
", at " << ctime(&timevalLeaving.tv_sec)
884 template <
class P_V,
class P_M>
888 const P_M& unifiedCovMatrix,
891 unsigned int indexOfFirstWeight,
898 double& cumulativeRunTime,
899 unsigned int& cumulativeRejections,
903 m_env.fullComm().Barrier();
905 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
906 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
907 <<
": unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.
unbLinkedChains.size()
908 <<
", indexOfFirstWeight = " << indexOfFirstWeight
912 P_V auxInitialPosition(m_vectorSpace.zeroVector());
913 double auxInitialLogPrior;
914 double auxInitialLogLikelihood;
916 unsigned int chainIdMax = 0;
917 if (m_env.inter0Rank() >= 0) {
922 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
923 "failed MPI.Bcast() for chainIdMax");
925 struct timeval timevalEntering;
927 iRC = gettimeofday(&timevalEntering, NULL);
930 if (m_env.inter0Rank() >= 0) {
931 unsigned int numberOfPositions = 0;
932 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
933 numberOfPositions += unbalancedLinkControl.
unbLinkedChains[chainId].numberOfPositions;
936 std::vector<unsigned int> auxBuf(1,0);
938 unsigned int minNumberOfPositions = 0;
939 auxBuf[0] = numberOfPositions;
940 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &minNumberOfPositions, (int) auxBuf.size(),
RawValue_MPI_MIN,
941 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
942 "failed MPI.Allreduce() for min");
944 unsigned int maxNumberOfPositions = 0;
945 auxBuf[0] = numberOfPositions;
946 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &maxNumberOfPositions, (int) auxBuf.size(),
RawValue_MPI_MAX,
947 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
948 "failed MPI.Allreduce() for max");
950 unsigned int sumNumberOfPositions = 0;
951 auxBuf[0] = numberOfPositions;
952 m_env.inter0Comm().template Allreduce<unsigned int>(&auxBuf[0], &sumNumberOfPositions, (int) auxBuf.size(),
RawValue_MPI_SUM,
953 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
954 "failed MPI.Allreduce() for sum");
956 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
957 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
959 <<
", step " << m_currStep
960 <<
": chainIdMax = " << chainIdMax
961 <<
", numberOfPositions = " << numberOfPositions
962 <<
", at " << ctime(&timevalEntering.tv_sec)
964 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
966 <<
", step " << m_currStep
967 <<
": minNumberOfPositions = " << minNumberOfPositions
968 <<
", avgNumberOfPositions = " << ((double) sumNumberOfPositions)/((double) m_env.inter0Comm().NumProc())
969 <<
", maxNumberOfPositions = " << maxNumberOfPositions
973 if ((m_debugExponent == 1.) &&
974 (m_currStep == 10)) {
977 double expRatio = currExponent;
978 if (prevExponent > 0.0) {
979 expRatio /= prevExponent;
981 unsigned int cumulativeNumPositions = 0;
982 for (
unsigned int chainId = 0; chainId < chainIdMax; ++chainId) {
983 unsigned int tmpChainSize = 0;
984 if (m_env.inter0Rank() >= 0) {
985 unsigned int auxIndex = unbalancedLinkControl.
unbLinkedChains[chainId].initialPositionIndexInPreviousChain - indexOfFirstWeight;
987 auxInitialLogPrior = prevLogTargetValues[auxIndex] - prevLogLikelihoodValues[auxIndex];
988 auxInitialLogLikelihood = expRatio * prevLogLikelihoodValues[auxIndex];
989 tmpChainSize = unbalancedLinkControl.
unbLinkedChains[chainId].numberOfPositions+1;
990 if ((m_env.subDisplayFile() ) &&
991 (m_env.displayVerbosity() >= 3)) {
992 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
994 <<
", step " << m_currStep
995 <<
", chainId = " << chainId
996 <<
" < " << chainIdMax
997 <<
": begin generating " << tmpChainSize
998 <<
" chain positions"
1002 auxInitialPosition.mpiBcast(0, m_env.subComm());
1004 #if 0 // For debug only
1005 for (
int r = 0; r < m_env.subComm().NumProc(); ++r) {
1006 if (r == m_env.subComm().MyPID()) {
1007 std::cout <<
"Vector 'auxInitialPosition at rank " << r
1008 <<
" has contents " << auxInitialPosition
1011 m_env.subComm().Barrier();
1018 "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
1019 "failed MPI.Bcast() for tmpChainSize");
1024 m_options.m_prefix+
"tmp_chain");
1032 "MLSamplingClass<P_V,P_M>::generateUnbLinkedChains_all()",
1033 "failed MPI.Bcast() for auxInitialLogPrior");
1035 "MLSamplingClass<P_V,P_M>::generateUnbLinkedChains_all",
1036 "failed MPI.Bcast() for auxInitialLogLikelihood");
1041 auxInitialLogLikelihood,
1044 &tmpLogLikelihoodValues,
1045 &tmpLogTargetValues);
1054 &tmpLogLikelihoodValues,
1055 &tmpLogTargetValues);
1059 cumulativeRunTime += mcRawInfo.
runTime;
1062 if (m_env.inter0Rank() >= 0) {
1063 if (m_env.exceptionalCircumstance()) {
1064 if ((m_env.subDisplayFile() ) &&
1065 (m_env.displayVerbosity() >= 0)) {
1066 P_V tmpVec(m_vectorSpace.zeroVector());
1067 for (
unsigned int i = 0; i < tmpLogLikelihoodValues.
subSequenceSize(); ++i) {
1069 *m_env.subDisplayFile() <<
"DEBUG finalChain[" << cumulativeNumPositions+i <<
"] "
1070 <<
"= tmpChain[" << i <<
"] = " << tmpVec
1071 <<
", tmpLogLikelihoodValues[" << i <<
"] = " << tmpLogLikelihoodValues[i]
1072 <<
", tmpLogTargetValues[" << i <<
"] = " << tmpLogTargetValues[i]
1078 cumulativeNumPositions += tmpChainSize;
1079 if (cumulativeNumPositions > 100) m_env.setExceptionalCircumstance(
false);
1081 if ((m_env.subDisplayFile() ) &&
1082 (m_env.displayVerbosity() >= 3)) {
1083 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1085 <<
", step " << m_currStep
1086 <<
", chainId = " << chainId
1087 <<
" < " << chainIdMax
1089 <<
" chain positions"
1095 if (currLogLikelihoodValues) {
1096 currLogLikelihoodValues->
append(tmpLogLikelihoodValues,1,tmpLogLikelihoodValues.
subSequenceSize()-1);
1097 if ((m_env.subDisplayFile() ) &&
1098 (m_env.displayVerbosity() >= 99) &&
1100 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1102 <<
", step " << m_currStep
1103 <<
", chainId = " << chainId
1104 <<
", tmpLogLikelihoodValues.subSequenceSize() = " << tmpLogLikelihoodValues.
subSequenceSize()
1105 <<
", tmpLogLikelihoodValues[0] = " << tmpLogLikelihoodValues[0]
1106 <<
", tmpLogLikelihoodValues[1] = " << tmpLogLikelihoodValues[1]
1107 <<
", currLogLikelihoodValues[0] = " << (*currLogLikelihoodValues)[0]
1111 if (currLogTargetValues) {
1117 struct timeval timevalBarrier;
1118 iRC = gettimeofday(&timevalBarrier, NULL);
1121 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1122 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1124 <<
", step " << m_currStep
1125 <<
": ended chain loop after " << loopTime <<
" seconds"
1126 <<
", calling fullComm().Barrier() at " << ctime(&timevalBarrier.tv_sec)
1130 m_env.fullComm().Barrier();
1132 struct timeval timevalLeaving;
1133 iRC = gettimeofday(&timevalLeaving, NULL);
1136 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1137 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1139 <<
", step " << m_currStep
1140 <<
": after " << barrierTime <<
" seconds in fullComm().Barrier()"
1141 <<
", at " << ctime(&timevalLeaving.tv_sec)
1148 #ifdef QUESO_HAS_GLPK
1149 template <
class P_V,
class P_M>
1152 std::vector<ExchangeInfoStruct>& exchangeStdVec)
1154 if (m_env.inter0Rank() != 0)
return;
1157 struct timeval timevalBIP;
1158 iRC = gettimeofday(&timevalBIP, NULL);
1161 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
1162 unsigned int Nc = exchangeStdVec.size();
1164 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1165 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::solveBIP_proc0()"
1167 <<
", step " << m_currStep
1177 lp = glp_create_prob();
1178 glp_set_prob_name(lp,
"sample");
1183 unsigned int m = Nc+Np-1;
1184 unsigned int n = Nc*Np;
1185 unsigned int ne = Nc*Np + 2*Nc*(Np -1);
1187 glp_add_rows(lp, m);
1188 for (
int i = 1; i <= (int) Nc; ++i) {
1189 glp_set_row_bnds(lp, i, GLP_FX, 1.0, 1.0);
1190 glp_set_row_name(lp, i,
"");
1192 for (
int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1193 glp_set_row_bnds(lp, i, GLP_UP, 0.0, 0.0);
1194 glp_set_row_name(lp, i,
"");
1197 glp_add_cols(lp, n);
1198 for (
int j = 1; j <= (int) n; ++j) {
1200 glp_set_col_kind(lp, j, GLP_IV);
1201 glp_set_col_bnds(lp, j, GLP_DB, 0.0, 1.0);
1202 glp_set_col_name(lp, j,
"");
1205 glp_set_obj_dir(lp, GLP_MIN);
1206 for (
int chainId = 0; chainId <= (int) (Nc-1); ++chainId) {
1207 glp_set_obj_coef(lp, (chainId*Np)+1, exchangeStdVec[chainId].numberOfPositions);
1210 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1211 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1213 <<
", step " << m_currStep
1214 <<
": finished setting BIP rows and cols"
1224 std::vector<int > iVec(ne+1,0);
1225 std::vector<int > jVec(ne+1,0);
1226 std::vector<double> aVec(ne+1,0.);
1228 for (
int i = 1; i <= (int) Nc; ++i) {
1229 for (
int j = 1; j <= (int) Np; ++j) {
1231 jVec[coefId] = (i-1)*Np + j;
1236 for (
int i = 1; i <= (int) (Np-1); ++i) {
1237 for (
int j = 1; j <= (int) Nc; ++j) {
1238 iVec[coefId] = Nc+i;
1239 jVec[coefId] = (j-1)*Np + i;
1240 aVec[coefId] = -((double) exchangeStdVec[j-1].numberOfPositions);
1243 iVec[coefId] = Nc+i;
1244 jVec[coefId] = (j-1)*Np + i + 1;
1245 aVec[coefId] = exchangeStdVec[j-1].numberOfPositions;
1249 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1250 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1252 <<
", step " << m_currStep
1253 <<
": finished setting BIP constraint matrix"
1255 <<
", coefId = " << coefId
1261 glp_load_matrix(lp, ne, &iVec[0], &jVec[0], &aVec[0]);
1263 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1264 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1266 <<
", step " << m_currStep
1267 <<
": finished loading BIP constraint matrix"
1268 <<
", glp_get_num_rows(lp) = " << glp_get_num_rows(lp)
1269 <<
", glp_get_num_cols(lp) = " << glp_get_num_cols(lp)
1270 <<
", glp_get_num_nz(lp) = " << glp_get_num_nz(lp)
1271 <<
", glp_get_num_int(lp) = " << glp_get_num_int(lp)
1272 <<
", glp_get_num_bin(lp) = " << glp_get_num_bin(lp)
1292 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1293 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1294 int j = chainId*Np + nodeId + 1;
1296 glp_set_col_stat(lp, j, GLP_BS);
1299 glp_set_col_stat(lp, j, GLP_BS);
1304 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1305 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1306 int j = chainId*Np + nodeId + 1;
1307 int initialState = glp_mip_col_val(lp, j);
1317 for (
int i = 1; i <= (int) Nc; ++i) {
1318 glp_set_row_stat(lp, i, GLP_NS);
1320 for (
int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1321 glp_set_row_stat(lp, i, GLP_BS);
1331 BIP_routine_struct BIP_routine_info;
1332 BIP_routine_info.env = &m_env;
1333 BIP_routine_info.currLevel = m_currLevel;
1335 glp_iocp BIP_params;
1336 glp_init_iocp(&BIP_params);
1337 BIP_params.presolve = GLP_ON;
1342 int BIP_rc = glp_intopt(lp, &BIP_params);
1344 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1345 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1347 <<
", step " << m_currStep
1348 <<
": finished solving BIP"
1349 <<
", BIP_rc = " << BIP_rc
1358 int BIP_Status = glp_mip_status(lp);
1359 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1360 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1362 <<
", step " << m_currStep
1363 <<
": BIP_Status = " << BIP_Status
1367 switch (BIP_Status) {
1370 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1371 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1373 <<
", step " << m_currStep
1374 <<
": BIP solution is optimal"
1381 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1382 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1384 <<
", step " << m_currStep
1385 <<
": BIP solution is guaranteed to be 'only' feasible"
1395 for (
int i = 1; i <= (int) Nc; ++i) {
1398 for (
int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1405 std::vector<unsigned int> finalNumChainsPerNode (Np,0);
1406 std::vector<unsigned int> finalNumPositionsPerNode(Np,0);
1407 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1408 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1409 int j = chainId*Np + nodeId + 1;
1410 if (glp_mip_col_val(lp, j) == 0) {
1413 else if (glp_mip_col_val(lp, j) == 1) {
1415 exchangeStdVec[chainId].finalNodeOfInitialPosition = nodeId;
1416 finalNumChainsPerNode [nodeId] += 1;
1417 finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
1425 unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1426 unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1428 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1429 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::solveBIP_proc0()"
1431 <<
", step " << m_currStep
1432 <<
": finished preparing output information"
1439 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1440 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1442 <<
", step " << m_currStep
1443 <<
": solution gives the following redistribution"
1445 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1446 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1448 <<
", step " << m_currStep
1449 <<
", finalNumChainsPerNode[" << nodeId <<
"] = " << finalNumChainsPerNode[nodeId]
1450 <<
", finalNumPositionsPerNode[" << nodeId <<
"] = " << finalNumPositionsPerNode[nodeId]
1453 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1455 <<
", step " << m_currStep
1456 <<
", finalRatioOfPosPerNode = " << ((double) finalMaxPosPerNode) / ((double)finalMinPosPerNode)
1465 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
1466 queso_require_greater_equal_msg(finalNumPositionsPerNode[nodeId-1], finalNumPositionsPerNode[nodeId],
"Next node should have a number of positions equal or less than the current node");
1469 for (
int i = (
int) (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1470 unsigned int nodeId = i - Nc;
1471 int diff = ((int) finalNumPositionsPerNode[nodeId]) - ((int) finalNumPositionsPerNode[nodeId-1]);
1478 glp_delete_prob(lp);
1481 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1482 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::solveBIP_proc0()"
1484 <<
", step " << m_currStep
1485 <<
", after " << bipRunTime <<
" seconds"
1491 #endif // QUESO_HAS_GLPK
1493 template <
class P_V,
class P_M>
1497 std::vector<ExchangeInfoStruct>& exchangeStdVec)
1499 if (m_env.inter0Rank() != 0)
return;
1502 struct timeval timevalBal;
1503 iRC = gettimeofday(&timevalBal, NULL);
1506 unsigned int Np = m_env.numSubEnvironments();
1507 unsigned int Nc = exchangeStdVec.size();
1509 std::vector<ExchangeInfoStruct> currExchangeStdVec(Nc);
1510 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1511 currExchangeStdVec[chainId] = exchangeStdVec[chainId];
1512 currExchangeStdVec[chainId].finalNodeOfInitialPosition = currExchangeStdVec[chainId].originalNodeOfInitialPosition;
1518 unsigned int iterIdMax = 0;
1519 std::vector<unsigned int> currNumChainsPerNode (Np,0);
1520 std::vector<unsigned int> currNumPositionsPerNode(Np,0);
1521 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1522 unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1523 currNumChainsPerNode [nodeId] += 1;
1524 currNumPositionsPerNode[nodeId] += currExchangeStdVec[chainId].numberOfPositions;
1525 iterIdMax += currExchangeStdVec[chainId].numberOfPositions;
1527 unsigned int currMinPosPerNode = *std::min_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1528 unsigned int currMaxPosPerNode = *std::max_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1529 double currRatioOfPosPerNode = ((double) currMaxPosPerNode ) / ((double) currMinPosPerNode);
1536 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1537 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1539 <<
", step " << m_currStep
1540 <<
", iter " << iterId
1541 <<
", currRatioOfPosPerNode = " << currRatioOfPosPerNode
1543 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1544 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1546 <<
", step " << m_currStep
1547 <<
", iter " << iterId
1548 <<
", currNumChainsPerNode[" << nodeId <<
"] = " << currNumChainsPerNode[nodeId]
1549 <<
", currNumPositionsPerNode[" << nodeId <<
"] = " << currNumPositionsPerNode[nodeId]
1554 std::vector<std::vector<double> > vectorOfChainSizesPerNode(Np);
1555 while ((iterId < (
int) iterIdMax ) &&
1562 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1563 vectorOfChainSizesPerNode[nodeId].clear();
1565 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1566 unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1567 vectorOfChainSizesPerNode[nodeId].push_back(currExchangeStdVec[chainId].numberOfPositions);
1570 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1571 std::sort(vectorOfChainSizesPerNode[nodeId].begin(), vectorOfChainSizesPerNode[nodeId].end());
1572 queso_require_equal_to_msg(vectorOfChainSizesPerNode[nodeId].size(), currNumChainsPerNode[nodeId],
"inconsistent number of chains in node");
1578 unsigned int currBiggestAmountOfPositionsPerNode = currNumPositionsPerNode[0];
1579 unsigned int currSmallestAmountOfPositionsPerNode = currNumPositionsPerNode[0];
1580 unsigned int currNodeWithMostPositions = 0;
1581 unsigned int currNodeWithLeastPositions = 0;
1582 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1583 if (currNumPositionsPerNode[nodeId] > currBiggestAmountOfPositionsPerNode) {
1584 currBiggestAmountOfPositionsPerNode = currNumPositionsPerNode[nodeId];
1585 currNodeWithMostPositions = nodeId;
1587 if (currNumPositionsPerNode[nodeId] < currSmallestAmountOfPositionsPerNode) {
1588 currSmallestAmountOfPositionsPerNode = currNumPositionsPerNode[nodeId];
1589 currNodeWithLeastPositions = nodeId;
1593 queso_require_equal_to_msg(currMinPosPerNode, currNumPositionsPerNode[currNodeWithLeastPositions],
"inconsistent currMinPosPerNode");
1595 queso_require_equal_to_msg(currMaxPosPerNode, currNumPositionsPerNode[currNodeWithMostPositions],
"inconsistent currMaxPosPerNode");
1597 unsigned int numberOfPositionsToMove = vectorOfChainSizesPerNode[currNodeWithMostPositions][0];
1599 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1600 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1602 <<
", step " << m_currStep
1603 <<
", iter " << iterId
1604 <<
", before update"
1605 <<
", node w/ most pos is "
1606 << currNodeWithMostPositions <<
"(cs=" << currNumChainsPerNode[currNodeWithMostPositions ] <<
", ps=" << currNumPositionsPerNode[currNodeWithMostPositions ] <<
")"
1607 <<
", node w/ least pos is "
1608 << currNodeWithLeastPositions <<
"(cs=" << currNumChainsPerNode[currNodeWithLeastPositions] <<
", ps=" << currNumPositionsPerNode[currNodeWithLeastPositions] <<
")"
1609 <<
", number of pos to move = " << numberOfPositionsToMove
1616 std::vector<ExchangeInfoStruct> newExchangeStdVec(Nc);
1617 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1618 newExchangeStdVec[chainId] = currExchangeStdVec[chainId];
1621 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1622 if ((newExchangeStdVec[chainId].finalNodeOfInitialPosition == (
int) currNodeWithMostPositions) &&
1623 (newExchangeStdVec[chainId].numberOfPositions == numberOfPositionsToMove )) {
1624 newExchangeStdVec[chainId].finalNodeOfInitialPosition = currNodeWithLeastPositions;
1632 std::vector<unsigned int> newNumChainsPerNode (Np,0);
1633 std::vector<unsigned int> newNumPositionsPerNode(Np,0);
1634 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1635 unsigned int nodeId = newExchangeStdVec[chainId].finalNodeOfInitialPosition;
1636 newNumChainsPerNode [nodeId] += 1;
1637 newNumPositionsPerNode[nodeId] += newExchangeStdVec[chainId].numberOfPositions;
1640 unsigned int newBiggestAmountOfPositionsPerNode = newNumPositionsPerNode[0];
1641 unsigned int newSmallestAmountOfPositionsPerNode = newNumPositionsPerNode[0];
1642 unsigned int newNodeWithMostPositions = 0;
1643 unsigned int newNodeWithLeastPositions = 0;
1644 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1645 if (newNumPositionsPerNode[nodeId] > newBiggestAmountOfPositionsPerNode) {
1646 newBiggestAmountOfPositionsPerNode = newNumPositionsPerNode[nodeId];
1647 newNodeWithMostPositions = nodeId;
1649 if (newNumPositionsPerNode[nodeId] < newSmallestAmountOfPositionsPerNode) {
1650 newSmallestAmountOfPositionsPerNode = newNumPositionsPerNode[nodeId];
1651 newNodeWithLeastPositions = nodeId;
1655 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1656 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1658 <<
", step " << m_currStep
1659 <<
", iter " << iterId
1661 <<
", node w/ most pos is "
1662 << newNodeWithMostPositions <<
"(cs=" << newNumChainsPerNode[newNodeWithMostPositions ] <<
", ps=" << newNumPositionsPerNode[newNodeWithMostPositions ] <<
")"
1663 <<
", node w/ least pos is "
1664 << newNodeWithLeastPositions <<
"(cs=" << newNumChainsPerNode[newNodeWithLeastPositions] <<
", ps=" << newNumPositionsPerNode[newNodeWithLeastPositions] <<
")"
1668 unsigned int newMinPosPerNode = *std::min_element(newNumPositionsPerNode.begin(), newNumPositionsPerNode.end());
1669 unsigned int newMaxPosPerNode = *std::max_element(newNumPositionsPerNode.begin(), newNumPositionsPerNode.end());
1670 double newRatioOfPosPerNode = ((double) newMaxPosPerNode ) / ((double) newMinPosPerNode);
1672 queso_require_equal_to_msg(newMinPosPerNode, newNumPositionsPerNode[newNodeWithLeastPositions],
"inconsistent newMinPosPerNode");
1674 queso_require_equal_to_msg(newMaxPosPerNode, newNumPositionsPerNode[newNodeWithMostPositions],
"inconsistent newMaxPosPerNode");
1676 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1677 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::justBalance_proc0()"
1679 <<
", step " << m_currStep
1680 <<
", iter " << iterId
1681 <<
", newMaxPosPerNode = " << newMaxPosPerNode
1682 <<
", newMinPosPerNode = " << newMinPosPerNode
1683 <<
", newRatioOfPosPerNode = " << newRatioOfPosPerNode
1685 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1686 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1688 <<
", step " << m_currStep
1689 <<
", iter " << iterId
1690 <<
", newNumChainsPerNode[" << nodeId <<
"] = " << newNumChainsPerNode [nodeId]
1691 <<
", newNumPositionsPerNode[" << nodeId <<
"] = " << newNumPositionsPerNode[nodeId]
1699 if (newRatioOfPosPerNode > currRatioOfPosPerNode) {
1706 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1707 currNumChainsPerNode [nodeId] = 0;
1708 currNumPositionsPerNode[nodeId] = 0;
1710 currRatioOfPosPerNode = newRatioOfPosPerNode;
1711 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1712 currExchangeStdVec[chainId] = newExchangeStdVec[chainId];
1713 unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1714 currNumChainsPerNode [nodeId] += 1;
1715 currNumPositionsPerNode[nodeId] += currExchangeStdVec[chainId].numberOfPositions;
1717 currMinPosPerNode = *std::min_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1718 currMaxPosPerNode = *std::max_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1719 currRatioOfPosPerNode = ((double) currMaxPosPerNode ) / ((double) currMinPosPerNode);
1725 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1726 exchangeStdVec[chainId].finalNodeOfInitialPosition = currExchangeStdVec[chainId].finalNodeOfInitialPosition;
1732 std::vector<unsigned int> finalNumChainsPerNode (Np,0);
1733 std::vector<unsigned int> finalNumPositionsPerNode(Np,0);
1734 for (
unsigned int chainId = 0; chainId < Nc; ++chainId) {
1735 unsigned int nodeId = exchangeStdVec[chainId].finalNodeOfInitialPosition;
1736 finalNumChainsPerNode [nodeId] += 1;
1737 finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
1739 unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1740 unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1741 double finalRatioOfPosPerNode = ((double) finalMaxPosPerNode ) / ((double) finalMinPosPerNode);
1743 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1744 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1746 <<
", step " << m_currStep
1747 <<
": solution gives the following redistribution"
1749 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1750 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1752 <<
", step " << m_currStep
1753 <<
", finalNumChainsPerNode[" << nodeId <<
"] = " << finalNumChainsPerNode[nodeId]
1754 <<
", finalNumPositionsPerNode[" << nodeId <<
"] = " << finalNumPositionsPerNode[nodeId]
1757 *m_env.subDisplayFile() <<
" KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1759 <<
", step " << m_currStep
1760 <<
", finalRatioOfPosPerNode = " << finalRatioOfPosPerNode
1768 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1769 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::justBalance_proc0()"
1771 <<
", step " << m_currStep
1772 <<
", iterId = " << iterId
1773 <<
", currRatioOfPosPerNode = " << currRatioOfPosPerNode
1774 <<
", after " << balRunTime <<
" seconds"
1781 template <
class P_V,
class P_M>
1785 double prevExponent,
1786 double currExponent,
1789 const std::vector<ExchangeInfoStruct>& exchangeStdVec,
1790 const std::vector<unsigned int>& finalNumChainsPerNode,
1791 const std::vector<unsigned int>& finalNumPositionsPerNode,
1794 if (m_env.inter0Rank() < 0) {
1798 unsigned int Np = (
unsigned int) m_env.inter0Comm().NumProc();
1799 unsigned int Nc = exchangeStdVec.size();
1801 double expRatio = currExponent;
1802 if (prevExponent > 0.0) {
1803 expRatio /= prevExponent;
1806 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1807 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1809 <<
", step " << m_currStep
1820 for (
unsigned int r = 0; r < Np; ++r) {
1824 unsigned int numberOfInitialPositionsNodeRAlreadyHas = 0;
1825 std::vector<unsigned int> numberOfInitialPositionsNodeRHasToReceiveFromNode(Np,0);
1826 std::vector<unsigned int> indexesOfInitialPositionsNodeRHasToReceiveFromMe(0);
1828 unsigned int sumOfChainLenghtsNodeRAlreadyHas = 0;
1829 std::vector<unsigned int> chainLenghtsNodeRHasToInherit(0);
1831 for (
unsigned int i = 0; i < Nc; ++i) {
1832 if (exchangeStdVec[i].finalNodeOfInitialPosition == (
int) r) {
1833 if (exchangeStdVec[i].originalNodeOfInitialPosition == (
int) r) {
1834 numberOfInitialPositionsNodeRAlreadyHas++;
1835 sumOfChainLenghtsNodeRAlreadyHas += exchangeStdVec[i].numberOfPositions;
1838 numberOfInitialPositionsNodeRHasToReceiveFromNode[exchangeStdVec[i].originalNodeOfInitialPosition]++;
1839 chainLenghtsNodeRHasToInherit.push_back(exchangeStdVec[i].numberOfPositions);
1840 if (m_env.inter0Rank() == exchangeStdVec[i].originalNodeOfInitialPosition) {
1841 indexesOfInitialPositionsNodeRHasToReceiveFromMe.push_back(exchangeStdVec[i].originalIndexOfInitialPosition);
1847 unsigned int totalNumberOfInitialPositionsNodeRHasToReceive = 0;
1848 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1849 totalNumberOfInitialPositionsNodeRHasToReceive += numberOfInitialPositionsNodeRHasToReceiveFromNode[nodeId];
1852 unsigned int totalNumberOfChainLenghtsNodeRHasToInherit = chainLenghtsNodeRHasToInherit.size();
1853 unsigned int totalSumOfChainLenghtsNodeRHasToInherit = 0;
1854 for (
unsigned int i = 0; i < totalNumberOfChainLenghtsNodeRHasToInherit; ++i) {
1855 totalSumOfChainLenghtsNodeRHasToInherit += chainLenghtsNodeRHasToInherit[i];
1861 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1862 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1864 <<
", step " << m_currStep
1866 <<
", finalNumChainsPerNode[r] = " << finalNumChainsPerNode[r]
1867 <<
", totalNumberOfInitialPositionsNodeRHasToReceive = " << totalNumberOfInitialPositionsNodeRHasToReceive
1868 <<
", numberOfInitialPositionsNodeRAlreadyHas = " << numberOfInitialPositionsNodeRAlreadyHas
1870 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1872 <<
", step " << m_currStep
1874 <<
", finalNumPositionsPerNode[r] = " << finalNumPositionsPerNode[r]
1875 <<
", totalSumOfChainLenghtsNodeRHasToInherit = " << totalSumOfChainLenghtsNodeRHasToInherit
1876 <<
", sumOfChainLenghtsNodeRAlreadyHas = " << sumOfChainLenghtsNodeRAlreadyHas
1883 queso_require_equal_to_msg(indexesOfInitialPositionsNodeRHasToReceiveFromMe.size(), numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()],
"inconsistent number of initial positions to send to node 'r'");
1885 queso_require_equal_to_msg(finalNumChainsPerNode[r], (totalNumberOfInitialPositionsNodeRHasToReceive + numberOfInitialPositionsNodeRAlreadyHas),
"inconsistent number of chains in node 'r'");
1887 queso_require_equal_to_msg(finalNumPositionsPerNode[r], (totalSumOfChainLenghtsNodeRHasToInherit + sumOfChainLenghtsNodeRAlreadyHas),
"inconsistent sum of chain lenghts in node 'r'");
1889 queso_require_equal_to_msg(totalNumberOfInitialPositionsNodeRHasToReceive, totalNumberOfChainLenghtsNodeRHasToInherit,
"inconsistent on total number of initial positions to receive in node 'r'");
1892 indexesOfInitialPositionsNodeRHasToReceiveFromMe.resize(numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()]);
1893 chainLenghtsNodeRHasToInherit.resize (totalSumOfChainLenghtsNodeRHasToInherit);
1898 unsigned int dimSize = m_vectorSpace.dimLocal();
1899 unsigned int nValuesPerInitialPosition = dimSize + 2;
1900 P_V auxInitialPosition(m_vectorSpace.zeroVector());
1901 std::vector<double> sendbuf(0);
1902 unsigned int sendcnt = 0;
1903 if (m_env.inter0Rank() != (int) r) {
1904 sendcnt = numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()] * nValuesPerInitialPosition;
1905 sendbuf.resize(sendcnt);
1906 for (
unsigned int i = 0; i < numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()]; ++i) {
1907 unsigned int auxIndex = indexesOfInitialPositionsNodeRHasToReceiveFromMe[i];
1909 for (
unsigned int j = 0; j < dimSize; ++j) {
1910 sendbuf[i*nValuesPerInitialPosition + j] = auxInitialPosition[j];
1912 sendbuf[i*nValuesPerInitialPosition + dimSize] = prevLogLikelihoodValues[auxIndex];
1913 sendbuf[i*nValuesPerInitialPosition + dimSize + 1] = prevLogTargetValues[auxIndex];
1917 std::vector<double> recvbuf(0);
1918 std::vector<int> recvcnts(Np,0);
1919 if (m_env.inter0Rank() == (int) r) {
1920 recvbuf.resize(totalNumberOfInitialPositionsNodeRHasToReceive * nValuesPerInitialPosition);
1921 for (
unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1922 recvcnts[nodeId] = numberOfInitialPositionsNodeRHasToReceiveFromNode[nodeId] * nValuesPerInitialPosition;
1926 std::vector<int> displs(Np,0);
1927 for (
unsigned int nodeId = 1; nodeId < Np; ++nodeId) {
1928 displs[nodeId] = displs[nodeId-1] + recvcnts[nodeId-1];
1932 if (m_env.inter0Rank() == r) {
1934 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0(1)",
1935 "failed MPI.Gatherv()");
1939 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0(2)",
1940 "failed MPI.Gatherv()");
1943 m_env.inter0Comm().template Gatherv<double>(&sendbuf[0], (int) sendcnt,
1944 &recvbuf[0], (
int *) &recvcnts[0], (
int *) &displs[0],
1946 "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
1947 "failed MPI.Gatherv()");
1959 if (m_env.inter0Rank() == (int) r) {
1961 unsigned int auxIndex = 0;
1963 for (
unsigned int i = 0; i < Nc; ++i) {
1964 if ((exchangeStdVec[i].finalNodeOfInitialPosition == (
int) r) &&
1965 (exchangeStdVec[i].originalNodeOfInitialPosition == (
int) r)) {
1966 unsigned int originalIndex = exchangeStdVec[i].originalIndexOfInitialPosition;
1968 balancedLinkControl.
balLinkedChains[auxIndex].initialPosition =
new P_V(auxInitialPosition);
1969 balancedLinkControl.
balLinkedChains[auxIndex].initialLogPrior = prevLogTargetValues[originalIndex] - prevLogLikelihoodValues[originalIndex];
1970 balancedLinkControl.
balLinkedChains[auxIndex].initialLogLikelihood = expRatio*prevLogLikelihoodValues[originalIndex];
1971 balancedLinkControl.
balLinkedChains[auxIndex].numberOfPositions = exchangeStdVec[i].numberOfPositions;
1976 for (
unsigned int i = 0; i < totalNumberOfInitialPositionsNodeRHasToReceive; ++i) {
1977 for (
unsigned int j = 0; j < dimSize; ++j) {
1978 auxInitialPosition[j] = recvbuf[i*nValuesPerInitialPosition + j];
1980 balancedLinkControl.
balLinkedChains[auxIndex].initialPosition =
new P_V(auxInitialPosition);
1981 double prevLogLikelihood = recvbuf[i*nValuesPerInitialPosition + dimSize];
1982 double prevLogTarget = recvbuf[i*nValuesPerInitialPosition + dimSize + 1];
1983 balancedLinkControl.
balLinkedChains[auxIndex].initialLogPrior = prevLogTarget - prevLogLikelihood;
1984 balancedLinkControl.
balLinkedChains[auxIndex].initialLogLikelihood = expRatio*prevLogLikelihood;
1985 balancedLinkControl.
balLinkedChains[auxIndex].numberOfPositions = chainLenghtsNodeRHasToInherit[i];
1990 m_env.inter0Comm().Barrier();
1993 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1994 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1996 <<
", step " << m_currStep
2004 template <
class P_V,
class P_M>
2007 double currExponent,
2013 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2014 *m_env.subDisplayFile() <<
"\n CHECKPOINTING initiating at level " << m_currLevel
2015 <<
"\n" << std::endl;
2022 unsigned int quantity2 = currLogLikelihoodValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2023 unsigned int quantity3 = currLogTargetValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2024 if (m_env.inter0Rank() >= 0) {
2030 if (m_env.fullRank() == 0) {
2031 std::ofstream* ofsVar =
new std::ofstream((m_options.m_restartOutput_baseNameForFiles +
"Control.txt").c_str(),
2032 std::ofstream::out | std::ofstream::trunc);
2033 *ofsVar << m_currLevel << std::endl
2034 << m_vectorSpace.dimGlobal() << std::endl
2035 << currExponent << std::endl
2036 << currEta << std::endl
2037 << quantity1 << std::endl;
2038 unsigned int savedPrecision = ofsVar->precision();
2039 ofsVar->precision(16);
2040 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2041 *ofsVar << m_logEvidenceFactors[i] << std::endl;
2043 ofsVar->precision(savedPrecision);
2044 *ofsVar <<
"COMPLETE" << std::endl;
2048 m_env.fullComm().Barrier();
2053 char levelSufix[256];
2056 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2057 *m_env.subDisplayFile() <<
"\n CHECKPOINTING chain at level " << m_currLevel
2058 <<
"\n" << std::endl;
2060 currChain.
unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles +
"Chain_l" + levelSufix,
2061 m_options.m_restartOutput_fileType);
2062 m_env.fullComm().Barrier();
2064 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2065 *m_env.subDisplayFile() <<
"\n CHECKPOINTING like at level " << m_currLevel
2066 <<
"\n" << std::endl;
2068 currLogLikelihoodValues.
unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles +
"LogLike_l" + levelSufix,
2069 m_options.m_restartOutput_fileType);
2070 m_env.fullComm().Barrier();
2072 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2073 *m_env.subDisplayFile() <<
"\n CHECKPOINTING target at level " << m_currLevel
2074 <<
"\n" << std::endl;
2076 currLogTargetValues.
unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles +
"LogTarget_l" + levelSufix,
2077 m_options.m_restartOutput_fileType);
2078 m_env.fullComm().Barrier();
2083 if (m_env.fullRank() == 0) {
2084 std::ofstream* ofsVar =
new std::ofstream((m_options.m_restartOutput_baseNameForFiles +
"Control_l" + levelSufix +
".txt").c_str(),
2085 std::ofstream::out | std::ofstream::trunc);
2086 *ofsVar << m_currLevel << std::endl
2087 << m_vectorSpace.dimGlobal() << std::endl
2088 << currExponent << std::endl
2089 << currEta << std::endl
2090 << quantity1 << std::endl;
2091 unsigned int savedPrecision = ofsVar->precision();
2092 ofsVar->precision(16);
2093 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2094 *ofsVar << m_logEvidenceFactors[i] << std::endl;
2096 ofsVar->precision(savedPrecision);
2097 *ofsVar <<
"COMPLETE" << std::endl;
2101 m_env.fullComm().Barrier();
2103 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2104 *m_env.subDisplayFile() <<
"\n CHECKPOINTING done at level " << m_currLevel
2105 <<
"\n" << std::endl;
2111 template <
class P_V,
class P_M>
2114 double& currExponent,
2120 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2121 *m_env.subDisplayFile() <<
"\n RESTARTING initiating at level " << m_currLevel
2122 <<
"\n" << std::endl;
2128 unsigned int vectorSpaceDim = 0;
2129 unsigned int quantity1 = 0;
2130 std::string checkingString(
"");
2131 if (m_env.fullRank() == 0) {
2132 std::ifstream* ifsVar =
new std::ifstream((m_options.m_restartInput_baseNameForFiles +
"Control.txt").c_str(),
2138 unsigned int numLines = std::count(std::istreambuf_iterator<char>(*ifsVar),
2139 std::istreambuf_iterator<char>(),
2141 ifsVar->seekg(0,std::ios_base::beg);
2142 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2143 *m_env.subDisplayFile() <<
"Restart input file has " << numLines
2151 *ifsVar >> m_currLevel;
2154 m_logEvidenceFactors.clear();
2155 m_logEvidenceFactors.resize(m_currLevel,0.);
2156 *ifsVar >> vectorSpaceDim
2160 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2161 *ifsVar >> m_logEvidenceFactors[i];
2163 *ifsVar >> checkingString;
2166 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2167 *m_env.subDisplayFile() <<
"Restart input file has the following information:"
2168 <<
"\n m_currLevel = " << m_currLevel
2169 <<
"\n vectorSpaceDim = " << vectorSpaceDim
2170 <<
"\n currExponent = " << currExponent
2171 <<
"\n currEta = " << currEta
2172 <<
"\n quantity1 = " << quantity1;
2173 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2174 *m_env.subDisplayFile() <<
"\n [" << i <<
"] = " << m_logEvidenceFactors[i];
2176 *m_env.subDisplayFile() << std::endl;
2179 #if 0 // For debug only
2180 std::string tmpString;
2181 for (
unsigned int i = 0; i < 2; ++i) {
2182 *ifsVar >> tmpString;
2183 std::cout <<
"Just read '" << tmpString <<
"'" << std::endl;
2185 while ((lineId < numLines) && (ifsVar->eof() ==
false)) {
2187 ifsVar->ignore(maxCharsPerLine,
'\n');
2192 m_env.fullComm().Barrier();
2197 unsigned int tmpUint = (
unsigned int) m_currLevel;
2199 "MLSampling<P_V,P_M>::restartML()",
2200 "failed MPI.Bcast() for m_currLevel");
2201 if (m_env.fullRank() != 0) {
2202 m_currLevel = tmpUint;
2209 if (m_env.fullRank() == 0) {
2210 tmpData[0] = vectorSpaceDim;
2211 tmpData[1] = currExponent;
2212 tmpData[2] = currEta;
2213 tmpData[3] = quantity1;
2214 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2215 tmpData[4+i] = m_logEvidenceFactors[i];
2219 m_logEvidenceFactors.clear();
2220 m_logEvidenceFactors.resize(m_currLevel,0.);
2222 m_env.fullComm().Bcast((
void *) &tmpData[0], (
int) tmpData.size(),
RawValue_MPI_DOUBLE, 0,
2223 "MLSampling<P_V,P_M>::restartML()",
2224 "failed MPI.Bcast() for rest of information read from input file");
2225 if (m_env.fullRank() != 0) {
2226 vectorSpaceDim = tmpData[0];
2227 currExponent = tmpData[1];
2228 currEta = tmpData[2];
2229 quantity1 = tmpData[3];
2230 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2231 m_logEvidenceFactors[i] = tmpData[4+i];
2239 queso_require_msg(!((currExponent < 0.) || (currExponent > 1.)),
"read currExponent is not consistent");
2240 queso_require_equal_to_msg((quantity1 % m_env.numSubEnvironments()), 0,
"read size of chain should be a multiple of the number of subenvironments");
2241 unsigned int subSequenceSize = 0;
2242 subSequenceSize = ((double) quantity1) / ((double) m_env.numSubEnvironments());
2244 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2245 *m_env.subDisplayFile() <<
"Restart input file has the following information"
2246 <<
": subSequenceSize = " << subSequenceSize
2253 char levelSufix[256];
2256 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2257 *m_env.subDisplayFile() <<
"\n RESTARTING chain at level " << m_currLevel
2258 <<
"\n" << std::endl;
2260 currChain.
unifiedReadContents(m_options.m_restartInput_baseNameForFiles +
"Chain_l" + levelSufix,
2261 m_options.m_restartInput_fileType,
2263 m_env.fullComm().Barrier();
2265 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2266 *m_env.subDisplayFile() <<
"\n RESTARTING like at level " << m_currLevel
2267 <<
"\n" << std::endl;
2269 currLogLikelihoodValues.
unifiedReadContents(m_options.m_restartInput_baseNameForFiles +
"LogLike_l" + levelSufix,
2270 m_options.m_restartInput_fileType,
2272 m_env.fullComm().Barrier();
2274 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2275 *m_env.subDisplayFile() <<
"\n RESTARTING target at level " << m_currLevel
2276 <<
"\n" << std::endl;
2278 currLogTargetValues.
unifiedReadContents(m_options.m_restartInput_baseNameForFiles +
"LogTarget_l" + levelSufix,
2279 m_options.m_restartInput_fileType,
2281 m_env.fullComm().Barrier();
2283 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2284 *m_env.subDisplayFile() <<
"\n RESTARTING done at level " << m_currLevel
2285 <<
"\n" << std::endl;
2291 template <
class P_V,
class P_M>
2295 unsigned int& unifiedRequestedNumSamples,
2300 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2301 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateSequence()"
2303 <<
", currOptions.m_rawChainSize = " << currOptions.
m_rawChainSize
2308 struct timeval timevalLevel;
2309 iRC = gettimeofday(&timevalLevel, NULL);
2312 if (m_env.inter0Rank() >= 0) {
2314 m_env.inter0Comm().template Allreduce<unsigned int>(&tmpSize, &unifiedRequestedNumSamples, (int) 1,
RawValue_MPI_SUM,
2315 "MLSampling<P_V,P_M>::generateSequence()",
2316 "failed MPI.Allreduce() for requested num samples in level 0");
2323 currLogLikelihoodValues.
setName(currOptions.
m_prefix +
"rawLogLikelihood");
2330 P_V auxVec(m_vectorSpace.zeroVector());
2334 bool outOfSupport =
true;
2337 m_priorRv.realizer().realization(auxVec);
2338 if (m_numDisabledParameters > 0) {
2339 unsigned int disabledCounter = 0;
2340 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2341 if (m_parameterEnabledStatus[paramId] ==
false) {
2347 auxVec.mpiBcast(0, m_env.subComm());
2349 outOfSupport = !(m_targetDomain->contains(auxVec));
2350 }
while (outOfSupport);
2354 #if 1 // prudencio 2010-08-01
2355 currLogLikelihoodValues[i] = likelihoodSynchronizer.callFunction(&auxVec,NULL,NULL,NULL,NULL,NULL,NULL);
2357 currLogLikelihoodValues[i] = m_likelihoodFunction.lnValue(auxVec,NULL,NULL,NULL,NULL);
2359 currLogTargetValues[i] = m_priorRv.pdf().lnValue(auxVec,NULL,NULL,NULL,NULL) + currLogLikelihoodValues[i];
2363 if (m_env.inter0Rank() >= 0) {
2364 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2365 if (currOptions.m_rawChainComputeStats) {
2374 currChain.computeStatistics(*currOptions.m_rawChainStatisticalOptionsObj,
2391 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2392 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2395 <<
" chain positions"
2411 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2412 if (currOptions.m_filteredChainComputeStats) {
2428 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2429 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2431 <<
", total level time = " << levelRunTime <<
" seconds"
2438 template <
class P_V,
class P_M>
2442 unsigned int& unifiedRequestedNumSamples)
2445 struct timeval timevalStep;
2446 iRC = gettimeofday(&timevalStep, NULL);
2449 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2450 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2452 <<
", step " << m_currStep
2453 <<
": beginning step 1 of 11"
2460 m_env.inter0Comm().template Allreduce<unsigned int>(&tmpSize, &unifiedRequestedNumSamples, (int) 1,
RawValue_MPI_SUM,
2461 "MLSampling<P_V,P_M>::generateSequence()",
2462 "failed MPI.Allreduce() for requested num samples in step 1");
2464 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2465 *m_env.subDisplayFile() <<
"KEY In MLSampling<P_V,P_M>::generateSequence()"
2467 <<
", step " << m_currStep
2468 <<
", currOptions->m_rawChainSize = " << currOptions->
m_rawChainSize
2473 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2474 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2476 <<
", step " << m_currStep
2477 <<
", after " << stepRunTime <<
" seconds"
2484 template <
class P_V,
class P_M>
2494 unsigned int& indexOfFirstWeight,
2495 unsigned int& indexOfLastWeight)
2498 struct timeval timevalStep;
2499 iRC = gettimeofday(&timevalStep, NULL);
2502 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2503 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2505 <<
", step " << m_currStep
2506 <<
": beginning step 2 of 11"
2510 prevChain = currChain;
2514 prevLogLikelihoodValues = currLogLikelihoodValues;
2515 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2516 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
2518 <<
", step " << m_currStep
2519 <<
", prevLogLikelihoodValues[0] = " << prevLogLikelihoodValues[0]
2522 prevLogTargetValues = currLogTargetValues;
2524 currLogLikelihoodValues.
clear();
2525 currLogLikelihoodValues.
setName(currOptions->
m_prefix +
"rawLogLikelihood");
2527 currLogTargetValues.
clear();
2530 #if 0 // For debug only
2531 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2532 P_V prevPosition(m_vectorSpace.zeroVector());
2533 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2535 <<
", step " << m_currStep
2540 *m_env.subDisplayFile() <<
" prevChain[" << i
2541 <<
"] = " << prevPosition
2542 <<
", prevLogLikelihoodValues[" << i
2543 <<
"] = " << prevLogLikelihoodValues[i]
2544 <<
", prevLogTargetValues[" << i
2545 <<
"] = " << prevLogTargetValues[i]
2553 unsigned int quantity3 = prevLogLikelihoodValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2554 unsigned int quantity4 = currLogLikelihoodValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2555 unsigned int quantity5 = prevLogTargetValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2556 unsigned int quantity6 = currLogTargetValues.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2557 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2558 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2560 <<
", step " << m_currStep
2561 <<
": prevChain.unifiedSequenceSize() = " << quantity1
2562 <<
", currChain.unifiedSequenceSize() = " << quantity2
2563 <<
", prevLogLikelihoodValues.unifiedSequenceSize() = " << quantity3
2564 <<
", currLogLikelihoodValues.unifiedSequenceSize() = " << quantity4
2565 <<
", prevLogTargetValues.unifiedSequenceSize() = " << quantity5
2566 <<
", currLogTargetValues.unifiedSequenceSize() = " << quantity6
2575 indexOfFirstWeight = 0;
2576 indexOfLastWeight = indexOfFirstWeight + prevChain.
subSequenceSize()-1;
2579 int r = m_env.inter0Rank();
2581 m_env.inter0Comm().Barrier();
2582 unsigned int auxUint = 0;
2587 "MLSampling<P_V,P_M>::generateSequence()",
2588 "failed MPI.Recv()");
2590 indexOfFirstWeight = auxUint;
2591 indexOfLastWeight = indexOfFirstWeight + prevChain.
subSequenceSize()-1;
2593 if (r < (m_env.inter0Comm().NumProc()-1)) {
2594 auxUint = indexOfLastWeight + 1;
2597 "MLSampling<P_V,P_M>::generateSequence()",
2598 "failed MPI.Send()");
2601 m_env.inter0Comm().Barrier();
2605 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2606 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2608 <<
", step " << m_currStep
2609 <<
", after " << stepRunTime <<
" seconds"
2616 template <
class P_V,
class P_M>
2621 double prevExponent,
2622 double failedExponent,
2623 double& currExponent,
2627 struct timeval timevalStep;
2628 iRC = gettimeofday(&timevalStep, NULL);
2631 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2632 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2634 <<
", step " << m_currStep
2635 <<
", failedExponent = " << failedExponent
2636 <<
": beginning step 3 of 11"
2640 std::vector<double> exponents(2,0.);
2641 exponents[0] = prevExponent;
2644 double nowExponent = 1.;
2645 double nowEffectiveSizeRatio = 0.;
2647 unsigned int nowAttempt = 0;
2648 bool testResult =
false;
2652 double nowUnifiedEvidenceLnFactor = 0.;
2654 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2655 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2657 <<
", step " << m_currStep
2658 <<
", failedExponent = " << failedExponent
2659 <<
": entering loop for computing next exponent"
2660 <<
", with nowAttempt = " << nowAttempt
2664 if (failedExponent > 0.) {
2665 nowExponent = .5*(prevExponent+failedExponent);
2668 if (nowAttempt > 0) {
2669 if (nowEffectiveSizeRatio > meanEffectiveSizeRatio) {
2670 exponents[0] = nowExponent;
2673 exponents[1] = nowExponent;
2675 nowExponent = .5*(exponents[0] + exponents[1]);
2678 double auxExponent = nowExponent;
2679 if (prevExponent != 0.) {
2680 auxExponent /= prevExponent;
2683 double subWeightRatioSum = 0.;
2684 double unifiedWeightRatioSum = 0.;
2687 omegaLnDiffSequence[i] = prevLogLikelihoodValues[i]*auxExponent;
2690 #if 1 // prudenci-2012-07-06
2692 double unifiedOmegaLnMax = omegaLnDiffSequence.
unifiedMaxPlain(m_vectorSpace.numOfProcsForStorage() == 1);
2694 double unifiedOmegaLnMin = 0.;
2695 double unifiedOmegaLnMax = 0.;
2696 omegaLnDiffSequence.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
2698 omegaLnDiffSequence.subSequenceSize(),
2703 omegaLnDiffSequence[i] -= unifiedOmegaLnMax;
2704 weightSequence[i] = exp(omegaLnDiffSequence[i]);
2705 subWeightRatioSum += weightSequence[i];
2706 #if 0 // For debug only
2707 if ((m_currLevel == 1) && (nowAttempt == 6)) {
2708 if (m_env.subDisplayFile() && (m_env.displayVerbosity() >= 99)) {
2709 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2711 <<
", step " << m_currStep
2713 <<
", prevLogLikelihoodValues[i] = " << prevLogLikelihoodValues[i]
2714 <<
", omegaLnDiffSequence[i] = " << omegaLnDiffSequence[i]
2715 <<
", weightSequence[i] = " << weightSequence[i]
2722 m_env.inter0Comm().template Allreduce<double>(&subWeightRatioSum, &unifiedWeightRatioSum, (int) 1,
RawValue_MPI_SUM,
2723 "MLSampling<P_V,P_M>::generateSequence()",
2724 "failed MPI.Allreduce() for weight ratio sum");
2726 unsigned int auxQuantity = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2727 nowUnifiedEvidenceLnFactor = log(unifiedWeightRatioSum) + unifiedOmegaLnMax - log(auxQuantity);
2729 double effectiveSampleSize = 0.;
2731 weightSequence[i] /= unifiedWeightRatioSum;
2732 effectiveSampleSize += weightSequence[i]*weightSequence[i];
2743 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2744 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2746 <<
", step " << m_currStep
2747 <<
": nowAttempt = " << nowAttempt
2748 <<
", prevExponent = " << prevExponent
2749 <<
", exponents[0] = " << exponents[0]
2750 <<
", nowExponent = " << nowExponent
2751 <<
", exponents[1] = " << exponents[1]
2752 <<
", subWeightRatioSum = " << subWeightRatioSum
2753 <<
", unifiedWeightRatioSum = " << unifiedWeightRatioSum
2754 <<
", unifiedOmegaLnMax = " << unifiedOmegaLnMax
2755 <<
", weightSequence.unifiedSequenceSize() = " << auxQuantity
2756 <<
", nowUnifiedEvidenceLnFactor = " << nowUnifiedEvidenceLnFactor
2757 <<
", effectiveSampleSize = " << effectiveSampleSize
2761 #if 0 // For debug only
2762 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2763 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2765 <<
", step " << m_currStep
2770 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2771 *m_env.subDisplayFile() <<
" weightSequence[" << i
2772 <<
"] = " << weightSequence[i]
2778 double subQuantity = effectiveSampleSize;
2779 effectiveSampleSize = 0.;
2780 m_env.inter0Comm().template Allreduce<double>(&subQuantity, &effectiveSampleSize, (int) 1,
RawValue_MPI_SUM,
2781 "MLSampling<P_V,P_M>::generateSequence()",
2782 "failed MPI.Allreduce() for effective sample size");
2784 effectiveSampleSize = 1./effectiveSampleSize;
2785 nowEffectiveSizeRatio = effectiveSampleSize/((double) weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1));
2792 if (failedExponent > 0.) {
2797 bool aux2 = (nowExponent == 1. )
2799 (nowEffectiveSizeRatio > meanEffectiveSizeRatio);
2802 (nowEffectiveSizeRatio <= currOptions->m_maxEffectiveSizeRatio);
2803 testResult = aux2 || aux3;
2806 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2807 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2809 <<
", step " << m_currStep
2810 <<
": nowAttempt = " << nowAttempt
2811 <<
", prevExponent = " << prevExponent
2812 <<
", failedExponent = " << failedExponent
2813 <<
", exponents[0] = " << exponents[0]
2814 <<
", nowExponent = " << nowExponent
2815 <<
", exponents[1] = " << exponents[1]
2816 <<
", effectiveSampleSize = " << effectiveSampleSize
2819 <<
", nowEffectiveSizeRatio = " << nowEffectiveSizeRatio
2823 <<
", testResult = " << testResult
2832 "MLSampling<P_V,P_M>::generateSequence(), step 3, nowExponent") ==
false) {
2833 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2834 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2836 <<
", step " << m_currStep
2837 <<
": nowAttempt = " << nowAttempt
2838 <<
", MiscCheck for 'nowExponent' detected a problem"
2847 "MLSampling<P_V,P_M>::generateSequence(), step 3, testResult") ==
false) {
2848 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2849 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2851 <<
", step " << m_currStep
2852 <<
": nowAttempt = " << nowAttempt
2853 <<
", MiscCheck for 'testResult' detected a problem"
2857 }
while (testResult ==
false);
2858 currExponent = nowExponent;
2859 if (failedExponent > 0.) {
2860 m_logEvidenceFactors[m_logEvidenceFactors.size()-1] = nowUnifiedEvidenceLnFactor;
2863 m_logEvidenceFactors.push_back(nowUnifiedEvidenceLnFactor);
2866 unsigned int quantity1 = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2867 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2868 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2870 <<
", step " << m_currStep
2871 <<
": weightSequence.subSequenceSize() = " << weightSequence.
subSequenceSize()
2872 <<
", weightSequence.unifiedSequenceSize() = " << quantity1
2873 <<
", failedExponent = " << failedExponent
2874 <<
", currExponent = " << currExponent
2875 <<
", effective ratio = " << nowEffectiveSizeRatio
2876 <<
", log(evidence factor) = " << m_logEvidenceFactors[m_logEvidenceFactors.size()-1]
2877 <<
", evidence factor = " << exp(m_logEvidenceFactors[m_logEvidenceFactors.size()-1])
2895 "MLSampling<P_V,P_M>::generateSequence(), step 3, logEvidenceFactor") ==
false) {
2896 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2897 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2899 <<
", step " << m_currStep
2900 <<
", failedExponent = " << failedExponent
2901 <<
": nowAttempt = " << nowAttempt
2902 <<
", MiscCheck for 'logEvidenceFactor' detected a problem"
2908 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2909 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2911 <<
", step " << m_currStep
2912 <<
", failedExponent = " << failedExponent
2913 <<
", after " << stepRunTime <<
" seconds"
2920 template <
class P_V,
class P_M>
2925 P_M& unifiedCovMatrix)
2928 struct timeval timevalStep;
2929 iRC = gettimeofday(&timevalStep, NULL);
2932 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2933 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2935 <<
", step " << m_currStep
2936 <<
": beginning step 4 of 11"
2940 P_V auxVec(m_vectorSpace.zeroVector());
2941 P_V subWeightedMeanVec(m_vectorSpace.zeroVector());
2944 subWeightedMeanVec += weightSequence[i]*auxVec;
2948 P_V unifiedWeightedMeanVec(m_vectorSpace.zeroVector());
2949 if (m_env.inter0Rank() >= 0) {
2950 subWeightedMeanVec.mpiAllReduce(
RawValue_MPI_SUM,m_env.inter0Comm(),unifiedWeightedMeanVec);
2953 unifiedWeightedMeanVec = subWeightedMeanVec;
2956 P_V diffVec(m_vectorSpace.zeroVector());
2957 P_M subCovMatrix(m_vectorSpace.zeroVector());
2960 diffVec = auxVec - unifiedWeightedMeanVec;
2961 subCovMatrix += weightSequence[i]*
matrixProduct(diffVec,diffVec);
2964 for (
unsigned int i = 0; i < unifiedCovMatrix.numRowsLocal(); ++i) {
2965 for (
unsigned int j = 0; j < unifiedCovMatrix.numCols(); ++j) {
2966 double localValue = subCovMatrix(i,j);
2967 double sumValue = 0.;
2968 if (m_env.inter0Rank() >= 0) {
2969 m_env.inter0Comm().template Allreduce<double>(&localValue, &sumValue, (int) 1,
RawValue_MPI_SUM,
2970 "MLSampling<P_V,P_M>::generateSequence()",
2971 "failed MPI.Allreduce() for cov matrix");
2974 sumValue = localValue;
2976 unifiedCovMatrix(i,j) = sumValue;
2980 if (m_numDisabledParameters > 0) {
2981 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2982 if (m_parameterEnabledStatus[paramId] ==
false) {
2983 for (
unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2984 unifiedCovMatrix(i,paramId) = 0.;
2986 for (
unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2987 unifiedCovMatrix(paramId,j) = 0.;
2989 unifiedCovMatrix(paramId,paramId) = 1.;
2994 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2995 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
2997 <<
", step " << m_currStep
2998 <<
": unifiedCovMatrix = " << unifiedCovMatrix
3003 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3004 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3006 <<
", step " << m_currStep
3007 <<
", after " << stepRunTime <<
" seconds"
3014 template <
class P_V,
class P_M>
3017 unsigned int unifiedRequestedNumSamples,
3019 std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
3020 std::vector<double>& unifiedWeightStdVectorAtProc0Only)
3023 struct timeval timevalStep;
3024 iRC = gettimeofday(&timevalStep, NULL);
3027 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3028 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3030 <<
", step " << m_currStep
3031 <<
": beginning step 5 of 11"
3035 #if 0 // For debug only
3036 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3037 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3039 <<
", step " << m_currStep
3040 <<
", before weightSequence.getUnifiedContentsAtProc0Only()"
3045 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3046 *m_env.subDisplayFile() <<
", weightSequence[" << i
3047 <<
"] = " << weightSequence[i]
3054 unifiedWeightStdVectorAtProc0Only);
3056 #if 0 // For debug only
3057 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3058 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3060 <<
", step " << m_currStep
3061 <<
", after weightSequence.getUnifiedContentsAtProc0Only()"
3065 for (
unsigned int i = 0; i < unifiedWeightStdVectorAtProc0Only.size(); ++i) {
3066 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3067 *m_env.subDisplayFile() <<
" unifiedWeightStdVectorAtProc0Only[" << i
3068 <<
"] = " << unifiedWeightStdVectorAtProc0Only[i]
3073 sampleIndexes_proc0(unifiedRequestedNumSamples,
3074 unifiedWeightStdVectorAtProc0Only,
3075 unifiedIndexCountersAtProc0Only);
3077 unsigned int auxUnifiedSize = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3078 if (m_env.inter0Rank() == 0) {
3079 queso_require_equal_to_msg(unifiedIndexCountersAtProc0Only.size(), auxUnifiedSize,
"wrong output from sampleIndexesAtProc0() in step 5");
3083 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3084 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3086 <<
", step " << m_currStep
3087 <<
", after " << stepRunTime <<
" seconds"
3094 template <
class P_V,
class P_M>
3098 unsigned int indexOfFirstWeight,
3099 unsigned int indexOfLastWeight,
3100 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
3101 bool& useBalancedChains,
3102 std::vector<ExchangeInfoStruct>& exchangeStdVec)
3105 struct timeval timevalStep;
3106 iRC = gettimeofday(&timevalStep, NULL);
3109 useBalancedChains = decideOnBalancedChains_all(currOptions,
3112 unifiedIndexCountersAtProc0Only,
3116 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3117 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3119 <<
", step " << m_currStep
3120 <<
", after " << stepRunTime <<
" seconds"
3127 template <
class P_V,
class P_M>
3130 bool useBalancedChains,
3131 unsigned int indexOfFirstWeight,
3132 unsigned int indexOfLastWeight,
3133 const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only,
3137 double prevExponent,
3138 double currExponent,
3141 std::vector<ExchangeInfoStruct>& exchangeStdVec,
3145 struct timeval timevalStep;
3146 iRC = gettimeofday(&timevalStep, NULL);
3149 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3150 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3152 <<
", step " << m_currStep
3153 <<
": beginning step 7 of 11"
3157 if (useBalancedChains) {
3158 prepareBalLinkedChains_inter0(currOptions,
3162 prevLogLikelihoodValues,
3163 prevLogTargetValues,
3165 balancedLinkControl);
3168 prepareUnbLinkedChains_inter0(indexOfFirstWeight,
3170 unifiedIndexCountersAtProc0Only,
3171 unbalancedLinkControl);
3174 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3175 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3177 <<
", step " << m_currStep
3178 <<
": balancedLinkControl.balLinkedChains.size() = " << balancedLinkControl.
balLinkedChains.size()
3179 <<
", unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.
unbLinkedChains.size()
3184 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3185 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3187 <<
", step " << m_currStep
3188 <<
", after " << stepRunTime <<
" seconds"
3195 template <
class P_V,
class P_M>
3202 struct timeval timevalStep;
3203 iRC = gettimeofday(&timevalStep, NULL);
3206 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3207 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3209 <<
", step " << m_currStep
3210 <<
": beginning step 8 of 11"
3217 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3218 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3220 <<
", step " << m_currStep
3221 <<
", after " << stepRunTime <<
" seconds"
3228 template <
class P_V,
class P_M>
3232 double prevExponent,
3233 double currExponent,
3236 unsigned int indexOfFirstWeight,
3237 unsigned int indexOfLastWeight,
3238 const std::vector<double>& unifiedWeightStdVectorAtProc0Only,
3243 P_M& unifiedCovMatrix,
3247 struct timeval timevalStep;
3248 iRC = gettimeofday(&timevalStep, NULL);
3252 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3253 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3255 <<
", step " << m_currStep
3256 <<
": skipping step 9 of 11"
3261 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3262 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3264 <<
", step " << m_currStep
3265 <<
": beginning step 9 of 11"
3269 double beforeEta = prevEta;
3270 double beforeRejectionRate = 0.;
3271 bool beforeRejectionRateIsBelowRange =
true;
3273 double nowEta = prevEta;
3274 double nowRejectionRate = 0.;
3275 bool nowRejectionRateIsBelowRange =
true;
3277 std::vector<double> etas(2,0.);
3278 etas[0] = beforeEta;
3281 std::vector<double> rejs(2,0.);
3285 unsigned int nowAttempt = 0;
3286 bool testResult =
false;
3288 bool useMiddlePointLogicForEta =
false;
3289 P_M nowCovMatrix(unifiedCovMatrix);
3290 #if 0 // KAUST, to check
3291 std::vector<double> unifiedWeightStdVectorAtProc0Only(0);
3293 unifiedWeightStdVectorAtProc0Only);
3296 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3297 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3299 <<
", step " << m_currStep
3300 <<
": entering loop for assessing rejection rate"
3301 <<
", with nowAttempt = " << nowAttempt
3302 <<
", nowRejectionRate = " << nowRejectionRate
3305 nowCovMatrix = unifiedCovMatrix;
3307 if (nowRejectionRate < currOptions->m_minRejectionRate) {
3308 nowRejectionRateIsBelowRange =
true;
3311 nowRejectionRateIsBelowRange =
false;
3314 queso_error_msg(
"nowRejectionRate should be out of the requested range at this point of the logic");
3317 if (m_env.inter0Rank() >= 0) {
3318 if (nowAttempt > 0) {
3319 if (useMiddlePointLogicForEta ==
false) {
3320 if (nowAttempt == 1) {
3323 else if ((beforeRejectionRateIsBelowRange ==
true) &&
3324 (nowRejectionRateIsBelowRange ==
true)) {
3327 else if ((beforeRejectionRateIsBelowRange ==
false) &&
3328 (nowRejectionRateIsBelowRange ==
false)) {
3331 else if ((beforeRejectionRateIsBelowRange ==
true ) &&
3332 (nowRejectionRateIsBelowRange ==
false)) {
3333 useMiddlePointLogicForEta =
true;
3336 etas[0] = std::min(beforeEta,nowEta);
3337 etas[1] = std::max(beforeEta,nowEta);
3339 if (etas[0] == beforeEta) {
3340 rejs[0] = beforeRejectionRate;
3341 rejs[1] = nowRejectionRate;
3344 rejs[0] = nowRejectionRate;
3345 rejs[1] = beforeRejectionRate;
3348 else if ((beforeRejectionRateIsBelowRange ==
false) &&
3349 (nowRejectionRateIsBelowRange ==
true )) {
3350 useMiddlePointLogicForEta =
true;
3353 etas[0] = std::min(beforeEta,nowEta);
3354 etas[1] = std::max(beforeEta,nowEta);
3362 beforeRejectionRate = nowRejectionRate;
3363 beforeRejectionRateIsBelowRange = nowRejectionRateIsBelowRange;
3364 if (useMiddlePointLogicForEta ==
false) {
3365 if (beforeRejectionRateIsBelowRange) nowEta *= 4.;
3367 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3368 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3370 <<
", step " << m_currStep
3371 <<
": in loop for assessing rejection rate"
3372 <<
", with nowAttempt = " << nowAttempt
3373 <<
", useMiddlePointLogicForEta = false"
3374 <<
", nowEta just updated to value (to be tested) " << nowEta
3379 if (nowRejectionRate > meanRejectionRate) {
3380 if (rejs[0] > meanRejectionRate) {
3390 if (rejs[0] < meanRejectionRate) {
3399 nowEta = .5*(etas[0] + etas[1]);
3400 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3401 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3403 <<
", step " << m_currStep
3404 <<
": in loop for assessing rejection rate"
3405 <<
", with nowAttempt = " << nowAttempt
3406 <<
", useMiddlePointLogicForEta = true"
3407 <<
", nowEta just updated to value (to be tested) " << nowEta
3408 <<
", etas[0] = " << etas[0]
3409 <<
", etas[1] = " << etas[1]
3416 nowCovMatrix *= nowEta;
3420 unsigned int originalSubNumSamples = 1 + (
unsigned int) (doubSubNumSamples);
3421 double auxDouble = (double) originalSubNumSamples;
3422 if ((auxDouble - doubSubNumSamples) < 1.e-8) {
3423 originalSubNumSamples += 1;
3426 if (m_env.inter0Rank() >= 0) {
3427 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3428 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3430 <<
", step " << m_currStep
3431 <<
": in loop for assessing rejection rate"
3432 <<
", about to sample " << originalSubNumSamples <<
" indexes"
3433 <<
", meanRejectionRate = " << meanRejectionRate
3439 std::vector<unsigned int> nowUnifiedIndexCountersAtProc0Only(0);
3440 if (m_env.inter0Rank() >= 0) {
3441 unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3442 sampleIndexes_proc0(tmpUnifiedNumSamples,
3443 unifiedWeightStdVectorAtProc0Only,
3444 nowUnifiedIndexCountersAtProc0Only);
3446 unsigned int auxUnifiedSize = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3447 if (m_env.inter0Rank() == 0) {
3448 queso_require_equal_to_msg(nowUnifiedIndexCountersAtProc0Only.size(), auxUnifiedSize,
"wrong output from sampleIndexesAtProc0() in step 9");
3451 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3452 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3454 <<
", step " << m_currStep
3455 <<
": in loop for assessing rejection rate"
3456 <<
", about to distribute sampled assessment indexes"
3461 std::vector<ExchangeInfoStruct> exchangeStdVec(0);
3466 bool useBalancedChains = decideOnBalancedChains_all(currOptions,
3469 nowUnifiedIndexCountersAtProc0Only,
3472 if (m_env.inter0Rank() >= 0) {
3473 if (useBalancedChains) {
3474 prepareBalLinkedChains_inter0(currOptions,
3478 prevLogLikelihoodValues,
3479 prevLogTargetValues,
3484 prepareUnbLinkedChains_inter0(indexOfFirstWeight,
3486 nowUnifiedIndexCountersAtProc0Only,
3491 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3492 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3494 <<
", step " << m_currStep
3495 <<
": in loop for assessing rejection rate"
3496 <<
", about to generate assessment chain"
3502 m_options.m_prefix+
"now_chain");
3503 double nowRunTime = 0.;
3504 unsigned int nowRejections = 0;
3509 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3510 bool savedRawChainComputeStats = currOptions->m_rawChainComputeStats;
3517 if (m_env.displayVerbosity() >= 999999) {
3521 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3522 currOptions->m_rawChainComputeStats =
false;
3529 if (useBalancedChains) {
3530 generateBalLinkedChains_all(*currOptions,
3541 generateUnbLinkedChains_all(*currOptions,
3549 prevLogLikelihoodValues,
3550 prevLogTargetValues,
3561 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3562 currOptions->m_rawChainComputeStats = savedRawChainComputeStats;
3568 for (
unsigned int i = 0; i < nowBalLinkControl.
balLinkedChains.size(); ++i) {
3575 if (m_env.inter0Rank() >= 0) {
3577 unsigned int nowUnifiedRejections = 0;
3578 m_env.inter0Comm().template Allreduce<unsigned int>(&nowRejections, &nowUnifiedRejections, (int) 1,
RawValue_MPI_SUM,
3579 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3580 "failed MPI.Allreduce() for now rejections");
3582 #if 0 // Round Rock 2009 12 29
3583 unsigned int tmpUnifiedNumSamples = 0;
3585 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3586 "failed MPI.Allreduce() for num samples in step 9");
3589 unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3590 nowRejectionRate = ((double) nowUnifiedRejections) / ((double) tmpUnifiedNumSamples);
3595 (nowRejectionRate <= currOptions->m_maxRejectionRate);
3602 "MLSampling<P_V,P_M>::generateSequence_Step09_all(), step 9, testResult") ==
false) {
3603 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3604 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3606 <<
", step " << m_currStep
3607 <<
": nowAttempt = " << nowAttempt
3608 <<
", MiscCheck for 'testResult' detected a problem"
3615 unsigned int tmpUint = (
unsigned int) testResult;
3617 "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3618 "failed MPI.Bcast() for testResult");
3619 testResult = (bool) tmpUint;
3621 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3622 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3624 <<
", step " << m_currStep
3625 <<
": in loop for assessing rejection rate"
3626 <<
", nowAttempt = " << nowAttempt
3627 <<
", beforeEta = " << beforeEta
3628 <<
", etas[0] = " << etas[0]
3629 <<
", nowEta = " << nowEta
3630 <<
", etas[1] = " << etas[1]
3632 <<
", nowRejectionRate = " << nowRejectionRate
3638 if (m_env.inter0Rank() >= 0) {
3643 "MLSampling<P_V,P_M>::generateSequence_Step09_all(), step 9, nowEta") ==
false) {
3644 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3645 *m_env.subDisplayFile() <<
"WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3647 <<
", step " << m_currStep
3648 <<
": nowAttempt = " << nowAttempt
3649 <<
", MiscCheck for 'nowEta' detected a problem"
3654 }
while (testResult ==
false);
3656 if (currEta != 1.) {
3657 unifiedCovMatrix *= currEta;
3658 if (m_numDisabledParameters > 0) {
3659 for (
unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
3660 if (m_parameterEnabledStatus[paramId] ==
false) {
3661 for (
unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
3662 unifiedCovMatrix(i,paramId) = 0.;
3664 for (
unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
3665 unifiedCovMatrix(paramId,j) = 0.;
3667 unifiedCovMatrix(paramId,paramId) = 1.;
3673 unsigned int quantity1 = weightSequence.
unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3674 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3675 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3677 <<
", step " << m_currStep
3678 <<
": weightSequence.subSequenceSize() = " << weightSequence.
subSequenceSize()
3679 <<
", weightSequence.unifiedSequenceSize() = " << quantity1
3680 <<
", currEta = " << currEta
3681 <<
", assessed rejection rate = " << nowRejectionRate
3687 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3688 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3690 <<
", step " << m_currStep
3691 <<
", after " << stepRunTime <<
" seconds"
3698 template <
class P_V,
class P_M>
3702 const P_M& unifiedCovMatrix,
3704 bool useBalancedChains,
3706 unsigned int indexOfFirstWeight,
3708 double prevExponent,
3709 double currExponent,
3714 double& cumulativeRawChainRunTime,
3715 unsigned int& cumulativeRawChainRejections,
3720 struct timeval timevalStep;
3721 iRC = gettimeofday(&timevalStep, NULL);
3724 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3725 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3727 <<
", step " << m_currStep
3728 <<
": beginning step 10 of 11"
3729 <<
", currLogLikelihoodValues = " << currLogLikelihoodValues
3736 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3737 bool savedRawChainComputeStats = currOptions.m_rawChainComputeStats;
3742 if (m_env.displayVerbosity() >= 999999) {
3746 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3747 currOptions.m_rawChainComputeStats =
false;
3752 if (useBalancedChains) {
3753 generateBalLinkedChains_all(currOptions,
3756 balancedLinkControl,
3758 cumulativeRawChainRunTime,
3759 cumulativeRawChainRejections,
3760 currLogLikelihoodValues,
3761 currLogTargetValues);
3764 generateUnbLinkedChains_all(currOptions,
3767 unbalancedLinkControl,
3772 prevLogLikelihoodValues,
3773 prevLogTargetValues,
3775 cumulativeRawChainRunTime,
3776 cumulativeRawChainRejections,
3777 currLogLikelihoodValues,
3778 currLogTargetValues);
3781 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3782 double tmpValue = INFINITY;
3783 if (currLogLikelihoodValues) tmpValue = (*currLogLikelihoodValues)[0];
3784 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
3786 <<
", step " << m_currStep
3787 <<
", after chain generatrion"
3788 <<
", currLogLikelihoodValues[0] = " << tmpValue
3795 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3796 currOptions.m_rawChainComputeStats = savedRawChainComputeStats;
3801 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3802 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3804 <<
", step " << m_currStep
3805 <<
", after " << stepRunTime <<
" seconds"
3812 template <
class P_V,
class P_M>
3816 unsigned int unifiedRequestedNumSamples,
3817 unsigned int cumulativeRawChainRejections,
3821 unsigned int& unifiedNumberOfRejections)
3824 struct timeval timevalStep;
3825 iRC = gettimeofday(&timevalStep, NULL);
3828 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3829 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
3831 <<
", step " << m_currStep
3832 <<
": beginning step 11 of 11"
3838 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3839 if (currOptions->m_rawChainComputeStats) {
3847 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
3848 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
3850 <<
", step " << m_currStep
3851 <<
", calling computeStatistics for raw chain"
3852 <<
". Ofstream pointer value = " << filePtrSet.
ofsVar
3853 <<
", statistical options are"
3854 <<
"\n" << *currOptions->m_rawChainStatisticalOptionsObj
3858 currChain.computeStatistics(*currOptions->m_rawChainStatisticalOptionsObj,
3869 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3870 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
3872 <<
", step " << m_currStep
3873 <<
", before calling currLogLikelihoodValues.unifiedWriteContents()"
3874 <<
", currLogLikelihoodValues[0] = " << currLogLikelihoodValues[0]
3893 if (filterSpacing == 0) {
3900 currChain.
filter(filterInitialPos,
3904 currLogLikelihoodValues.
filter(filterInitialPos,
3906 currLogLikelihoodValues.
setName(currOptions->
m_prefix +
"filtLogLikelihood");
3908 currLogTargetValues.
filter(filterInitialPos,
3912 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3913 if (currOptions->m_filteredChainComputeStats) {
3914 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
3915 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence_Step()"
3917 <<
", step " << m_currStep
3918 <<
", calling computeStatistics for filtered chain"
3919 <<
". Ofstream pointer value = " << filePtrSet.
ofsVar
3920 <<
", statistical options are"
3921 <<
"\n" << *currOptions->m_rawChainStatisticalOptionsObj
3926 currChain.computeStatistics(*currOptions->m_filteredChainStatisticalOptionsObj,
3951 unsigned int unifiedGeneratedNumSamples = 0;
3952 m_env.inter0Comm().template Allreduce<unsigned int>(&tmpSize, &unifiedGeneratedNumSamples, (int) 1,
RawValue_MPI_SUM,
3953 "MLSampling<P_V,P_M>::generateSequence()",
3954 "failed MPI.Allreduce() for generated num samples in step 11");
3958 queso_require_equal_to_msg(unifiedGeneratedNumSamples, unifiedRequestedNumSamples,
"currChain (linked one) has been generated with invalid size");
3962 m_env.inter0Comm().template Allreduce<unsigned int>(&cumulativeRawChainRejections, &unifiedNumberOfRejections, (int) 1,
RawValue_MPI_SUM,
3963 "MLSampling<P_V,P_M>::generateSequence()",
3964 "failed MPI.Allreduce() for number of rejections");
3967 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3968 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3970 <<
", step " << m_currStep
3971 <<
", after " << stepRunTime <<
" seconds"
3979 template<
class P_V,
class P_M>
3985 m_env (priorRv.env()),
3986 m_priorRv (priorRv),
3987 m_likelihoodFunction(likelihoodFunction),
3988 m_vectorSpace (m_priorRv.imageSet().vectorSpace()),
3990 m_numDisabledParameters (0),
3991 m_parameterEnabledStatus(m_vectorSpace.dimLocal(),true),
3992 m_options (m_env,prefix),
3995 m_debugExponent (0.),
3996 m_logEvidenceFactors(0),
3998 m_meanLogLikelihood (0.),
4012 template<
class P_V,
class P_M>
4015 m_numDisabledParameters = 0;
4016 m_parameterEnabledStatus.clear();
4017 if (m_targetDomain)
delete m_targetDomain;
4023 template <
class P_V,
class P_M>
4030 struct timeval timevalRoutineBegin;
4032 iRC = gettimeofday(&timevalRoutineBegin, NULL);
4035 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4036 *m_env.subDisplayFile() <<
"Entering MLSampling<P_V,P_M>::generateSequence()"
4037 <<
", at " << ctime(&timevalRoutineBegin.tv_sec)
4038 <<
", after " << timevalRoutineBegin.tv_sec - m_env.timevalBegin().tv_sec
4039 <<
" seconds from queso environment instatiation..."
4046 double currExponent = 0.;
4047 double currEta = 1.;
4048 unsigned int currUnifiedRequestedNumSamples = 0;
4051 m_options.m_prefix+
"curr_chain");
4059 bool stopAtEndOfLevel =
false;
4060 char levelPrefix[256];
4071 if (m_options.m_restartInput_baseNameForFiles !=
".") {
4072 restartML(currExponent,
4075 currLogLikelihoodValues,
4076 currLogTargetValues);
4078 if (currExponent == 1.) {
4079 if (lastLevelOptions.m_parameterDisabledSet.size() > 0) {
4080 for (std::set<unsigned int>::iterator setIt = lastLevelOptions.m_parameterDisabledSet.begin(); setIt != lastLevelOptions.m_parameterDisabledSet.end(); ++setIt) {
4081 unsigned int paramId = *setIt;
4082 if (paramId < m_vectorSpace.dimLocal()) {
4083 m_numDisabledParameters++;
4084 m_parameterEnabledStatus[paramId] =
false;
4089 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4090 if (lastLevelOptions.m_rawChainComputeStats) {
4092 m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4094 lastLevelOptions.m_dataOutputAllowedSet,
4098 currChain.computeStatistics(*lastLevelOptions.m_rawChainStatisticalOptionsObj,
4108 currChain.
unifiedWriteContents (lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4109 currLogLikelihoodValues.
unifiedWriteContents(lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4110 currLogTargetValues.
unifiedWriteContents (lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4113 if (lastLevelOptions.m_filteredChainGenerate) {
4115 m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4117 lastLevelOptions.m_dataOutputAllowedSet,
4121 unsigned int filterInitialPos = (
unsigned int) (lastLevelOptions.m_filteredChainDiscardedPortion * (
double) currChain.
subSequenceSize());
4122 unsigned int filterSpacing = lastLevelOptions.m_filteredChainLag;
4123 if (filterSpacing == 0) {
4130 currChain.
filter(filterInitialPos,
4132 currChain.
setName(lastLevelOptions.m_prefix +
"filtChain");
4134 currLogLikelihoodValues.
filter(filterInitialPos,
4136 currLogLikelihoodValues.
setName(lastLevelOptions.m_prefix +
"filtLogLikelihood");
4138 currLogTargetValues.
filter(filterInitialPos,
4140 currLogTargetValues.
setName(lastLevelOptions.m_prefix +
"filtLogTarget");
4142 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4143 if (lastLevelOptions.m_filteredChainComputeStats) {
4144 currChain.computeStatistics(*lastLevelOptions.m_filteredChainStatisticalOptionsObj,
4153 currChain.
unifiedWriteContents (lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4154 currLogLikelihoodValues.
unifiedWriteContents(lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4155 currLogTargetValues.
unifiedWriteContents (lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4165 if (currOptions.m_parameterDisabledSet.size() > 0) {
4166 for (std::set<unsigned int>::iterator setIt = currOptions.m_parameterDisabledSet.begin(); setIt != currOptions.m_parameterDisabledSet.end(); ++setIt) {
4167 unsigned int paramId = *setIt;
4168 if (paramId < m_vectorSpace.dimLocal()) {
4169 m_numDisabledParameters++;
4170 m_parameterEnabledStatus[paramId] =
false;
4175 generateSequence_Level0_all(currOptions,
4176 currUnifiedRequestedNumSamples,
4178 currLogLikelihoodValues,
4179 currLogTargetValues);
4181 stopAtEndOfLevel = currOptions.m_stopAtEnd;
4182 bool performCheckpoint = stopAtEndOfLevel;
4183 if (m_options.m_restartOutput_levelPeriod > 0) {
4184 performCheckpoint = performCheckpoint || ( ((m_currLevel + 1) % m_options.m_restartOutput_levelPeriod) == 0 );
4186 if (performCheckpoint) {
4187 checkpointML(currExponent,
4190 currLogLikelihoodValues,
4191 currLogTargetValues);
4197 double minLogLike = -INFINITY;
4198 double maxLogLike = INFINITY;
4203 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4204 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4206 <<
", sub minLogLike = " << minLogLike
4207 <<
", sub maxLogLike = " << maxLogLike
4211 m_env.fullComm().Barrier();
4213 minLogLike = -INFINITY;
4214 maxLogLike = INFINITY;
4215 currLogLikelihoodValues.
unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4220 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4221 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4223 <<
", unified minLogLike = " << minLogLike
4224 <<
", unified maxLogLike = " << maxLogLike
4231 while ((currExponent < 1. ) &&
4232 (stopAtEndOfLevel ==
false)) {
4235 struct timeval timevalLevelBegin;
4237 iRC = gettimeofday(&timevalLevelBegin, NULL);
4239 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4240 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4242 <<
", at " << ctime(&timevalLevelBegin.tv_sec)
4243 <<
", after " << timevalLevelBegin.tv_sec - timevalRoutineBegin.tv_sec
4244 <<
" seconds from entering the routine"
4245 <<
", after " << timevalLevelBegin.tv_sec - m_env.timevalBegin().tv_sec
4246 <<
" seconds from queso environment instatiation"
4251 struct timeval timevalLevel;
4252 iRC = gettimeofday(&timevalLevel, NULL);
4253 double cumulativeRawChainRunTime = 0.;
4254 unsigned int cumulativeRawChainRejections = 0;
4256 bool tryExponentEta =
true;
4257 double failedExponent = 0.;
4258 double failedEta = 0.;
4264 double prevExponent = 0.;
4265 unsigned int indexOfFirstWeight = 0;
4266 unsigned int indexOfLastWeight = 0;
4267 P_M* unifiedCovMatrix = NULL;
4268 bool useBalancedChains =
false;
4274 unsigned int exponentEtaTriedAmount = 0;
4275 while (tryExponentEta) {
4276 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4277 *m_env.subDisplayFile() <<
"In IMLSampling<P_V,P_M>::generateSequence()"
4279 <<
", beginning 'do-while(tryExponentEta)"
4280 <<
": failedExponent = " << failedExponent
4281 <<
", failedEta = " << failedEta
4295 unsigned int paramId = *setIt;
4296 if (paramId < m_vectorSpace.dimLocal()) {
4297 m_numDisabledParameters++;
4298 m_parameterEnabledStatus[paramId] =
false;
4303 if (m_env.inter0Rank() >= 0) {
4304 generateSequence_Step01_inter0(currOptions,
4305 currUnifiedRequestedNumSamples);
4312 prevExponent = currExponent;
4313 double prevEta = currEta;
4314 unsigned int prevUnifiedRequestedNumSamples = currUnifiedRequestedNumSamples;
4317 m_options.m_prefix+
"prev_chain");
4319 indexOfFirstWeight = 0;
4320 indexOfLastWeight = 0;
4323 if (m_env.inter0Rank() >= 0) {
4324 generateSequence_Step02_inter0(currOptions,
4326 currLogLikelihoodValues,
4328 currLogTargetValues,
4330 prevLogLikelihoodValues,
4331 prevLogTargetValues,
4342 if (m_env.inter0Rank() >= 0) {
4343 generateSequence_Step03_inter0(currOptions,
4344 prevLogLikelihoodValues,
4353 "MLSampling<P_V,P_M>::generateSequence()",
4354 "failed MPI.Bcast() for currExponent");
4355 m_debugExponent = currExponent;
4357 if (currExponent == 1.) {
4358 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4359 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4361 <<
", step " << m_currStep
4362 <<
": copying 'last' level options to current options"
4366 currOptions = &lastLevelOptions;
4368 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4369 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4371 <<
", step " << m_currStep
4372 <<
": after copying 'last' level options to current options, the current options are"
4373 <<
"\n" << *currOptions
4377 if (m_env.inter0Rank() >= 0) {
4381 m_env.inter0Comm().template Allreduce<unsigned int>(&tmpSize, &currUnifiedRequestedNumSamples, (int) 1,
RawValue_MPI_SUM,
4382 "MLSampling<P_V,P_M>::generateSequence()",
4383 "failed MPI.Allreduce() for requested num samples in step 3");
4391 P_V oneVec(m_vectorSpace.zeroVector());
4394 unifiedCovMatrix = NULL;
4395 if (m_env.inter0Rank() >= 0) {
4396 unifiedCovMatrix = m_vectorSpace.newMatrix();
4399 unifiedCovMatrix =
new P_M(oneVec);
4402 if (m_env.inter0Rank() >= 0) {
4403 generateSequence_Step04_inter0(*prevChain,
4412 std::vector<unsigned int> unifiedIndexCountersAtProc0Only(0);
4413 std::vector<double> unifiedWeightStdVectorAtProc0Only(0);
4414 if (m_env.inter0Rank() >= 0) {
4415 generateSequence_Step05_inter0(currUnifiedRequestedNumSamples,
4417 unifiedIndexCountersAtProc0Only,
4418 unifiedWeightStdVectorAtProc0Only);
4425 useBalancedChains =
false;
4426 std::vector<ExchangeInfoStruct> exchangeStdVec(0);
4428 generateSequence_Step06_all(currOptions,
4431 unifiedIndexCountersAtProc0Only,
4442 if (m_env.inter0Rank() >= 0) {
4443 generateSequence_Step07_inter0(useBalancedChains,
4446 unifiedIndexCountersAtProc0Only,
4447 *unbalancedLinkControl,
4452 prevLogLikelihoodValues,
4453 prevLogTargetValues,
4455 *balancedLinkControl);
4466 m_likelihoodFunction,
4474 generateSequence_Step08_all(*currPdf,
4481 generateSequence_Step09_all(*prevChain,
4484 prevLogLikelihoodValues,
4485 prevLogTargetValues,
4488 unifiedWeightStdVectorAtProc0Only,
4496 tryExponentEta =
false;
4498 (currEta < currOptions->m_minAcceptableEta)) {
4499 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4500 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4502 <<
", preparing to retry ExponentEta"
4503 <<
": currExponent = " << currExponent
4504 <<
", currEta = " << currEta
4507 exponentEtaTriedAmount++;
4508 tryExponentEta =
true;
4509 failedExponent = currExponent;
4510 failedEta = currEta;
4518 delete balancedLinkControl;
4519 balancedLinkControl = NULL;
4520 delete unbalancedLinkControl;
4521 unbalancedLinkControl = NULL;
4523 delete unifiedCovMatrix;
4524 unifiedCovMatrix = NULL;
4526 currExponent = prevExponent;
4528 currUnifiedRequestedNumSamples = prevUnifiedRequestedNumSamples;
4531 currChain = (*prevChain);
4535 currLogLikelihoodValues = prevLogLikelihoodValues;
4536 currLogTargetValues = prevLogTargetValues;
4540 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4541 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4543 <<
", exited 'do-while(tryExponentEta)"
4544 <<
", failedExponent = " << failedExponent
4545 <<
", failedEta = " << failedEta
4554 generateSequence_Step10_all(*currOptions,
4558 *unbalancedLinkControl,
4563 prevLogLikelihoodValues,
4564 prevLogTargetValues,
4565 *balancedLinkControl,
4567 cumulativeRawChainRunTime,
4568 cumulativeRawChainRejections,
4569 &currLogLikelihoodValues,
4570 &currLogTargetValues);
4576 bool performCheckpoint = stopAtEndOfLevel;
4577 if (m_options.m_restartOutput_levelPeriod > 0) {
4578 performCheckpoint = performCheckpoint || ( ((m_currLevel + 1) % m_options.m_restartOutput_levelPeriod) == 0 );
4579 if (currExponent == 1.) {
4580 performCheckpoint =
true;
4583 if (performCheckpoint) {
4584 checkpointML(currExponent,
4587 currLogLikelihoodValues,
4588 currLogTargetValues);
4595 delete unifiedCovMatrix;
4597 for (
unsigned int i = 0; i < balancedLinkControl->
balLinkedChains.size(); ++i) {
4609 unsigned int unifiedNumberOfRejections = 0;
4610 if (m_env.inter0Rank() >= 0) {
4611 generateSequence_Step11_inter0(currOptions,
4612 currUnifiedRequestedNumSamples,
4613 cumulativeRawChainRejections,
4615 currLogLikelihoodValues,
4616 currLogTargetValues,
4617 unifiedNumberOfRejections);
4620 minLogLike = -INFINITY;
4621 maxLogLike = INFINITY;
4626 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4627 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4629 <<
", sub minLogLike = " << minLogLike
4630 <<
", sub maxLogLike = " << maxLogLike
4634 m_env.fullComm().Barrier();
4636 minLogLike = -INFINITY;
4637 maxLogLike = INFINITY;
4638 currLogLikelihoodValues.
unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4643 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4644 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4646 <<
", unified minLogLike = " << minLogLike
4647 <<
", unified maxLogLike = " << maxLogLike
4655 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4656 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4659 <<
" chain positions"
4660 <<
", cumulativeRawChainRunTime = " << cumulativeRawChainRunTime <<
" seconds"
4661 <<
", total level time = " << levelRunTime <<
" seconds"
4662 <<
", cumulativeRawChainRejections = " << cumulativeRawChainRejections
4663 <<
" (" << 100.*((double) cumulativeRawChainRejections)/((double) currOptions->
m_rawChainSize)
4664 <<
"% at this processor)"
4665 <<
" (" << 100.*((double) unifiedNumberOfRejections)/((double) currUnifiedRequestedNumSamples)
4666 <<
"% over all processors)"
4667 <<
", stopAtEndOfLevel = " << stopAtEndOfLevel
4671 if (m_env.inter0Rank() >= 0) {
4672 double minCumulativeRawChainRunTime = 0.;
4673 m_env.inter0Comm().template Allreduce<double>(&cumulativeRawChainRunTime, &minCumulativeRawChainRunTime, (int) 1,
RawValue_MPI_MIN,
4674 "MLSampling<P_V,P_M>::generateSequence()",
4675 "failed MPI.Allreduce() for min cumulative raw chain run time");
4677 double maxCumulativeRawChainRunTime = 0.;
4678 m_env.inter0Comm().template Allreduce<double>(&cumulativeRawChainRunTime, &maxCumulativeRawChainRunTime, (int) 1,
RawValue_MPI_MAX,
4679 "MLSampling<P_V,P_M>::generateSequence()",
4680 "failed MPI.Allreduce() for max cumulative raw chain run time");
4682 double avgCumulativeRawChainRunTime = 0.;
4683 m_env.inter0Comm().template Allreduce<double>(&cumulativeRawChainRunTime, &avgCumulativeRawChainRunTime, (int) 1,
RawValue_MPI_SUM,
4684 "MLSampling<P_V,P_M>::generateSequence()",
4685 "failed MPI.Allreduce() for sum cumulative raw chain run time");
4686 avgCumulativeRawChainRunTime /= ((double) m_env.inter0Comm().NumProc());
4688 double minLevelRunTime = 0.;
4689 m_env.inter0Comm().template Allreduce<double>(&levelRunTime, &minLevelRunTime, (int) 1,
RawValue_MPI_MIN,
4690 "MLSampling<P_V,P_M>::generateSequence()",
4691 "failed MPI.Allreduce() for min level run time");
4693 double maxLevelRunTime = 0.;
4694 m_env.inter0Comm().template Allreduce<double>(&levelRunTime, &maxLevelRunTime, (int) 1,
RawValue_MPI_MAX,
4695 "MLSampling<P_V,P_M>::generateSequence()",
4696 "failed MPI.Allreduce() for max level run time");
4698 double avgLevelRunTime = 0.;
4699 m_env.inter0Comm().template Allreduce<double>(&levelRunTime, &avgLevelRunTime, (int) 1,
RawValue_MPI_SUM,
4700 "MLSampling<P_V,P_M>::generateSequence()",
4701 "failed MPI.Allreduce() for sum level run time");
4702 avgLevelRunTime /= ((double) m_env.inter0Comm().NumProc());
4704 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4705 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4707 <<
": min cumul seconds = " << minCumulativeRawChainRunTime
4708 <<
", avg cumul seconds = " << avgCumulativeRawChainRunTime
4709 <<
", max cumul seconds = " << maxCumulativeRawChainRunTime
4710 <<
", min level seconds = " << minLevelRunTime
4711 <<
", avg level seconds = " << avgLevelRunTime
4712 <<
", max level seconds = " << maxLevelRunTime
4717 if (currExponent != 1.)
delete currOptions;
4719 struct timeval timevalLevelEnd;
4721 iRC = gettimeofday(&timevalLevelEnd, NULL);
4723 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4724 *m_env.subDisplayFile() <<
"Getting at the end of level " << m_currLevel+
LEVEL_REF_ID
4725 <<
", as part of a 'while' on levels"
4726 <<
", at " << ctime(&timevalLevelEnd.tv_sec)
4727 <<
", after " << timevalLevelEnd.tv_sec - timevalRoutineBegin.tv_sec
4728 <<
" seconds from entering the routine"
4729 <<
", after " << timevalLevelEnd.tv_sec - m_env.timevalBegin().tv_sec
4730 <<
" seconds from queso environment instatiation"
4744 if (m_env.inter0Rank() >= 0) {
4747 for (
unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
4748 m_logEvidence += m_logEvidenceFactors[i];
4751 #if 1 // prudenci-2012-07-06
4752 m_meanLogLikelihood = currLogLikelihoodValues.
unifiedMeanPlain(m_vectorSpace.numOfProcsForStorage() == 1);
4754 m_meanLogLikelihood = currLogLikelihoodValues.
unifiedMeanExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4759 m_eig = m_meanLogLikelihood - m_logEvidence;
4761 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4762 *m_env.subDisplayFile() <<
"In MLSampling<P_V,P_M>::generateSequence()"
4763 <<
", log(evidence) = " << m_logEvidence
4764 <<
", evidence = " << exp(m_logEvidence)
4765 <<
", meanLogLikelihood = " << m_meanLogLikelihood
4766 <<
", eig = " << m_eig
4772 "MLSampling<P_V,P_M>::generateSequence()",
4773 "failed MPI.Bcast() for m_logEvidence");
4776 "MLSampling<P_V,P_M>::generateSequence()",
4777 "failed MPI.Bcast() for m_meanLogLikelihood");
4780 "MLSampling<P_V,P_M>::generateSequence()",
4781 "failed MPI.Bcast() for m_eig");
4786 workingChain.
clear();
4788 P_V auxVec(m_vectorSpace.zeroVector());
4790 if (m_env.inter0Rank() >= 0) {
4796 if (workingLogLikelihoodValues) *workingLogLikelihoodValues = currLogLikelihoodValues;
4797 if (workingLogTargetValues ) *workingLogTargetValues = currLogTargetValues;
4799 struct timeval timevalRoutineEnd;
4801 iRC = gettimeofday(&timevalRoutineEnd, NULL);
4803 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4804 *m_env.subDisplayFile() <<
"Leaving MLSampling<P_V,P_M>::generateSequence()"
4805 <<
", at " << ctime(&timevalRoutineEnd.tv_sec)
4806 <<
", after " << timevalRoutineEnd.tv_sec - timevalRoutineBegin.tv_sec
4807 <<
" seconds from entering the routine"
4808 <<
", after " << timevalRoutineEnd.tv_sec - m_env.timevalBegin().tv_sec
4809 <<
" seconds from queso environment instatiation"
4816 template<
class P_V,
class P_M>
4817 std::ostream& operator<<(std::ostream& os, const MLSampling<P_V,P_M>& obj)
4824 template <
class P_V,
class P_M>
4827 return m_logEvidence;
4830 template <
class P_V,
class P_M>
4833 return m_meanLogLikelihood;
4836 template <
class P_V,
class P_M>
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
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)...
bool m_stopAtEnd
Stop at end of such level.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
void sampleIndexes_proc0(unsigned int unifiedRequestedNumSamples, const std::vector< double > &unifiedWeightStdVectorAtProc0Only, std::vector< unsigned int > &unifiedIndexCountersAtProc0Only)
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).
#define RawValue_MPI_UNSIGNED
int finalNodeOfInitialPosition
std::string m_filteredChainDataOutputFileType
Type of output file for filtered chain.
void append(const BaseVectorSequence< V, M > &src, unsigned int initialPos, unsigned int numPos)
Appends the vector src to this vector.
void getUnifiedContentsAtProc0Only(bool useOnlyInter0Comm, std::vector< T > &outputVec) const
Gets the unified contents of processor of rank equals to 0.
MLSampling(const char *prefix, const BaseVectorRV< P_V, P_M > &priorRv, const BaseScalarFunction< P_V, P_M > &likelihoodFunction)
Constructor.
#define queso_error_msg(msg)
double m_minEffectiveSizeRatio
Minimum allowed effective size ratio wrt previous level.
bool m_initialPositionUsePreviousLevelLikelihood
Use previous level likelihood for initial chain position instead of re-computing it from target pdf...
unsigned int m_filteredChainLag
Spacing for chain filtering.
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
unsigned int initialPositionIndexInPreviousChain
double meanLogLikelihood() const
Method to calculate the mean of the logarithm of the likelihood.
std::vector< double > m_initialValuesOfDisabledParameters
void clear()
Clears the sequence of scalars.
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.
unsigned int numberOfPositions
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
#define queso_require_less_msg(expr1, expr2, msg)
double m_loadBalanceTreshold
Perform load balancing if load unbalancing ratio > threshold.
std::string m_filteredChainDataOutputFileName
Name of output file for filtered chain.
std::set< unsigned int > m_dataOutputAllowedSet
subEnvs that will write to generic output file.
void justBalance_proc0(const MLSamplingLevelOptions *currOptions, std::vector< ExchangeInfoStruct > &exchangeStdVec)
double m_filteredChainDiscardedPortion
Initial discarded portion for chain filtering.
double m_covRejectionRate
c.o.v. for judging attempted rejection rate at current level.
void restartML(double &currExponent, double &currEta, SequenceOfVectors< P_V, P_M > &currChain, ScalarSequence< double > &currLogLikelihoodValues, ScalarSequence< double > &currLogTargetValues)
Restarts ML algorithm.
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).
std::string m_rawChainDataOutputFileType
Type of output file for raw chain.
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
#define ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
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).
A templated class for a finite distribution.
double logEvidence() const
Method to calculate the logarithm of the evidence.
std::string m_dataOutputFileName
Name of generic output file.
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)
double MiscGetEllapsedSeconds(struct timeval *timeval0)
unsigned int numberOfPositions
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
#define queso_require_less_equal_msg(expr1, expr2, msg)
#define queso_require_msg(asserted, msg)
#define queso_require_equal_to_msg(expr1, expr2, msg)
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
double eig() const
Calculates the expected information gain value, EIG.
#define queso_require_greater_equal_msg(expr1, expr2, msg)
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...
double m_maxEffectiveSizeRatio
Maximum allowed effective size ratio wrt previous level.
void clear()
Reset the values and the size of the sequence of vectors.
unsigned int numRejections
const BaseEnvironment & m_env
Queso enviroment.
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.
std::vector< UnbalancedLinkedChainControlStruct > unbLinkedChains
unsigned int m_drMaxNumExtraStages
'dr' maximum number of extra stages.
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
bool m_scaleCovMatrix
Whether or not scale proposal covariance matrix.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
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 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...
Struct for handling data input and output from files.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
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).
void append(const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
Appends the scalar sequence src to this sequence.
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
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...
bool decideOnBalancedChains_all(const MLSamplingLevelOptions *currOptions, unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, std::vector< ExchangeInfoStruct > &exchangeStdVec)
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)
unsigned int sample() const
Samples.
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
bool MiscCheckForSameValueInAllNodes(T &inputValue, double acceptableTreshold, const MpiComm &comm, const char *whereString)
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)...
std::set< unsigned int > m_parameterDisabledSet
void generateSequence_Step01_inter0(const MLSamplingLevelOptions *currOptions, unsigned int &unifiedRequestedNumSamples)
Reads options for the ML algorithm (Step 01 from ML algorithm).
A templated class that represents a Metropolis-Hastings generator of samples.
if the work is an executable linked with the with the complete machine readable work that uses the as object code and or source so that the user can modify the Library and then relink to produce a modified executable containing the modified rather than copying library functions into the if the user installs as long as the modified version is interface compatible with the version that the work was made with c Accompany the work with a written valid for at least three to give the same user the materials specified in for a charge no more than the cost of performing this distribution d If distribution of the work is made by offering access to copy from a designated offer equivalent access to copy the above specified materials from the same place e Verify that the user has already received a copy of these materials or that you have already sent this user a copy For an the required form of the work that uses the Library must include any data and utility programs needed for reproducing the executable from it as a special the materials to be distributed need not include anything that is normally and so on of the operating system on which the executable unless that component itself accompanies the executable It may happen that this requirement contradicts the license restrictions of other proprietary libraries that do not normally accompany the operating system Such a contradiction means you cannot use both them and the Library together in an executable that you distribute You may place library facilities that are a work based on the Library side by side in a single library together with other library facilities not covered by this and distribute such a combined provided that the separate distribution of the work based on the Library and of the other library facilities is otherwise and provided that you do these two uncombined with any other library facilities This must be distributed under the terms of the Sections above b Give prominent notice with the combined library of the fact that part of it is a work based on the and explaining where to find the accompanying uncombined form of the same work You may not link or distribute the Library except as expressly provided under this License Any attempt otherwise to link or distribute the Library is and will automatically terminate your rights under this License parties who have received or from you under this License will not have their licenses terminated so long as such parties remain in full compliance You are not required to accept this since you have not signed it nothing else grants you permission to modify or distribute the Library or its derivative works These actions are prohibited by law if you do not accept this License by modifying or distributing the you indicate your acceptance of this License to do and all its terms and conditions for distributing or modifying the Library or works based on it Each time you redistribute the the recipient automatically receives a license from the original licensor to link with or modify the Library subject to these terms and conditions You may not impose any further restrictions on the recipients exercise of the rights granted herein You are not responsible for enforcing compliance by third parties with this License as a consequence of a court judgment or allegation of patent infringement or for any other reason(not limited to patent issues)
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
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).
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
A struct that represents a Metropolis-Hastings sample.
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
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...
bool m_filteredChainGenerate
Whether or not to generate filtered chain.
void scanOptionsValues(const MLSamplingLevelOptions *defaultOptions)
unsigned int originalIndexOfInitialPosition
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
This class provides options for each level of the Multilevel sequence generator if no input file is a...
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 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)
double m_minRejectionRate
minimum allowed attempted rejection rate at current level
bool m_totallyMute
Whether or not to be totally mute (no printout message).
void prepareUnbLinkedChains_inter0(unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, UnbalancedLinkedChainsPerNodeStruct &unbalancedLinkControl)
unsigned int m_rawChainSize
Size of raw chain.
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
int originalNodeOfInitialPosition
double m_minAcceptableEta
minimum acceptable eta,
unsigned int m_loadBalanceAlgorithmId
Perform load balancing with chosen algorithm (0 = no balancing).
std::string m_rawChainDataOutputFileName
Name of output file for raw chain.
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
#define RawValue_MPI_CHAR
double m_maxRejectionRate
maximum allowed attempted rejection rate at current level.
std::vector< BalancedLinkedChainControlStruct< P_V > > balLinkedChains
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
#define RawValue_MPI_DOUBLE
unsigned int unifiedSequenceSize() const
Calculates the size of the unified 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.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
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)
unsigned int m_amAdaptInterval
'am' adaptation interval.
A templated class that represents a Multilevel generator of samples.