queso-0.55.0
MLSampling.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/MLSampling.h>
26 #include <queso/InstantiateIntersection.h>
27 #include <queso/GslVector.h>
28 #include <queso/GslMatrix.h>
29 
30 namespace QUESO {
31 
32 #ifdef QUESO_HAS_GLPK
33 
34 void BIP_routine(glp_tree *tree, void *info)
35 {
36  const BaseEnvironment& env = *(((BIP_routine_struct *) info)->env);
37  unsigned int currLevel = ((BIP_routine_struct *) info)->currLevel;
38 
39  int reason = glp_ios_reason(tree);
40 
41  if ((env.subDisplayFile()) && (env.displayVerbosity() >= 1)) {
42  *env.subDisplayFile() << "In BIP_routine()"
43  << ", level " << currLevel+LEVEL_REF_ID
44  << ": glp_ios_reason() = " << reason
45  << std::endl;
46  }
47  std::cout << "In BIP_routine: reason = " << reason << std::endl;
48 
49  switch (reason) {
50  case GLP_IROWGEN: // 0x01 /* request for row generation */
51  sleep(1);
52  break;
53 
54  case GLP_IBINGO: // 0x02 /* better integer solution found */
55  sleep(1);
56  break;
57 
58  case GLP_IHEUR: // 0x03 /* request for heuristic solution */
59  // Do nothing
60  break;
61 
62  case GLP_ICUTGEN: // 0x04 /* request for cut generation */
63  // Do nothing
64  break;
65 
66  case GLP_IBRANCH: // 0x05 /* request for branching */
67  // Do nothing
68  break;
69 
70  case GLP_ISELECT: // 0x06 /* request for subproblem selection */
71  // Do nothing
72  break;
73 
74  case GLP_IPREPRO: // 0x07 /* request for preprocessing */
75  // Do nothing
76  break;
77 
78  default:
79  queso_error_msg("invalid glp_ios_readon");
80  break;
81  }
82 
83  return;
84 }
85 
86 #endif // QUESO_HAS_GLPK
87 
88 template <class P_V,class P_M>
89 void
90 MLSampling<P_V,P_M>::sampleIndexes_proc0(
91  unsigned int unifiedRequestedNumSamples, // input
92  const std::vector<double>& unifiedWeightStdVectorAtProc0Only, // input
93  std::vector<unsigned int>& unifiedIndexCountersAtProc0Only) // output
94 {
95  if (m_env.inter0Rank() != 0) return;
96 
97  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
98  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::sampleIndexes_proc0()"
99  << ", level " << m_currLevel+LEVEL_REF_ID
100  << ", step " << m_currStep
101  << ": unifiedRequestedNumSamples = " << unifiedRequestedNumSamples
102  << ", unifiedWeightStdVectorAtProc0Only.size() = " << unifiedWeightStdVectorAtProc0Only.size()
103  << std::endl;
104  }
105 
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()"
109  << ", level " << m_currLevel+LEVEL_REF_ID
110  << ", step " << m_currStep
111  << ":"
112  << std::endl;
113  unsigned int numZeros = 0;
114  for (unsigned int i = 0; i < unifiedWeightStdVectorAtProc0Only.size(); ++i) {
115  *m_env.subDisplayFile() << " unifiedWeightStdVectorAtProc0Only[" << i
116  << "] = " << unifiedWeightStdVectorAtProc0Only[i]
117  << std::endl;
118  if (unifiedWeightStdVectorAtProc0Only[i] == 0.) numZeros++;
119  }
120  *m_env.subDisplayFile() << "Number of zeros in unifiedWeightStdVectorAtProc0Only = " << numZeros
121  << std::endl;
122  }
123 #endif
124 
125  if (m_env.inter0Rank() == 0) {
126  unsigned int resizeSize = unifiedWeightStdVectorAtProc0Only.size();
127  unifiedIndexCountersAtProc0Only.resize(resizeSize,0);
128 
129  // Generate 'unifiedRequestedNumSamples' samples from 'tmpFD'
130  FiniteDistribution tmpFd(m_env,
131  "",
132  unifiedWeightStdVectorAtProc0Only);
133  for (unsigned int i = 0; i < unifiedRequestedNumSamples; ++i) {
134  unsigned int index = tmpFd.sample();
135  unifiedIndexCountersAtProc0Only[index] += 1;
136  }
137  }
138 
139  return;
140 }
141 
142 template <class P_V,class P_M>
143 bool
144 MLSampling<P_V,P_M>::decideOnBalancedChains_all(
145  const MLSamplingLevelOptions* currOptions, // input
146  unsigned int indexOfFirstWeight, // input
147  unsigned int indexOfLastWeight, // input
148  const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // input
149  std::vector<ExchangeInfoStruct>& exchangeStdVec) // output
150 {
151  bool result = false;
152 
153  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
154  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
155  << ", level " << m_currLevel+LEVEL_REF_ID
156  << ", step " << m_currStep
157  << ": indexOfFirstWeight = " << indexOfFirstWeight
158  << ", indexOfLastWeight = " << indexOfLastWeight
159  << std::endl;
160  }
161 
162  if (true) {
163  unsigned int Np = 0;
164  if (m_env.inter0Rank() >= 0) { // Yes, '>= 0'
165  Np = (unsigned int) m_env.inter0Comm().NumProc();
166  }
167  std::vector<unsigned int> allFirstIndexes(Np,0); // '0' is already the correct value for recvcnts[0]
168  std::vector<unsigned int> allLastIndexes(Np,0); // '0' is NOT the correct value for recvcnts[0]
169 
170  if (m_env.inter0Rank() >= 0) { // Yes, '>= 0'
172  // Gather information at proc 0: number of chains and positions per node
174  unsigned int auxUInt = indexOfFirstWeight;
175  m_env.inter0Comm().template Gather<unsigned int>(&auxUInt, 1, &allFirstIndexes[0], (int) 1, 0, // LOAD BALANCE
176  "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
177  "failed MPI.Gather() for first indexes");
178 
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");
181  }
182 
183  auxUInt = indexOfLastWeight;
184  m_env.inter0Comm().template Gather<unsigned int>(&auxUInt, 1, &allLastIndexes[0], (int) 1, 0, // LOAD BALANCE
185  "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
186  "failed MPI.Gather() for last indexes");
187 
188  if (m_env.inter0Rank() == 0) { // Yes, '== 0'
189  //allLastIndexes[0] = indexOfLastWeight; // FIX ME: really necessary????
190  queso_require_equal_to_msg(allLastIndexes[0], indexOfLastWeight, "failed MPI.Gather() result for last indexes, at proc 0");
191  }
192  }
193 
195  // Proc 0 prepares information to decide if load balancing is needed
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()"
200  << ", level " << m_currLevel+LEVEL_REF_ID
201  << ", step " << m_currStep
202  << ": original distribution of unified indexes in 'inter0Comm' is as follows"
203  << std::endl;
204  for (unsigned int r = 0; r < Np; ++r) {
205  *m_env.subDisplayFile() << " allFirstIndexes[" << r << "] = " << allFirstIndexes[r]
206  << " allLastIndexes[" << r << "] = " << allLastIndexes[r]
207  << std::endl;
208  }
209  }
210  for (unsigned int r = 0; r < (Np-1); ++r) { // Yes, '-1'
211  queso_require_equal_to_msg(allFirstIndexes[r+1], (allLastIndexes[r]+1), "wrong indexes");
212  }
213 
214  for (unsigned int r = 0; r < (Np-1); ++r) { // Yes, '-1'
215  queso_require_equal_to_msg(allFirstIndexes[r+1], (allLastIndexes[r]+1), "wrong indexes");
216  }
217 
218  std::vector<unsigned int> origNumChainsPerNode (Np,0);
219  std::vector<unsigned int> origNumPositionsPerNode(Np,0);
220  int r = 0;
221  for (unsigned int i = 0; i < unifiedIndexCountersAtProc0Only.size(); ++i) {
222  if ((allFirstIndexes[r] <= i) && // FIX ME: not a robust logic
223  (i <= allLastIndexes[r] )) {
224  // Ok
225  }
226  else {
227  r++;
228  if ((r < (int) Np ) &&
229  (allFirstIndexes[r] <= i) &&
230  (i <= allLastIndexes[r] )) {
231  // Ok
232  }
233  else {
234  std::cerr << "In MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
235  << ": i = " << i
236  << ", r = " << r
237  << ", allFirstIndexes[r] = " << allFirstIndexes[r]
238  << ", allLastIndexes[r] = " << allLastIndexes[r]
239  << std::endl;
240  queso_error_msg("wrong indexes or 'r' got too large");
241  }
242  }
243  if (unifiedIndexCountersAtProc0Only[i] != 0) {
244  origNumChainsPerNode [r] += 1;
245  origNumPositionsPerNode[r] += unifiedIndexCountersAtProc0Only[i];
246 
247  ExchangeInfoStruct auxInfo;
248  auxInfo.originalNodeOfInitialPosition = r;
249  auxInfo.originalIndexOfInitialPosition = i - allFirstIndexes[r];
250  auxInfo.finalNodeOfInitialPosition = -1; // Yes, '-1' for now, important
251  auxInfo.numberOfPositions = unifiedIndexCountersAtProc0Only[i];
252  exchangeStdVec.push_back(auxInfo);
253  }
254  // FIX ME: swap trick to save memory
255  }
256 
257  // Check if number of procs is too large
258  unsigned int totalNumberOfChains = 0;
259  for (unsigned int r = 0; r < Np; ++r) {
260  totalNumberOfChains += origNumChainsPerNode[r];
261  }
262  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
263  *m_env.subDisplayFile() << " KEY"
264  << ", level " << m_currLevel+LEVEL_REF_ID
265  << ", step " << m_currStep
266  << ", Np = " << Np
267  << ", totalNumberOfChains = " << totalNumberOfChains
268  << std::endl;
269  }
270 
271  // Check if ratio max/min justifies optimization
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"
277  << ", level " << m_currLevel+LEVEL_REF_ID
278  << ", step " << m_currStep
279  << ", origNumChainsPerNode[" << nodeId << "] = " << origNumChainsPerNode[nodeId]
280  << ", origNumPositionsPerNode[" << nodeId << "] = " << origNumPositionsPerNode[nodeId]
281  << std::endl;
282  }
283  }
284  double origRatioOfPosPerNode = ((double) origMaxPosPerNode ) / ((double) origMinPosPerNode);
285  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
286  *m_env.subDisplayFile() << " KEY"
287  << ", level " << m_currLevel+LEVEL_REF_ID
288  << ", step " << m_currStep
289  << ", origRatioOfPosPerNode = " << origRatioOfPosPerNode
290  << ", option loadBalanceTreshold = " << currOptions->m_loadBalanceTreshold
291  << std::endl;
292  }
293 
294  // At this point, only proc 0 is running...
295  // Set boolean 'result' for good
296  if ((currOptions->m_loadBalanceAlgorithmId > 0 ) &&
297  (m_env.numSubEnvironments() > 1 ) && // Cannot use 'm_env.inter0Comm().NumProc()': not all nodes at this point of the code belong to 'inter0Comm'
298  (Np < totalNumberOfChains ) &&
299  (origRatioOfPosPerNode > currOptions->m_loadBalanceTreshold)) {
300  result = true;
301  }
302  } // if (m_env.inter0Rank() == 0)
303  } // if (true)
304 
305  m_env.fullComm().Barrier();
306  unsigned int tmpValue = result;
307  m_env.fullComm().Bcast((void *) &tmpValue, (int) 1, RawValue_MPI_UNSIGNED, 0, // LOAD BALANCE
308  "MLSampling<P_V,P_M>::decideOnBalancedChains_all()",
309  "failed MPI.Bcast() for 'result'");
310  if (m_env.inter0Rank() != 0) result = tmpValue;
311 
312  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
313  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::decideOnBalancedChains_all()"
314  << ", level " << m_currLevel+LEVEL_REF_ID
315  << ", step " << m_currStep
316  << ": result = " << result
317  << std::endl;
318  }
319 
320  return result;
321 }
322 
323 template <class P_V,class P_M>
324 void
325 MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0( // EXTRA FOR LOAD BALANCE
326  const MLSamplingLevelOptions* currOptions, // input
327  const SequenceOfVectors<P_V,P_M>& prevChain, // input
328  double prevExponent, // input
329  double currExponent, // input
330  const ScalarSequence<double>& prevLogLikelihoodValues, // input
331  const ScalarSequence<double>& prevLogTargetValues, // input
332  std::vector<ExchangeInfoStruct>& exchangeStdVec, // input/output
333  BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl) // output
334 {
335  if (m_env.inter0Rank() < 0) return;
336 
337  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
338  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
339  << ", level " << m_currLevel+LEVEL_REF_ID
340  << ", step " << m_currStep
341  << std::endl;
342  }
343 
344  unsigned int Np = (unsigned int) m_env.inter0Comm().NumProc();
345  if (m_env.inter0Rank() == 0) {
346  switch (currOptions->m_loadBalanceAlgorithmId) {
347  case 2:
348  justBalance_proc0(currOptions, // input
349  exchangeStdVec); // input/output
350  break;
351 
352  case 1:
353  default:
354 #ifdef QUESO_HAS_GLPK
355  // Get final node responsible for a linked chain by solving BIP at node zero only
356  solveBIP_proc0(exchangeStdVec); // input/output
357 #else
358  if (m_env.subDisplayFile()) {
359  *m_env.subDisplayFile() << "WARNING in MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
360  << ": algorithm id '" << currOptions->m_loadBalanceAlgorithmId
361  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
362  << ". Code will therefore process the algorithm id '" << 2
363  << "' instead..."
364  << std::endl;
365  }
366  if (m_env.subRank() == 0) {
367  std::cerr << "WARNING in MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()"
368  << ": algorithm id '" << currOptions->m_loadBalanceAlgorithmId
369  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
370  << ". Code will therefore process the algorithm id '" << 2
371  << "' instead..."
372  << std::endl;
373  }
374  justBalance_proc0(currOptions, // input
375  exchangeStdVec); // input/output
376 #endif
377  break;
378  }
379  } // if (m_env.inter0Rank() == 0)
380 
381  m_env.inter0Comm().Barrier();
382 
384  // Proc 0 now broadcasts the information on 'exchangeStdVec'
386  unsigned int exchangeStdVecSize = exchangeStdVec.size();
387  m_env.inter0Comm().Bcast((void *) &exchangeStdVecSize, (int) 1, RawValue_MPI_UNSIGNED, 0, // LOAD BALANCE
388  "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
389  "failed MPI.Bcast() for exchangeStdVec size");
390  if (m_env.inter0Rank() > 0) exchangeStdVec.resize(exchangeStdVecSize);
391 
392  m_env.inter0Comm().Bcast((void *) &exchangeStdVec[0], (int) (exchangeStdVecSize*sizeof(ExchangeInfoStruct)), RawValue_MPI_CHAR, 0, // LOAD BALANCE
393  "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
394  "failed MPI.Bcast() for exchangeStdVec data");
395 
397  // All "management" nodes update 'finalNumChainsPerNode' and 'finalNumPostionsPerNode'
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;
406  }
407 
409  // Sanity check
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);
414  //std::cout << m_env.worldRank() << ", finalRatioOfPosPerNode = " << finalRatioOfPosPerNode << std::endl;
415 
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, // LOAD BALANCE
420  "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
421  "failed MPI.Allreduce() for min");
422  //std::cout << m_env.worldRank() << ", minRatio = " << minRatio << std::endl;
423  queso_require_equal_to_msg(minRatio, finalRatioOfPosPerNode, "failed minRatio sanity check");
424 
425  double maxRatio = 0.;
426  auxBuf[0] = finalRatioOfPosPerNode;
427  m_env.inter0Comm().template Allreduce<double>(&auxBuf[0], &maxRatio, (int) auxBuf.size(), RawValue_MPI_MAX, // LOAD BALANCE
428  "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
429  "failed MPI.Allreduce() for max");
430  //std::cout << m_env.worldRank() << ", maxRatio = " << maxRatio << std::endl;
431  queso_require_equal_to_msg(maxRatio, finalRatioOfPosPerNode, "failed maxRatio sanity check");
432 
434  // Proc 0 now broadcasts the information on 'finalNumChainsPerNode'
436  unsigned int finalNumChainsPerNodeSize = finalNumChainsPerNode.size();
437  m_env.inter0Comm().Bcast((void *) &finalNumChainsPerNodeSize, (int) 1, RawValue_MPI_UNSIGNED, 0, // LOAD BALANCE
438  "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
439  "failed MPI.Bcast() for finalNumChainsPerNode size");
440  if (m_env.inter0Rank() > 0) finalNumChainsPerNode.resize(finalNumChainsPerNodeSize);
441 
442  m_env.inter0Comm().Bcast((void *) &finalNumChainsPerNode[0], (int) finalNumChainsPerNodeSize, RawValue_MPI_UNSIGNED, 0, // LOAD BALANCE
443  "MLSampling<P_V,P_M>::prepareBalLinkedChains_inter0()",
444  "failed MPI.Bcast() for finalNumChainsPerNode data");
445 
447  // Mpi exchange information between nodes and properly populate
448  // balancedLinkControl.linkedChains at each node
450  mpiExchangePositions_inter0(prevChain,
451  prevExponent,
452  currExponent,
453  prevLogLikelihoodValues,
454  prevLogTargetValues,
455  exchangeStdVec,
456  finalNumChainsPerNode,
457  finalNumPositionsPerNode, // It is already valid at all "management" nodes (not only at node 0) because of the sanity check above
458  balancedLinkControl);
459 
460  return;
461 }
462 
463 template <class P_V,class P_M>
464 void
465 MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0(
466  unsigned int indexOfFirstWeight, // input
467  unsigned int indexOfLastWeight, // input
468  const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // input
469  UnbalancedLinkedChainsPerNodeStruct& unbalancedLinkControl) // output
470 {
471  if (m_env.inter0Rank() < 0) return;
472 
473  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
474  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
475  << ", level " << m_currLevel+LEVEL_REF_ID
476  << ", step " << m_currStep
477  << ": indexOfFirstWeight = " << indexOfFirstWeight
478  << ", indexOfLastWeight = " << indexOfLastWeight
479  << std::endl;
480  }
481 
482  unsigned int subNumSamples = 0;
483  std::vector<unsigned int> unifiedIndexCountersAtAllProcs(0);
484 
485  // All nodes in 'inter0Comm' should resize to the same size // KAUST3
486  unsigned int resizeSize = unifiedIndexCountersAtProc0Only.size();
487  m_env.inter0Comm().Bcast((void *) &resizeSize, (int) 1, RawValue_MPI_UNSIGNED, 0,
488  "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
489  "failed MPI.Bcast() for resizeSize");
490  unifiedIndexCountersAtAllProcs.resize(resizeSize,0);
491 
492  if (m_env.inter0Rank() == 0) unifiedIndexCountersAtAllProcs = unifiedIndexCountersAtProc0Only;
493 
494  // Broadcast index counters to all nodes
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()"
501  << ", level " << m_currLevel+LEVEL_REF_ID
502  << ", step " << m_currStep
503  << ":"
504  << std::endl;
505  for (int r = 0; r < m_env.inter0Comm().NumProc(); ++r) {
506  *m_env.subDisplayFile() << " unifiedIndexCountersAtAllProcs[" << r << "] = " << unifiedIndexCountersAtAllProcs[r]
507  << std::endl;
508  }
509  }
510 #endif
511  //for (unsigned int i = 0; i < unifiedIndexCountersAtAllProcs.size(); ++i) {
512  // *m_env.subDisplayFile() << "unifiedIndexCountersAtAllProcs[" << i
513  // << "] = " << unifiedIndexCountersAtAllProcs[i]
514  // << std::endl;
515  //}
516 
517  // Use 'indexOfFirstWeight' and 'indexOfLastWeight' in order to update 'subNumSamples'
518  queso_require_less_msg(indexOfFirstWeight, unifiedIndexCountersAtAllProcs.size(), "invalid indexOfFirstWeight");
519  queso_require_less_msg(indexOfLastWeight, unifiedIndexCountersAtAllProcs.size(), "invalid indexOfLastWeight");
520  subNumSamples = 0;
521  for (unsigned int i = indexOfFirstWeight; i <= indexOfLastWeight; ++i) {
522  subNumSamples += unifiedIndexCountersAtAllProcs[i];
523  }
524 
525  std::vector<unsigned int> auxBuf(1,0);
526 
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");
532 
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");
538 
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");
544 
545 
546  // m_env.worldRank(),
547  // "MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()",
548  // "invalid state");
549 
550  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
551  *m_env.subDisplayFile() << "KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
552  << ", level " << m_currLevel+LEVEL_REF_ID
553  << ", step " << m_currStep
554  << ": subNumSamples = " << subNumSamples
555  << ", unifiedIndexCountersAtAllProcs.size() = " << unifiedIndexCountersAtAllProcs.size()
556  << std::endl;
557  *m_env.subDisplayFile() << "KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
558  << ", level " << m_currLevel+LEVEL_REF_ID
559  << ", step " << m_currStep
560  << ": minModifiedSubNumSamples = " << minModifiedSubNumSamples
561  << ", avgModifiedSubNumSamples = " << ((double) sumModifiedSubNumSamples)/((double) m_env.inter0Comm().NumProc())
562  << ", maxModifiedSubNumSamples = " << maxModifiedSubNumSamples
563  << std::endl;
564  }
565 
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()"
569  << ", level " << m_currLevel+LEVEL_REF_ID
570  << ", step " << m_currStep
571  << ": numberOfPositionsToGuaranteeForNode = " << numberOfPositionsToGuaranteeForNode
572  << std::endl;
573  }
574  for (unsigned int i = indexOfFirstWeight; i <= indexOfLastWeight; ++i) {
575 //for (unsigned int i = 0; i < unifiedIndexCountersAtAllProcs.size(); ++i) { // KAUST4: important
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]
581  << std::endl;
582  }
583  if (unifiedIndexCountersAtAllProcs[i] < numberOfPositionsToGuaranteeForNode) {
584  UnbalancedLinkedChainControlStruct auxControl;
585  auxControl.initialPositionIndexInPreviousChain = i;
586  auxControl.numberOfPositions = unifiedIndexCountersAtAllProcs[i];
587  unbalancedLinkControl.unbLinkedChains.push_back(auxControl);
588 
589  numberOfPositionsToGuaranteeForNode -= unifiedIndexCountersAtAllProcs[i];
590  unifiedIndexCountersAtAllProcs[i] = 0;
591  }
592  else if ((unifiedIndexCountersAtAllProcs[i] == numberOfPositionsToGuaranteeForNode) &&
593  (unifiedIndexCountersAtAllProcs[i] > 0 )) {
594  //else { // KAUST4
595  UnbalancedLinkedChainControlStruct auxControl;
596  auxControl.initialPositionIndexInPreviousChain = i;
597  auxControl.numberOfPositions = numberOfPositionsToGuaranteeForNode;
598  unbalancedLinkControl.unbLinkedChains.push_back(auxControl);
599 
600  unifiedIndexCountersAtAllProcs[i] -= numberOfPositionsToGuaranteeForNode;
601  numberOfPositionsToGuaranteeForNode = 0;
602  }
603  else if ((unifiedIndexCountersAtAllProcs[i] == numberOfPositionsToGuaranteeForNode) &&
604  (unifiedIndexCountersAtAllProcs[i] == 0 )) {
605  // Ok
606  }
607  else {
608  queso_error_msg("should never get here");
609  }
610  }
611  }
612  queso_require_equal_to_msg(numberOfPositionsToGuaranteeForNode, 0, "numberOfPositionsToGuaranteeForNode exited loop with wrong value");
613  // FIX ME: swap trick to save memory
614 
615  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
616  *m_env.subDisplayFile() << "KEY Leaving MLSampling<P_V,P_M>::prepareUnbLinkedChains_inter0()"
617  << ", level " << m_currLevel+LEVEL_REF_ID
618  << ", step " << m_currStep
619  << ": unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.unbLinkedChains.size()
620  << std::endl;
621  }
622 
623  return;
624 }
625 
626 template <class P_V,class P_M>
627 void
628 MLSampling<P_V,P_M>::generateBalLinkedChains_all( // EXTRA FOR LOAD BALANCE
629  MLSamplingLevelOptions& inputOptions, // input, only m_rawChainSize changes
630  const P_M& unifiedCovMatrix, // input
631  const GenericVectorRV <P_V,P_M>& rv, // input
632  const BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl, // input // Round Rock
633  SequenceOfVectors <P_V,P_M>& workingChain, // output
634  double& cumulativeRunTime, // output
635  unsigned int& cumulativeRejections, // output
636  ScalarSequence <double>* currLogLikelihoodValues, // output
637  ScalarSequence <double>* currLogTargetValues) // output
638 {
639  m_env.fullComm().Barrier();
640 
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()
644  << std::endl;
645  }
646 
647  P_V auxInitialPosition(m_vectorSpace.zeroVector());
648  double auxInitialLogPrior;
649  double auxInitialLogLikelihood;
650 
651  unsigned int chainIdMax = 0;
652  if (m_env.inter0Rank() >= 0) {
653  chainIdMax = balancedLinkControl.balLinkedChains.size();
654  }
655  // KAUST: all nodes in 'subComm' should have the same 'chainIdMax'
656  m_env.subComm().Bcast((void *) &chainIdMax, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'subComm', important // LOAD BALANCE
657  "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
658  "failed MPI.Bcast() for chainIdMax");
659 
660  struct timeval timevalEntering;
661  int iRC = 0;
662  iRC = gettimeofday(&timevalEntering, NULL);
663  if (iRC) {}; // just to remove compiler warning
664 
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;
669  }
670 
671  std::vector<unsigned int> auxBuf(1,0);
672 
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, // LOAD BALANCE
676  "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
677  "failed MPI.Allreduce() for min");
678 
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, // LOAD BALANCE
682  "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
683  "failed MPI.Allreduce() for max");
684 
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, // LOAD BALANCE
688  "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
689  "failed MPI.Allreduce() for sum");
690 
691  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
692  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
693  << ", level " << m_currLevel+LEVEL_REF_ID
694  << ", step " << m_currStep
695  << ": chainIdMax = " << chainIdMax
696  << ", numberOfPositions = " << numberOfPositions
697  << ", at " << ctime(&timevalEntering.tv_sec)
698  << std::endl;
699  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
700  << ", level " << m_currLevel+LEVEL_REF_ID
701  << ", step " << m_currStep
702  << ": minNumberOfPositions = " << minNumberOfPositions
703  << ", avgNumberOfPositions = " << ((double) sumNumberOfPositions)/((double) m_env.inter0Comm().NumProc())
704  << ", maxNumberOfPositions = " << maxNumberOfPositions
705  << std::endl;
706  }
707 
708  // 2013-02-23: print sizes, and expected final size
709  }
710  if ((m_debugExponent == 1.) &&
711  (m_currStep == 10)) {
712  //m_env.setExceptionalCircumstance(true);
713  }
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) {
718  // aqui 4
719  auxInitialPosition = *(balancedLinkControl.balLinkedChains[chainId].initialPosition); // Round Rock
720  auxInitialLogPrior = balancedLinkControl.balLinkedChains[chainId].initialLogPrior;
721  auxInitialLogLikelihood = balancedLinkControl.balLinkedChains[chainId].initialLogLikelihood;
722  tmpChainSize = balancedLinkControl.balLinkedChains[chainId].numberOfPositions+1; // IMPORTANT: '+1' in order to discard initial position afterwards
723  if ((m_env.subDisplayFile() ) &&
724  (m_env.displayVerbosity() >= 3)) {
725  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
726  << ", level " << m_currLevel+LEVEL_REF_ID
727  << ", step " << m_currStep
728  << ", chainId = " << chainId
729  << " < " << chainIdMax
730  << ": begin generating " << tmpChainSize
731  << " chain positions"
732  << std::endl;
733  }
734  }
735  auxInitialPosition.mpiBcast(0, m_env.subComm()); // Yes, 'subComm', important // KAUST
736 
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
742  << std::endl;
743  }
744  m_env.subComm().Barrier();
745  }
746  sleep(1);
747 #endif
748 
749  // KAUST: all nodes in 'subComm' should have the same 'tmpChainSize'
750  m_env.subComm().Bcast((void *) &tmpChainSize, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'subComm', important // LOAD BALANCE
751  "MLSampling<P_V,P_M>::generateBalLinkedChains_all()",
752  "failed MPI.Bcast() for tmpChainSize");
753 
754  inputOptions.m_rawChainSize = tmpChainSize;
755  SequenceOfVectors<P_V,P_M> tmpChain(m_vectorSpace,
756  0,
757  m_options.m_prefix+"tmp_chain");
758  ScalarSequence<double> tmpLogLikelihoodValues(m_env,0,"");
759  ScalarSequence<double> tmpLogTargetValues (m_env,0,"");
760 
761  MHRawChainInfoStruct mcRawInfo;
762  if (inputOptions.m_initialPositionUsePreviousLevelLikelihood) { // ml_likelihood_caching
763  m_env.subComm().Bcast((void *) &auxInitialLogPrior, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm', important
764  "MLSamplingClass<P_V,P_M>::generateBalLinkedChains_all()",
765  "failed MPI.Bcast() for auxInitialLogPrior");
766  m_env.subComm().Bcast((void *) &auxInitialLogLikelihood, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm', important
767  "MLSamplingClass<P_V,P_M>::generateBalLinkedChains_all()",
768  "failed MPI.Bcast() for auxInitialLogLikelihood");
769 
770  MetropolisHastingsSG<P_V,P_M> mcSeqGenerator(inputOptions,
771  rv,
772  auxInitialPosition, // KEY new: pass logPrior and logLikelihood
773  auxInitialLogPrior,
774  auxInitialLogLikelihood,
775  &unifiedCovMatrix);
776  mcSeqGenerator.generateSequence(tmpChain,
777  &tmpLogLikelihoodValues, // likelihood is IMPORTANT
778  &tmpLogTargetValues);
779  mcSeqGenerator.getRawChainInfo(mcRawInfo);
780  }
781  else {
782  MetropolisHastingsSG<P_V,P_M> mcSeqGenerator(inputOptions,
783  rv,
784  auxInitialPosition,
785  &unifiedCovMatrix);
786  mcSeqGenerator.generateSequence(tmpChain,
787  &tmpLogLikelihoodValues, // likelihood is IMPORTANT
788  &tmpLogTargetValues);
789  mcSeqGenerator.getRawChainInfo(mcRawInfo);
790  }
791 
792  cumulativeRunTime += mcRawInfo.runTime;
793  cumulativeRejections += mcRawInfo.numRejections;
794 
795  if (m_env.inter0Rank() >= 0) {
796  if (m_env.exceptionalCircumstance()) {
797  if ((m_env.subDisplayFile() ) &&
798  (m_env.displayVerbosity() >= 0)) { // detailed output debug
799  P_V tmpVec(m_vectorSpace.zeroVector());
800  for (unsigned int i = 0; i < tmpLogLikelihoodValues.subSequenceSize(); ++i) {
801  tmpChain.getPositionValues(i,tmpVec);
802  *m_env.subDisplayFile() << "DEBUG finalChain[" << cumulativeNumPositions+i << "] "
803  << "= tmpChain[" << i << "] = " << tmpVec
804  << ", tmpLogLikelihoodValues[" << i << "] = " << tmpLogLikelihoodValues[i]
805  << ", tmpLogTargetValues[" << i << "] = " << tmpLogTargetValues[i]
806  << std::endl;
807  }
808  }
809  } // exceptional
810 
811  cumulativeNumPositions += tmpChainSize;
812  if (cumulativeNumPositions > 100) m_env.setExceptionalCircumstance(false);
813 
814  if ((m_env.subDisplayFile() ) &&
815  (m_env.displayVerbosity() >= 3)) {
816  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
817  << ", level " << m_currLevel+LEVEL_REF_ID
818  << ", step " << m_currStep
819  << ", chainId = " << chainId
820  << " < " << chainIdMax
821  << ": finished generating " << tmpChain.subSequenceSize()
822  << " chain positions"
823  << std::endl;
824  }
825 
826  // KAUST5: what if workingChain ends up with different size in different nodes? Important
827  workingChain.append (tmpChain, 1,tmpChain.subSequenceSize()-1 ); // IMPORTANT: '1' in order to discard initial position
828  if (currLogLikelihoodValues) {
829  currLogLikelihoodValues->append(tmpLogLikelihoodValues,1,tmpLogLikelihoodValues.subSequenceSize()-1); // IMPORTANT: '1' in order to discard initial position
830  if ((m_env.subDisplayFile() ) &&
831  (m_env.displayVerbosity() >= 99) &&
832  (chainId == 0 )) {
833  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
834  << ", level " << m_currLevel+LEVEL_REF_ID
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]
841  << std::endl;
842  }
843  }
844  if (currLogTargetValues) {
845  currLogTargetValues->append (tmpLogTargetValues, 1,tmpLogTargetValues.subSequenceSize()-1 ); // IMPORTANT: '1' in order to discard initial position
846  }
847  // 2013-02-23: print size just appended
848  }
849  } // for 'chainId'
850 
851  // 2013-02-23: print final size
852 
853  struct timeval timevalBarrier;
854  iRC = gettimeofday(&timevalBarrier, NULL);
855  if (iRC) {}; // just to remove compiler warning
856  double loopTime = MiscGetEllapsedSeconds(&timevalEntering);
857  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
858  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
859  << ", level " << m_currLevel+LEVEL_REF_ID
860  << ", step " << m_currStep
861  << ": ended chain loop after " << loopTime << " seconds"
862  << ", calling fullComm().Barrier() at " << ctime(&timevalBarrier.tv_sec)
863  << std::endl;
864  }
865 
866  m_env.fullComm().Barrier(); // KAUST4
867 
868  struct timeval timevalLeaving;
869  iRC = gettimeofday(&timevalLeaving, NULL);
870  if (iRC) {}; // just to remove compiler warning
871  double barrierTime = MiscGetEllapsedSeconds(&timevalBarrier);
872  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
873  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateBalLinkedChains_all()"
874  << ", level " << m_currLevel+LEVEL_REF_ID
875  << ", step " << m_currStep
876  << ": after " << barrierTime << " seconds in fullComm().Barrier()"
877  << ", at " << ctime(&timevalLeaving.tv_sec)
878  << std::endl;
879  }
880 
881  return;
882 }
883 
884 template <class P_V,class P_M>
885 void
886 MLSampling<P_V,P_M>::generateUnbLinkedChains_all(
887  MLSamplingLevelOptions& inputOptions, // input, only m_rawChainSize changes
888  const P_M& unifiedCovMatrix, // input
889  const GenericVectorRV <P_V,P_M>& rv, // input
890  const UnbalancedLinkedChainsPerNodeStruct& unbalancedLinkControl, // input // Round Rock
891  unsigned int indexOfFirstWeight, // input // Round Rock
892  const SequenceOfVectors<P_V,P_M>& prevChain, // input // Round Rock
893  double prevExponent, // input
894  double currExponent, // input
895  const ScalarSequence<double>& prevLogLikelihoodValues, // input
896  const ScalarSequence<double>& prevLogTargetValues, // input
897  SequenceOfVectors <P_V,P_M>& workingChain, // output
898  double& cumulativeRunTime, // output
899  unsigned int& cumulativeRejections, // output
900  ScalarSequence <double>* currLogLikelihoodValues, // output
901  ScalarSequence <double>* currLogTargetValues) // output
902 {
903  m_env.fullComm().Barrier();
904 
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
909  << std::endl;
910  }
911 
912  P_V auxInitialPosition(m_vectorSpace.zeroVector());
913  double auxInitialLogPrior;
914  double auxInitialLogLikelihood;
915 
916  unsigned int chainIdMax = 0;
917  if (m_env.inter0Rank() >= 0) {
918  chainIdMax = unbalancedLinkControl.unbLinkedChains.size();
919  }
920  // KAUST: all nodes in 'subComm' should have the same 'chainIdMax'
921  m_env.subComm().Bcast((void *) &chainIdMax, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'subComm', important
922  "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
923  "failed MPI.Bcast() for chainIdMax");
924 
925  struct timeval timevalEntering;
926  int iRC = 0;
927  iRC = gettimeofday(&timevalEntering, NULL);
928  if (iRC) {}; // just to remove compiler warning
929 
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;
934  }
935 
936  std::vector<unsigned int> auxBuf(1,0);
937 
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");
943 
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");
949 
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");
955 
956  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
957  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
958  << ", level " << m_currLevel+LEVEL_REF_ID
959  << ", step " << m_currStep
960  << ": chainIdMax = " << chainIdMax
961  << ", numberOfPositions = " << numberOfPositions
962  << ", at " << ctime(&timevalEntering.tv_sec)
963  << std::endl;
964  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
965  << ", level " << m_currLevel+LEVEL_REF_ID
966  << ", step " << m_currStep
967  << ": minNumberOfPositions = " << minNumberOfPositions
968  << ", avgNumberOfPositions = " << ((double) sumNumberOfPositions)/((double) m_env.inter0Comm().NumProc())
969  << ", maxNumberOfPositions = " << maxNumberOfPositions
970  << std::endl;
971  }
972  }
973  if ((m_debugExponent == 1.) &&
974  (m_currStep == 10)) {
975  //m_env.setExceptionalCircumstance(true);
976  }
977  double expRatio = currExponent;
978  if (prevExponent > 0.0) {
979  expRatio /= prevExponent;
980  }
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; // KAUST4 // Round Rock
986  prevChain.getPositionValues(auxIndex,auxInitialPosition); // Round Rock
987  auxInitialLogPrior = prevLogTargetValues[auxIndex] - prevLogLikelihoodValues[auxIndex];
988  auxInitialLogLikelihood = expRatio * prevLogLikelihoodValues[auxIndex];
989  tmpChainSize = unbalancedLinkControl.unbLinkedChains[chainId].numberOfPositions+1; // IMPORTANT: '+1' in order to discard initial position afterwards
990  if ((m_env.subDisplayFile() ) &&
991  (m_env.displayVerbosity() >= 3)) {
992  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
993  << ", level " << m_currLevel+LEVEL_REF_ID
994  << ", step " << m_currStep
995  << ", chainId = " << chainId
996  << " < " << chainIdMax
997  << ": begin generating " << tmpChainSize
998  << " chain positions"
999  << std::endl;
1000  }
1001  }
1002  auxInitialPosition.mpiBcast(0, m_env.subComm()); // Yes, 'subComm', important // KAUST
1003 
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
1009  << std::endl;
1010  }
1011  m_env.subComm().Barrier();
1012  }
1013  sleep(1);
1014 #endif
1015 
1016  // KAUST: all nodes in 'subComm' should have the same 'tmpChainSize'
1017  m_env.subComm().Bcast((void *) &tmpChainSize, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'subComm', important
1018  "MLSampling<P_V,P_M>::generateUnbLinkedChains_all()",
1019  "failed MPI.Bcast() for tmpChainSize");
1020 
1021  inputOptions.m_rawChainSize = tmpChainSize;
1022  SequenceOfVectors<P_V,P_M> tmpChain(m_vectorSpace,
1023  0,
1024  m_options.m_prefix+"tmp_chain");
1025  ScalarSequence<double> tmpLogLikelihoodValues(m_env,0,"");
1026  ScalarSequence<double> tmpLogTargetValues (m_env,0,"");
1027 
1028  // KAUST: all nodes should call here
1029  MHRawChainInfoStruct mcRawInfo;
1030  if (inputOptions.m_initialPositionUsePreviousLevelLikelihood) { // ml_likelihood_caching
1031  m_env.subComm().Bcast((void *) &auxInitialLogPrior, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm', important
1032  "MLSamplingClass<P_V,P_M>::generateUnbLinkedChains_all()",
1033  "failed MPI.Bcast() for auxInitialLogPrior");
1034  m_env.subComm().Bcast((void *) &auxInitialLogLikelihood, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm', important
1035  "MLSamplingClass<P_V,P_M>::generateUnbLinkedChains_all",
1036  "failed MPI.Bcast() for auxInitialLogLikelihood");
1037  MetropolisHastingsSG<P_V,P_M> mcSeqGenerator(inputOptions,
1038  rv,
1039  auxInitialPosition, // KEY new: pass logPrior and logLikelihood
1040  auxInitialLogPrior,
1041  auxInitialLogLikelihood,
1042  &unifiedCovMatrix);
1043  mcSeqGenerator.generateSequence(tmpChain,
1044  &tmpLogLikelihoodValues, // likelihood is IMPORTANT
1045  &tmpLogTargetValues);
1046  mcSeqGenerator.getRawChainInfo(mcRawInfo);
1047  }
1048  else {
1049  MetropolisHastingsSG<P_V,P_M> mcSeqGenerator(inputOptions,
1050  rv,
1051  auxInitialPosition,
1052  &unifiedCovMatrix);
1053  mcSeqGenerator.generateSequence(tmpChain,
1054  &tmpLogLikelihoodValues, // likelihood is IMPORTANT
1055  &tmpLogTargetValues);
1056  mcSeqGenerator.getRawChainInfo(mcRawInfo);
1057  }
1058 
1059  cumulativeRunTime += mcRawInfo.runTime;
1060  cumulativeRejections += mcRawInfo.numRejections;
1061 
1062  if (m_env.inter0Rank() >= 0) {
1063  if (m_env.exceptionalCircumstance()) {
1064  if ((m_env.subDisplayFile() ) &&
1065  (m_env.displayVerbosity() >= 0)) { // detailed output debug
1066  P_V tmpVec(m_vectorSpace.zeroVector());
1067  for (unsigned int i = 0; i < tmpLogLikelihoodValues.subSequenceSize(); ++i) {
1068  tmpChain.getPositionValues(i,tmpVec);
1069  *m_env.subDisplayFile() << "DEBUG finalChain[" << cumulativeNumPositions+i << "] "
1070  << "= tmpChain[" << i << "] = " << tmpVec
1071  << ", tmpLogLikelihoodValues[" << i << "] = " << tmpLogLikelihoodValues[i]
1072  << ", tmpLogTargetValues[" << i << "] = " << tmpLogTargetValues[i]
1073  << std::endl;
1074  }
1075  }
1076  } // exceptional
1077 
1078  cumulativeNumPositions += tmpChainSize;
1079  if (cumulativeNumPositions > 100) m_env.setExceptionalCircumstance(false);
1080 
1081  if ((m_env.subDisplayFile() ) &&
1082  (m_env.displayVerbosity() >= 3)) {
1083  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1084  << ", level " << m_currLevel+LEVEL_REF_ID
1085  << ", step " << m_currStep
1086  << ", chainId = " << chainId
1087  << " < " << chainIdMax
1088  << ": finished generating " << tmpChain.subSequenceSize()
1089  << " chain positions"
1090  << std::endl;
1091  }
1092 
1093  // KAUST5: what if workingChain ends up with different size in different nodes? Important
1094  workingChain.append (tmpChain, 1,tmpChain.subSequenceSize()-1 ); // IMPORTANT: '1' in order to discard initial position
1095  if (currLogLikelihoodValues) {
1096  currLogLikelihoodValues->append(tmpLogLikelihoodValues,1,tmpLogLikelihoodValues.subSequenceSize()-1); // IMPORTANT: '1' in order to discard initial position
1097  if ((m_env.subDisplayFile() ) &&
1098  (m_env.displayVerbosity() >= 99) &&
1099  (chainId == 0 )) {
1100  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1101  << ", level " << m_currLevel+LEVEL_REF_ID
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]
1108  << std::endl;
1109  }
1110  }
1111  if (currLogTargetValues) {
1112  currLogTargetValues->append (tmpLogTargetValues, 1,tmpLogTargetValues.subSequenceSize()-1 ); // IMPORTANT: '1' in order to discard initial position
1113  }
1114  }
1115  } // for 'chainId'
1116 
1117  struct timeval timevalBarrier;
1118  iRC = gettimeofday(&timevalBarrier, NULL);
1119  if (iRC) {}; // just to remove compiler warning
1120  double loopTime = MiscGetEllapsedSeconds(&timevalEntering);
1121  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1122  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1123  << ", level " << m_currLevel+LEVEL_REF_ID
1124  << ", step " << m_currStep
1125  << ": ended chain loop after " << loopTime << " seconds"
1126  << ", calling fullComm().Barrier() at " << ctime(&timevalBarrier.tv_sec)
1127  << std::endl;
1128  }
1129 
1130  m_env.fullComm().Barrier(); // KAUST4
1131 
1132  struct timeval timevalLeaving;
1133  iRC = gettimeofday(&timevalLeaving, NULL);
1134  if (iRC) {}; // just to remove compiler warning
1135  double barrierTime = MiscGetEllapsedSeconds(&timevalBarrier);
1136  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1137  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateUnbLinkedChains_all()"
1138  << ", level " << m_currLevel+LEVEL_REF_ID
1139  << ", step " << m_currStep
1140  << ": after " << barrierTime << " seconds in fullComm().Barrier()"
1141  << ", at " << ctime(&timevalLeaving.tv_sec)
1142  << std::endl;
1143  }
1144 
1145  return;
1146 }
1147 
1148 #ifdef QUESO_HAS_GLPK
1149 template <class P_V,class P_M>
1150 void
1151 MLSampling<P_V,P_M>::solveBIP_proc0( // EXTRA FOR LOAD BALANCE
1152  std::vector<ExchangeInfoStruct>& exchangeStdVec) // input/output
1153 {
1154  if (m_env.inter0Rank() != 0) return;
1155 
1156  int iRC = UQ_OK_RC;
1157  struct timeval timevalBIP;
1158  iRC = gettimeofday(&timevalBIP, NULL);
1159  if (iRC) {}; // just to remove compiler warning
1160 
1161  unsigned int Np = (unsigned int) m_env.inter0Comm().NumProc();
1162  unsigned int Nc = exchangeStdVec.size();
1163 
1164  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1165  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::solveBIP_proc0()"
1166  << ", level " << m_currLevel+LEVEL_REF_ID
1167  << ", step " << m_currStep
1168  << ": Np = " << Np
1169  << ", Nc = " << Nc
1170  << std::endl;
1171  }
1172 
1174  // Instantiate BIP
1176  glp_prob *lp;
1177  lp = glp_create_prob();
1178  glp_set_prob_name(lp, "sample");
1179 
1181  // Set rows and colums of BIP constraint matrix
1183  unsigned int m = Nc+Np-1;
1184  unsigned int n = Nc*Np;
1185  unsigned int ne = Nc*Np + 2*Nc*(Np -1);
1186 
1187  glp_add_rows(lp, m); // Not 'm+1'
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, "");
1191  }
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, "");
1195  }
1196 
1197  glp_add_cols(lp, n); // Not 'n+1'
1198  for (int j = 1; j <= (int) n; ++j) {
1199  //glp_set_col_kind(lp, j, GLP_BV);
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, "");
1203  }
1204 
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);
1208  }
1209 
1210  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1211  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1212  << ", level " << m_currLevel+LEVEL_REF_ID
1213  << ", step " << m_currStep
1214  << ": finished setting BIP rows and cols"
1215  << ", m = " << m
1216  << ", n = " << n
1217  << ", ne = " << ne
1218  << std::endl;
1219  }
1220 
1222  // Load constraint matrix
1224  std::vector<int > iVec(ne+1,0);
1225  std::vector<int > jVec(ne+1,0);
1226  std::vector<double> aVec(ne+1,0.);
1227  int coefId = 1; // Yes, '1'
1228  for (int i = 1; i <= (int) Nc; ++i) {
1229  for (int j = 1; j <= (int) Np; ++j) {
1230  iVec[coefId] = i;
1231  jVec[coefId] = (i-1)*Np + j;
1232  aVec[coefId] = 1.;
1233  coefId++;
1234  }
1235  }
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);
1241  coefId++;
1242 
1243  iVec[coefId] = Nc+i;
1244  jVec[coefId] = (j-1)*Np + i + 1;
1245  aVec[coefId] = exchangeStdVec[j-1].numberOfPositions;
1246  coefId++;
1247  }
1248  }
1249  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1250  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1251  << ", level " << m_currLevel+LEVEL_REF_ID
1252  << ", step " << m_currStep
1253  << ": finished setting BIP constraint matrix"
1254  << ", ne = " << ne
1255  << ", coefId = " << coefId
1256  << std::endl;
1257  }
1258 
1259  queso_require_equal_to_msg(coefId, (int) (ne+1), "invalid final coefId");
1260 
1261  glp_load_matrix(lp, ne, &iVec[0], &jVec[0], &aVec[0]);
1262 
1263  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1264  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1265  << ", level " << m_currLevel+LEVEL_REF_ID
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)
1273  << std::endl;
1274  }
1275 
1277  // Check BIP before solving it
1279  queso_require_equal_to_msg(glp_get_num_rows(lp), (int) m, "invalid number of rows");
1280 
1281  queso_require_equal_to_msg(glp_get_num_cols(lp), (int) n, "invalid number of columnss");
1282 
1283  queso_require_equal_to_msg(glp_get_num_nz(lp), (int) ne, "invalid number of nonzero constraint coefficients");
1284 
1285  queso_require_equal_to_msg(glp_get_num_int(lp), (int) n, "invalid number of integer structural variables");
1286 
1287  queso_require_equal_to_msg(glp_get_num_bin(lp), (int) n, "invalid number of binary structural variables");
1288 
1290  // Set initial state
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;
1295  if (nodeId == 0) {
1296  glp_set_col_stat(lp, j, GLP_BS);
1297  }
1298  else {
1299  glp_set_col_stat(lp, j, GLP_BS);
1300  }
1301  }
1302  }
1303 #if 0
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);
1308  if (nodeId == 0) {
1309  queso_require_equal_to_msg(initialState, 1, "for nodeId = 0, initial state should be '1'");
1310  }
1311  else {
1312  queso_require_equal_to_msg(initialState, 0, "for nodeId > 0, initial state should be '0'");
1313  }
1314  }
1315  }
1316 #endif
1317  for (int i = 1; i <= (int) Nc; ++i) {
1318  glp_set_row_stat(lp, i, GLP_NS);
1319  }
1320  for (int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1321  glp_set_row_stat(lp, i, GLP_BS);
1322  }
1323 
1324  //glp_write_mps(lp, GLP_MPS_DECK, NULL, "nada.fixed_mps");
1325  //glp_write_mps(lp, GLP_MPS_FILE, NULL, "nada.free_mps" );
1326  //glp_write_lp (lp, NULL, "nada.cplex");
1327 
1329  // Solve BIP
1331  BIP_routine_struct BIP_routine_info;
1332  BIP_routine_info.env = &m_env;
1333  BIP_routine_info.currLevel = m_currLevel;
1334 
1335  glp_iocp BIP_params;
1336  glp_init_iocp(&BIP_params);
1337  BIP_params.presolve = GLP_ON;
1338  // aqui 2
1339  //BIP_params.binarize = GLP_ON;
1340  //BIP_params.cb_func = BIP_routine;
1341  //BIP_params.cb_info = (void *) (&BIP_routine_info);
1342  int BIP_rc = glp_intopt(lp, &BIP_params);
1343 
1344  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1345  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1346  << ", level " << m_currLevel+LEVEL_REF_ID
1347  << ", step " << m_currStep
1348  << ": finished solving BIP"
1349  << ", BIP_rc = " << BIP_rc
1350  << std::endl;
1351  }
1352 
1353  queso_require_equal_to_msg(BIP_rc, 0, "BIP returned rc != 0");
1354 
1356  // Check BIP status after solution
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()"
1361  << ", level " << m_currLevel+LEVEL_REF_ID
1362  << ", step " << m_currStep
1363  << ": BIP_Status = " << BIP_Status
1364  << std::endl;
1365  }
1366 
1367  switch (BIP_Status) {
1368  case GLP_OPT:
1369  // Ok
1370  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1371  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1372  << ", level " << m_currLevel+LEVEL_REF_ID
1373  << ", step " << m_currStep
1374  << ": BIP solution is optimal"
1375  << std::endl;
1376  }
1377  break;
1378 
1379  case GLP_FEAS:
1380  // Ok
1381  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1382  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1383  << ", level " << m_currLevel+LEVEL_REF_ID
1384  << ", step " << m_currStep
1385  << ": BIP solution is guaranteed to be 'only' feasible"
1386  << std::endl;
1387  }
1388  break;
1389 
1390  default:
1391  queso_error_msg("BIP has an undefined solution or has no solution");
1392  break;
1393  }
1394 
1395  for (int i = 1; i <= (int) Nc; ++i) {
1396  queso_require_equal_to_msg(glp_mip_row_val(lp, i), 1, "row should have value 1 at solution");
1397  }
1398  for (int i = (Nc+1); i <= (int) (Nc+Np-1); ++i) {
1399  queso_require_less_equal_msg(glp_mip_row_val(lp, i), 0, "row should have value 0 or should be negative at solution");
1400  }
1401 
1403  // Prepare output information, needed to for MPI distribution afterwards
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) {
1411  // Do nothing
1412  }
1413  else if (glp_mip_col_val(lp, j) == 1) {
1414  queso_require_equal_to_msg(exchangeStdVec[chainId].finalNodeOfInitialPosition, -1, "chain has already been taken care of");
1415  exchangeStdVec[chainId].finalNodeOfInitialPosition = nodeId;
1416  finalNumChainsPerNode [nodeId] += 1;
1417  finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
1418  }
1419  else {
1420  queso_error_msg("control variable should be either '0' or '1'");
1421  }
1422  }
1423  }
1424 
1425  unsigned int finalMinPosPerNode = *std::min_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1426  unsigned int finalMaxPosPerNode = *std::max_element(finalNumPositionsPerNode.begin(), finalNumPositionsPerNode.end());
1427 
1428  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1429  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::solveBIP_proc0()"
1430  << ", level " << m_currLevel+LEVEL_REF_ID
1431  << ", step " << m_currStep
1432  << ": finished preparing output information"
1433  << std::endl;
1434  }
1435 
1437  // Printout solution information
1439  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1440  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1441  << ", level " << m_currLevel+LEVEL_REF_ID
1442  << ", step " << m_currStep
1443  << ": solution gives the following redistribution"
1444  << std::endl;
1445  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1446  *m_env.subDisplayFile() << " KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1447  << ", level " << m_currLevel+LEVEL_REF_ID
1448  << ", step " << m_currStep
1449  << ", finalNumChainsPerNode[" << nodeId << "] = " << finalNumChainsPerNode[nodeId]
1450  << ", finalNumPositionsPerNode[" << nodeId << "] = " << finalNumPositionsPerNode[nodeId]
1451  << std::endl;
1452  }
1453  *m_env.subDisplayFile() << " KEY In MLSampling<P_V,P_M>::solveBIP_proc0()"
1454  << ", level " << m_currLevel+LEVEL_REF_ID
1455  << ", step " << m_currStep
1456  << ", finalRatioOfPosPerNode = " << ((double) finalMaxPosPerNode) / ((double)finalMinPosPerNode)
1457  << std::endl;
1458  }
1459 
1461  // Make sanity checks
1463  queso_require_equal_to_msg(glp_mip_obj_val(lp), (double) finalNumPositionsPerNode[0], "Invalid objective value");
1464 
1465  for (unsigned int nodeId = 1; nodeId < Np; ++nodeId) { // Yes, '1'
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");
1467  }
1468 
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]);
1472  queso_require_equal_to_msg(glp_mip_row_val(lp, i), diff, "wrong state");
1473  }
1474 
1476  // Free memory and return
1478  glp_delete_prob(lp);
1479 
1480  double bipRunTime = MiscGetEllapsedSeconds(&timevalBIP);
1481  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1482  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::solveBIP_proc0()"
1483  << ", level " << m_currLevel+LEVEL_REF_ID
1484  << ", step " << m_currStep
1485  << ", after " << bipRunTime << " seconds"
1486  << std::endl;
1487  }
1488 
1489  return;
1490 }
1491 #endif // QUESO_HAS_GLPK
1492 
1493 template <class P_V,class P_M>
1494 void
1495 MLSampling<P_V,P_M>::justBalance_proc0(
1496  const MLSamplingLevelOptions* currOptions, // input
1497  std::vector<ExchangeInfoStruct>& exchangeStdVec) // input/output
1498 {
1499  if (m_env.inter0Rank() != 0) return;
1500 
1501  int iRC = UQ_OK_RC;
1502  struct timeval timevalBal;
1503  iRC = gettimeofday(&timevalBal, NULL);
1504  if (iRC) {}; // just to remove compiler warning
1505 
1506  unsigned int Np = m_env.numSubEnvironments();
1507  unsigned int Nc = exchangeStdVec.size();
1508 
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; // final = original
1513  }
1514 
1516  // Compute original ratio of positions per node
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; // Yes, 'final'
1523  currNumChainsPerNode [nodeId] += 1;
1524  currNumPositionsPerNode[nodeId] += currExchangeStdVec[chainId].numberOfPositions;
1525  iterIdMax += currExchangeStdVec[chainId].numberOfPositions;
1526  }
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);
1530 
1532  // Loop
1534  //iterIdMax /= 2;
1535  int iterId = -1;
1536  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1537  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::justBalance_proc0()"
1538  << ", level " << m_currLevel+LEVEL_REF_ID
1539  << ", step " << m_currStep
1540  << ", iter " << iterId
1541  << ", currRatioOfPosPerNode = " << currRatioOfPosPerNode
1542  << std::endl;
1543  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1544  *m_env.subDisplayFile() << " KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1545  << ", level " << m_currLevel+LEVEL_REF_ID
1546  << ", step " << m_currStep
1547  << ", iter " << iterId
1548  << ", currNumChainsPerNode[" << nodeId << "] = " << currNumChainsPerNode[nodeId]
1549  << ", currNumPositionsPerNode[" << nodeId << "] = " << currNumPositionsPerNode[nodeId]
1550  << std::endl;
1551  }
1552  }
1553 
1554  std::vector<std::vector<double> > vectorOfChainSizesPerNode(Np);
1555  while ((iterId < (int) iterIdMax ) &&
1556  (currRatioOfPosPerNode > currOptions->m_loadBalanceTreshold)) {
1557  iterId++;
1558 
1560  // Initialize information
1562  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1563  vectorOfChainSizesPerNode[nodeId].clear(); // make sure vectors have size 0
1564  }
1565  for (unsigned int chainId = 0; chainId < Nc; ++chainId) {
1566  unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition; // Yes, 'final'
1567  vectorOfChainSizesPerNode[nodeId].push_back(currExchangeStdVec[chainId].numberOfPositions);
1568  }
1569  // FIX ME: swap to save memory
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");
1573  }
1574 
1576  // Find [node with most postions], [node with least positions] and [number of positions to move]
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;
1586  }
1587  if (currNumPositionsPerNode[nodeId] < currSmallestAmountOfPositionsPerNode) {
1588  currSmallestAmountOfPositionsPerNode = currNumPositionsPerNode[nodeId];
1589  currNodeWithLeastPositions = nodeId;
1590  }
1591  }
1592 
1593  queso_require_equal_to_msg(currMinPosPerNode, currNumPositionsPerNode[currNodeWithLeastPositions], "inconsistent currMinPosPerNode");
1594 
1595  queso_require_equal_to_msg(currMaxPosPerNode, currNumPositionsPerNode[currNodeWithMostPositions], "inconsistent currMaxPosPerNode");
1596 
1597  unsigned int numberOfPositionsToMove = vectorOfChainSizesPerNode[currNodeWithMostPositions][0];
1598 
1599  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1600  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::justBalance_proc0()"
1601  << ", level " << m_currLevel+LEVEL_REF_ID
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
1610  << std::endl;
1611  }
1612 
1614  // Update 'final' fields in the two nodes
1616  std::vector<ExchangeInfoStruct> newExchangeStdVec(Nc);
1617  for (unsigned int chainId = 0; chainId < Nc; ++chainId) {
1618  newExchangeStdVec[chainId] = currExchangeStdVec[chainId];
1619  }
1620 
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;
1625  break; // exit 'for'
1626  }
1627  }
1628 
1630  // Compute new ratio of positions per node
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; // Yes, 'final'
1636  newNumChainsPerNode [nodeId] += 1;
1637  newNumPositionsPerNode[nodeId] += newExchangeStdVec[chainId].numberOfPositions;
1638  }
1639 
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;
1648  }
1649  if (newNumPositionsPerNode[nodeId] < newSmallestAmountOfPositionsPerNode) {
1650  newSmallestAmountOfPositionsPerNode = newNumPositionsPerNode[nodeId];
1651  newNodeWithLeastPositions = nodeId;
1652  }
1653  }
1654 
1655  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1656  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::justBalance_proc0()"
1657  << ", level " << m_currLevel+LEVEL_REF_ID
1658  << ", step " << m_currStep
1659  << ", iter " << iterId
1660  << ", after update"
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] << ")"
1665  << std::endl;
1666  }
1667 
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);
1671 
1672  queso_require_equal_to_msg(newMinPosPerNode, newNumPositionsPerNode[newNodeWithLeastPositions], "inconsistent newMinPosPerNode");
1673 
1674  queso_require_equal_to_msg(newMaxPosPerNode, newNumPositionsPerNode[newNodeWithMostPositions], "inconsistent newMaxPosPerNode");
1675 
1676  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1677  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::justBalance_proc0()"
1678  << ", level " << m_currLevel+LEVEL_REF_ID
1679  << ", step " << m_currStep
1680  << ", iter " << iterId
1681  << ", newMaxPosPerNode = " << newMaxPosPerNode
1682  << ", newMinPosPerNode = " << newMinPosPerNode
1683  << ", newRatioOfPosPerNode = " << newRatioOfPosPerNode
1684  << std::endl;
1685  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1686  *m_env.subDisplayFile() << " KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1687  << ", level " << m_currLevel+LEVEL_REF_ID
1688  << ", step " << m_currStep
1689  << ", iter " << iterId
1690  << ", newNumChainsPerNode[" << nodeId << "] = " << newNumChainsPerNode [nodeId]
1691  << ", newNumPositionsPerNode[" << nodeId << "] = " << newNumPositionsPerNode[nodeId]
1692  << std::endl;
1693  }
1694  }
1695 
1697  // See if we need to exit 'while'
1699  if (newRatioOfPosPerNode > currRatioOfPosPerNode) {
1700  break; // exit 'while'
1701  }
1702 
1704  // Prepare for next iteration
1706  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1707  currNumChainsPerNode [nodeId] = 0;
1708  currNumPositionsPerNode[nodeId] = 0;
1709  }
1710  currRatioOfPosPerNode = newRatioOfPosPerNode;
1711  for (unsigned int chainId = 0; chainId < Nc; ++chainId) {
1712  currExchangeStdVec[chainId] = newExchangeStdVec[chainId];
1713  unsigned int nodeId = currExchangeStdVec[chainId].finalNodeOfInitialPosition; // Yes, 'final'
1714  currNumChainsPerNode [nodeId] += 1;
1715  currNumPositionsPerNode[nodeId] += currExchangeStdVec[chainId].numberOfPositions;
1716  }
1717  currMinPosPerNode = *std::min_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1718  currMaxPosPerNode = *std::max_element(currNumPositionsPerNode.begin(), currNumPositionsPerNode.end());
1719  currRatioOfPosPerNode = ((double) currMaxPosPerNode ) / ((double) currMinPosPerNode);
1720  }
1721 
1723  // Prepare output information
1725  for (unsigned int chainId = 0; chainId < Nc; ++chainId) {
1726  exchangeStdVec[chainId].finalNodeOfInitialPosition = currExchangeStdVec[chainId].finalNodeOfInitialPosition; // Yes, 'final' = 'final'
1727  }
1728 
1730  // Printout solution information
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; // Yes, 'final'
1736  finalNumChainsPerNode [nodeId] += 1;
1737  finalNumPositionsPerNode[nodeId] += exchangeStdVec[chainId].numberOfPositions;
1738  }
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);
1742 
1743  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1744  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1745  << ", level " << m_currLevel+LEVEL_REF_ID
1746  << ", step " << m_currStep
1747  << ": solution gives the following redistribution"
1748  << std::endl;
1749  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1750  *m_env.subDisplayFile() << " KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1751  << ", level " << m_currLevel+LEVEL_REF_ID
1752  << ", step " << m_currStep
1753  << ", finalNumChainsPerNode[" << nodeId << "] = " << finalNumChainsPerNode[nodeId]
1754  << ", finalNumPositionsPerNode[" << nodeId << "] = " << finalNumPositionsPerNode[nodeId]
1755  << std::endl;
1756  }
1757  *m_env.subDisplayFile() << " KEY In MLSampling<P_V,P_M>::justBalance_proc0()"
1758  << ", level " << m_currLevel+LEVEL_REF_ID
1759  << ", step " << m_currStep
1760  << ", finalRatioOfPosPerNode = " << finalRatioOfPosPerNode
1761  << std::endl;
1762  }
1763 
1765  // Measure time
1767  double balRunTime = MiscGetEllapsedSeconds(&timevalBal);
1768  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1769  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::justBalance_proc0()"
1770  << ", level " << m_currLevel+LEVEL_REF_ID
1771  << ", step " << m_currStep
1772  << ", iterId = " << iterId
1773  << ", currRatioOfPosPerNode = " << currRatioOfPosPerNode
1774  << ", after " << balRunTime << " seconds"
1775  << std::endl;
1776  }
1777 
1778  return;
1779 }
1780 
1781 template <class P_V,class P_M>
1782 void
1783 MLSampling<P_V,P_M>::mpiExchangePositions_inter0( // EXTRA FOR LOAD BALANCE
1784  const SequenceOfVectors<P_V,P_M>& prevChain, // input
1785  double prevExponent, // input
1786  double currExponent, // input
1787  const ScalarSequence<double>& prevLogLikelihoodValues, // input
1788  const ScalarSequence<double>& prevLogTargetValues, // input
1789  const std::vector<ExchangeInfoStruct>& exchangeStdVec, // input
1790  const std::vector<unsigned int>& finalNumChainsPerNode, // input
1791  const std::vector<unsigned int>& finalNumPositionsPerNode, // input
1792  BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl) // output
1793 {
1794  if (m_env.inter0Rank() < 0) {
1795  return;
1796  }
1797 
1798  unsigned int Np = (unsigned int) m_env.inter0Comm().NumProc();
1799  unsigned int Nc = exchangeStdVec.size();
1800 
1801  double expRatio = currExponent;
1802  if (prevExponent > 0.0) {
1803  expRatio /= prevExponent;
1804  }
1805 
1806  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1807  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1808  << ", level " << m_currLevel+LEVEL_REF_ID
1809  << ", step " << m_currStep
1810  << ": Np = " << Np
1811  << ", Nc = " << Nc
1812  << std::endl;
1813  }
1814 
1816  // Each node performs:
1817  // --> a 'gatherv' for collecting all necessary initial positions from other nodes
1818  // --> a 'gatherv' for collecting all necessary chain lenghts from other nodes
1820  for (unsigned int r = 0; r < Np; ++r) {
1822  // Prepare some counters
1824  unsigned int numberOfInitialPositionsNodeRAlreadyHas = 0;
1825  std::vector<unsigned int> numberOfInitialPositionsNodeRHasToReceiveFromNode(Np,0);
1826  std::vector<unsigned int> indexesOfInitialPositionsNodeRHasToReceiveFromMe(0);
1827 
1828  unsigned int sumOfChainLenghtsNodeRAlreadyHas = 0;
1829  std::vector<unsigned int> chainLenghtsNodeRHasToInherit(0);
1830 
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;
1836  }
1837  else {
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);
1842  }
1843  }
1844  }
1845  }
1846 
1847  unsigned int totalNumberOfInitialPositionsNodeRHasToReceive = 0;
1848  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) {
1849  totalNumberOfInitialPositionsNodeRHasToReceive += numberOfInitialPositionsNodeRHasToReceiveFromNode[nodeId];
1850  }
1851 
1852  unsigned int totalNumberOfChainLenghtsNodeRHasToInherit = chainLenghtsNodeRHasToInherit.size();
1853  unsigned int totalSumOfChainLenghtsNodeRHasToInherit = 0;
1854  for (unsigned int i = 0; i < totalNumberOfChainLenghtsNodeRHasToInherit; ++i) {
1855  totalSumOfChainLenghtsNodeRHasToInherit += chainLenghtsNodeRHasToInherit[i];
1856  }
1857 
1859  // Printout important information
1861  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1862  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1863  << ", level " << m_currLevel+LEVEL_REF_ID
1864  << ", step " << m_currStep
1865  << ": r = " << r
1866  << ", finalNumChainsPerNode[r] = " << finalNumChainsPerNode[r]
1867  << ", totalNumberOfInitialPositionsNodeRHasToReceive = " << totalNumberOfInitialPositionsNodeRHasToReceive
1868  << ", numberOfInitialPositionsNodeRAlreadyHas = " << numberOfInitialPositionsNodeRAlreadyHas
1869  << std::endl;
1870  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1871  << ", level " << m_currLevel+LEVEL_REF_ID
1872  << ", step " << m_currStep
1873  << ": r = " << r
1874  << ", finalNumPositionsPerNode[r] = " << finalNumPositionsPerNode[r]
1875  << ", totalSumOfChainLenghtsNodeRHasToInherit = " << totalSumOfChainLenghtsNodeRHasToInherit
1876  << ", sumOfChainLenghtsNodeRAlreadyHas = " << sumOfChainLenghtsNodeRAlreadyHas
1877  << std::endl;
1878  }
1879 
1881  // Make sanity checks
1883  queso_require_equal_to_msg(indexesOfInitialPositionsNodeRHasToReceiveFromMe.size(), numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()], "inconsistent number of initial positions to send to node 'r'");
1884 
1885  queso_require_equal_to_msg(finalNumChainsPerNode[r], (totalNumberOfInitialPositionsNodeRHasToReceive + numberOfInitialPositionsNodeRAlreadyHas), "inconsistent number of chains in node 'r'");
1886 
1887  queso_require_equal_to_msg(finalNumPositionsPerNode[r], (totalSumOfChainLenghtsNodeRHasToInherit + sumOfChainLenghtsNodeRAlreadyHas), "inconsistent sum of chain lenghts in node 'r'");
1888 
1889  queso_require_equal_to_msg(totalNumberOfInitialPositionsNodeRHasToReceive, totalNumberOfChainLenghtsNodeRHasToInherit, "inconsistent on total number of initial positions to receive in node 'r'");
1890 
1891  // Optimize use of memory (FIX ME: don't need to use swap here ????)
1892  indexesOfInitialPositionsNodeRHasToReceiveFromMe.resize(numberOfInitialPositionsNodeRHasToReceiveFromNode[m_env.inter0Rank()]);
1893  chainLenghtsNodeRHasToInherit.resize (totalSumOfChainLenghtsNodeRHasToInherit);
1894 
1896  // Prepare counters and buffers for gatherv of initial positions
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];
1908  prevChain.getPositionValues(auxIndex,auxInitialPosition);
1909  for (unsigned int j = 0; j < dimSize; ++j) {
1910  sendbuf[i*nValuesPerInitialPosition + j] = auxInitialPosition[j];
1911  }
1912  sendbuf[i*nValuesPerInitialPosition + dimSize] = prevLogLikelihoodValues[auxIndex];
1913  sendbuf[i*nValuesPerInitialPosition + dimSize + 1] = prevLogTargetValues[auxIndex];
1914  }
1915  }
1916 
1917  std::vector<double> recvbuf(0);
1918  std::vector<int> recvcnts(Np,0); // '0' is already the correct value for recvcnts[r]
1919  if (m_env.inter0Rank() == (int) r) {
1920  recvbuf.resize(totalNumberOfInitialPositionsNodeRHasToReceive * nValuesPerInitialPosition);
1921  for (unsigned int nodeId = 0; nodeId < Np; ++nodeId) { // Yes, from '0' on (for 'r', numberOf...ToReceiveFromNode[r] = 0 anyway)
1922  recvcnts[nodeId] = numberOfInitialPositionsNodeRHasToReceiveFromNode[nodeId/*m_env.inter0Rank()*/] * nValuesPerInitialPosition;
1923  }
1924  }
1925 
1926  std::vector<int> displs(Np,0);
1927  for (unsigned int nodeId = 1; nodeId < Np; ++nodeId) { // Yes, from '1' on
1928  displs[nodeId] = displs[nodeId-1] + recvcnts[nodeId-1];
1929  }
1930 
1931 #if 0
1932  if (m_env.inter0Rank() == r) {
1933  m_env.inter0Comm().Gatherv(RawValue_MPI_IN_PLACE, (int) sendcnt, RawValue_MPI_DOUBLE, (void *) &recvbuf[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, r, // LOAD BALANCE
1934  "MLSampling<P_V,P_M>::mpiExchangePositions_inter0(1)",
1935  "failed MPI.Gatherv()");
1936  }
1937  else {
1938  m_env.inter0Comm().Gatherv((void *) &sendbuf[0], (int) sendcnt, RawValue_MPI_DOUBLE, (void *) &recvbuf[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, r, // LOAD BALANCE
1939  "MLSampling<P_V,P_M>::mpiExchangePositions_inter0(2)",
1940  "failed MPI.Gatherv()");
1941  }
1942 #else
1943  m_env.inter0Comm().template Gatherv<double>(&sendbuf[0], (int) sendcnt,
1944  &recvbuf[0], (int *) &recvcnts[0], (int *) &displs[0],
1945  r, // LOAD BALANCE
1946  "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
1947  "failed MPI.Gatherv()");
1948 #endif
1949 
1951  // Make sanity checks
1953 
1955  // Transfer data from 'recvbuf' to 'balancedLinkControl'
1956  // Remember that finalNumChainsPerNode[r] = (totalNumberOfInitialPositionsNodeRHasToReceive + numberOfInitialPositionsNodeRAlreadyHas)
1957  // Remember that totalNumberOfInitialPositionsNodeRHasToReceive = totalNumberOfChainLenghtsNodeRHasToInherit
1959  if (m_env.inter0Rank() == (int) r) {
1960  balancedLinkControl.balLinkedChains.resize(finalNumChainsPerNode[r]);
1961  unsigned int auxIndex = 0;
1962 
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;
1967  prevChain.getPositionValues(originalIndex, auxInitialPosition);
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;
1972  auxIndex++;
1973  }
1974  }
1975 
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];
1979  }
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]; // aqui 3
1986  auxIndex++;
1987  }
1988  }
1989 
1990  m_env.inter0Comm().Barrier();
1991  } // for 'r'
1992 
1993  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1994  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1995  << ", level " << m_currLevel+LEVEL_REF_ID
1996  << ", step " << m_currStep
1997  << std::endl;
1998  }
1999 
2000  return;
2001 }
2002 
2003 // Statistical/private methods-----------------------
2004 template <class P_V,class P_M>
2005 void
2006 MLSampling<P_V,P_M>::checkpointML(
2007  double currExponent, // input
2008  double currEta, // input
2009  const SequenceOfVectors<P_V,P_M>& currChain, // input
2010  const ScalarSequence<double>& currLogLikelihoodValues, // input
2011  const ScalarSequence<double>& currLogTargetValues) // input
2012 {
2013  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2014  *m_env.subDisplayFile() << "\n CHECKPOINTING initiating at level " << m_currLevel
2015  << "\n" << std::endl;
2016  }
2017 
2018  //******************************************************************************
2019  // Write 'control' file without 'level' spefication in name
2020  //******************************************************************************
2021  unsigned int quantity1 = currChain.unifiedSequenceSize();
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) {
2025  queso_require_equal_to_msg(m_logEvidenceFactors.size(), m_currLevel, "number of evidence factors is not consistent");
2026  queso_require_equal_to_msg(quantity1, quantity2, "quantity2 is not consistent");
2027  queso_require_equal_to_msg(quantity1, quantity3, "quantity3 is not consistent");
2028  }
2029 
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 // 1
2034  << m_vectorSpace.dimGlobal() << std::endl // 2
2035  << currExponent << std::endl // 3
2036  << currEta << std::endl // 4
2037  << quantity1 << std::endl; // 5
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;
2042  }
2043  ofsVar->precision(savedPrecision);
2044  *ofsVar << "COMPLETE" << std::endl; // 6 = ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
2045 
2046  delete ofsVar;
2047  }
2048  m_env.fullComm().Barrier();
2049 
2050  //******************************************************************************
2051  // Write three 'data' files
2052  //******************************************************************************
2053  char levelSufix[256];
2054  sprintf(levelSufix,"%d",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
2055 
2056  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2057  *m_env.subDisplayFile() << "\n CHECKPOINTING chain at level " << m_currLevel
2058  << "\n" << std::endl;
2059  }
2060  currChain.unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles + "Chain_l" + levelSufix,
2061  m_options.m_restartOutput_fileType);
2062  m_env.fullComm().Barrier();
2063 
2064  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2065  *m_env.subDisplayFile() << "\n CHECKPOINTING like at level " << m_currLevel
2066  << "\n" << std::endl;
2067  }
2068  currLogLikelihoodValues.unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles + "LogLike_l" + levelSufix,
2069  m_options.m_restartOutput_fileType);
2070  m_env.fullComm().Barrier();
2071 
2072  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2073  *m_env.subDisplayFile() << "\n CHECKPOINTING target at level " << m_currLevel
2074  << "\n" << std::endl;
2075  }
2076  currLogTargetValues.unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles + "LogTarget_l" + levelSufix,
2077  m_options.m_restartOutput_fileType);
2078  m_env.fullComm().Barrier();
2079 
2080  //******************************************************************************
2081  // Write 'control' file *with* 'level' spefication in name
2082  //******************************************************************************
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 // 1
2087  << m_vectorSpace.dimGlobal() << std::endl // 2
2088  << currExponent << std::endl // 3
2089  << currEta << std::endl // 4
2090  << quantity1 << std::endl; // 5
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;
2095  }
2096  ofsVar->precision(savedPrecision);
2097  *ofsVar << "COMPLETE" << std::endl; // 6 = ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
2098 
2099  delete ofsVar;
2100  }
2101  m_env.fullComm().Barrier();
2102 
2103  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2104  *m_env.subDisplayFile() << "\n CHECKPOINTING done at level " << m_currLevel
2105  << "\n" << std::endl;
2106  }
2107 
2108  return;
2109 }
2110 //---------------------------------------------------
2111 template <class P_V,class P_M>
2112 void
2113 MLSampling<P_V,P_M>::restartML(
2114  double& currExponent, // output
2115  double& currEta, // output
2116  SequenceOfVectors<P_V,P_M>& currChain, // output
2117  ScalarSequence<double>& currLogLikelihoodValues, // output
2118  ScalarSequence<double>& currLogTargetValues) // output
2119 {
2120  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2121  *m_env.subDisplayFile() << "\n RESTARTING initiating at level " << m_currLevel
2122  << "\n" << std::endl;
2123  }
2124 
2125  //******************************************************************************
2126  // Read 'control' file
2127  //******************************************************************************
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(),
2133  std::ifstream::in);
2134 
2135  //******************************************************************************
2136  // Determine number of lines
2137  //******************************************************************************
2138  unsigned int numLines = std::count(std::istreambuf_iterator<char>(*ifsVar),
2139  std::istreambuf_iterator<char>(),
2140  '\n');
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
2144  << " lines"
2145  << std::endl;
2146  }
2147 
2148  //******************************************************************************
2149  // Read all values
2150  //******************************************************************************
2151  *ifsVar >> m_currLevel; // 1
2152  queso_require_equal_to_msg(numLines, (ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA + m_currLevel), "number of lines read is different than pre-established number of lines in control file");
2153 
2154  m_logEvidenceFactors.clear();
2155  m_logEvidenceFactors.resize(m_currLevel,0.);
2156  *ifsVar >> vectorSpaceDim // 2
2157  >> currExponent // 3
2158  >> currEta // 4
2159  >> quantity1; // 5
2160  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2161  *ifsVar >> m_logEvidenceFactors[i];
2162  }
2163  *ifsVar >> checkingString; // 6 = ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
2164  queso_require_equal_to_msg(checkingString, "COMPLETE", "control txt input file is not complete");
2165 
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];
2175  }
2176  *m_env.subDisplayFile() << std::endl;
2177  }
2178 
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;
2184  }
2185  while ((lineId < numLines) && (ifsVar->eof() == false)) {
2186  }
2187  ifsVar->ignore(maxCharsPerLine,'\n');
2188 #endif
2189 
2190  delete ifsVar;
2191  } // if (m_env.fullRank() == 0)
2192  m_env.fullComm().Barrier();
2193 
2194  //******************************************************************************
2195  // MPI_Bcast 'm_currLevel'
2196  //******************************************************************************
2197  unsigned int tmpUint = (unsigned int) m_currLevel;
2198  m_env.fullComm().Bcast((void *) &tmpUint, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'fullComm'
2199  "MLSampling<P_V,P_M>::restartML()",
2200  "failed MPI.Bcast() for m_currLevel");
2201  if (m_env.fullRank() != 0) {
2202  m_currLevel = tmpUint;
2203  }
2204 
2205  //******************************************************************************
2206  // MPI_Bcast the rest of the information just read
2207  //******************************************************************************
2208  std::vector<double> tmpData(ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA-1+m_currLevel,0.);
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];
2216  }
2217  }
2218  else {
2219  m_logEvidenceFactors.clear();
2220  m_logEvidenceFactors.resize(m_currLevel,0.);
2221  }
2222  m_env.fullComm().Bcast((void *) &tmpData[0], (int) tmpData.size(), RawValue_MPI_DOUBLE, 0, // Yes, 'fullComm'
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];
2232  }
2233  }
2234 
2235  //******************************************************************************
2236  // Process read data in all MPI nodes now
2237  //******************************************************************************
2238  queso_require_equal_to_msg(vectorSpaceDim, m_vectorSpace.dimGlobal(), "read vector space dimension is not consistent");
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());
2243 
2244  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2245  *m_env.subDisplayFile() << "Restart input file has the following information"
2246  << ": subSequenceSize = " << subSequenceSize
2247  << std::endl;
2248  }
2249 
2250  //******************************************************************************
2251  // Read three 'data' files
2252  //******************************************************************************
2253  char levelSufix[256];
2254  sprintf(levelSufix,"%d",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
2255 
2256  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2257  *m_env.subDisplayFile() << "\n RESTARTING chain at level " << m_currLevel
2258  << "\n" << std::endl;
2259  }
2260  currChain.unifiedReadContents(m_options.m_restartInput_baseNameForFiles + "Chain_l" + levelSufix,
2261  m_options.m_restartInput_fileType,
2262  subSequenceSize);
2263  m_env.fullComm().Barrier();
2264 
2265  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2266  *m_env.subDisplayFile() << "\n RESTARTING like at level " << m_currLevel
2267  << "\n" << std::endl;
2268  }
2269  currLogLikelihoodValues.unifiedReadContents(m_options.m_restartInput_baseNameForFiles + "LogLike_l" + levelSufix,
2270  m_options.m_restartInput_fileType,
2271  subSequenceSize);
2272  m_env.fullComm().Barrier();
2273 
2274  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2275  *m_env.subDisplayFile() << "\n RESTARTING target at level " << m_currLevel
2276  << "\n" << std::endl;
2277  }
2278  currLogTargetValues.unifiedReadContents(m_options.m_restartInput_baseNameForFiles + "LogTarget_l" + levelSufix,
2279  m_options.m_restartInput_fileType,
2280  subSequenceSize);
2281  m_env.fullComm().Barrier();
2282 
2283  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2284  *m_env.subDisplayFile() << "\n RESTARTING done at level " << m_currLevel
2285  << "\n" << std::endl;
2286  }
2287 
2288  return;
2289 }
2290 //---------------------------------------------------
2291 template <class P_V,class P_M>
2292 void
2293 MLSampling<P_V,P_M>::generateSequence_Level0_all(
2294  const MLSamplingLevelOptions& currOptions, // input
2295  unsigned int& unifiedRequestedNumSamples, // output
2296  SequenceOfVectors<P_V,P_M>& currChain, // output
2297  ScalarSequence<double>& currLogLikelihoodValues, // output
2298  ScalarSequence<double>& currLogTargetValues) // output
2299 {
2300  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2301  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateSequence()"
2302  << ": beginning level " << m_currLevel+LEVEL_REF_ID
2303  << ", currOptions.m_rawChainSize = " << currOptions.m_rawChainSize // Ok to use rawChainSize
2304  << std::endl;
2305  }
2306 
2307  int iRC = UQ_OK_RC;
2308  struct timeval timevalLevel;
2309  iRC = gettimeofday(&timevalLevel, NULL);
2310  if (iRC) {}; // just to remove compiler warning
2311 
2312  if (m_env.inter0Rank() >= 0) {
2313  unsigned int tmpSize = currOptions.m_rawChainSize;
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");
2317  }
2318  else {
2319  unifiedRequestedNumSamples = currOptions.m_rawChainSize;
2320  }
2321 
2322  currChain.setName (currOptions.m_prefix + "rawChain" );
2323  currLogLikelihoodValues.setName(currOptions.m_prefix + "rawLogLikelihood");
2324  currLogTargetValues.setName (currOptions.m_prefix + "rawLogTarget" );
2325 
2326  currChain.resizeSequence (currOptions.m_rawChainSize); // Ok to use rawChainSize
2327  currLogLikelihoodValues.resizeSequence(currOptions.m_rawChainSize); // Ok to use rawChainSize
2328  currLogTargetValues.resizeSequence (currOptions.m_rawChainSize); // Ok to use rawChainSize
2329 
2330  P_V auxVec(m_vectorSpace.zeroVector());
2331  ScalarFunctionSynchronizer<P_V,P_M> likelihoodSynchronizer(m_likelihoodFunction,auxVec); // prudencio 2010-08-01
2332  for (unsigned int i = 0; i < currChain.subSequenceSize(); ++i) {
2333  //std::cout << "In QUESO: before prior realizer with i = " << i << std::endl;
2334  bool outOfSupport = true;
2335  do {
2336 
2337  m_priorRv.realizer().realization(auxVec); // gpmsa2
2338  if (m_numDisabledParameters > 0) { // gpmsa2
2339  unsigned int disabledCounter = 0;
2340  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2341  if (m_parameterEnabledStatus[paramId] == false) {
2342  auxVec[paramId] = currOptions.m_initialValuesOfDisabledParameters[disabledCounter];
2343  disabledCounter++;
2344  }
2345  }
2346  }
2347  auxVec.mpiBcast(0, m_env.subComm()); // prudencio 2010-08-01
2348 
2349  outOfSupport = !(m_targetDomain->contains(auxVec));
2350  } while (outOfSupport); // prudenci 2011-Oct-04
2351 
2352  currChain.setPositionValues(i,auxVec);
2353  // KAUST: all nodes should call likelihood
2354 #if 1 // prudencio 2010-08-01
2355  currLogLikelihoodValues[i] = likelihoodSynchronizer.callFunction(&auxVec,NULL,NULL,NULL,NULL,NULL,NULL); // likelihood is important
2356 #else
2357  currLogLikelihoodValues[i] = m_likelihoodFunction.lnValue(auxVec,NULL,NULL,NULL,NULL); // likelihood is important
2358 #endif
2359  currLogTargetValues[i] = m_priorRv.pdf().lnValue(auxVec,NULL,NULL,NULL,NULL) + currLogLikelihoodValues[i];
2360  //std::cout << "In QUESO: currLogTargetValues[" << i << "] = " << currLogTargetValues[i] << std::endl;
2361  }
2362 
2363  if (m_env.inter0Rank() >= 0) { // KAUST
2364 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2365  if (currOptions.m_rawChainComputeStats) {
2366  FilePtrSetStruct filePtrSet;
2367  m_env.openOutputFile(currOptions.m_dataOutputFileName,
2368  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
2369  currOptions.m_dataOutputAllowedSet,
2370  false,
2371  filePtrSet);
2372 
2373  //m_env.syncPrintDebugMsg("At level 0, calling computeStatistics for chain",1,10,m_env.inter0Comm()); // output debug
2374  currChain.computeStatistics(*currOptions.m_rawChainStatisticalOptionsObj,
2375  filePtrSet.ofsVar);
2376 
2377  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
2378  }
2379  // Compute MLE and MAP
2380  // rr0
2381 #endif
2382  if (currOptions.m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
2383  currChain.unifiedWriteContents (currOptions.m_rawChainDataOutputFileName,
2384  currOptions.m_rawChainDataOutputFileType);
2385  currLogLikelihoodValues.unifiedWriteContents(currOptions.m_rawChainDataOutputFileName,
2386  currOptions.m_rawChainDataOutputFileType);
2387  currLogTargetValues.unifiedWriteContents (currOptions.m_rawChainDataOutputFileName,
2388  currOptions.m_rawChainDataOutputFileType);
2389  }
2390 
2391  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2392  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2393  << ", level " << m_currLevel+LEVEL_REF_ID
2394  << ": finished generating " << currChain.subSequenceSize()
2395  << " chain positions"
2396  << std::endl;
2397 
2398  //unsigned int numZeros = 0;
2399  //for (unsigned int i = 0; i < currTargetValues.subSequenceSize(); ++i) {
2400  // *m_env.subDisplayFile() << "currTargetValues[" << i
2401  // << "] = " << currTargetValues[i]
2402  // << std::endl;
2403  // if (currTargetValues[i] == 0.) numZeros++;
2404  //}
2405  //*m_env.subDisplayFile() << "Number of zeros in currTargetValues = " << numZeros
2406  // << std::endl;
2407  }
2408 
2409  if (currOptions.m_filteredChainGenerate) {
2410  // todo
2411 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2412  if (currOptions.m_filteredChainComputeStats) {
2413  // todo
2414 
2415  //currChain.computeStatistics(*currOptions.m_filteredChainStatisticalOptionsObj,
2416  // filePtrSet.ofsVar);
2417  }
2418  // Compute MLE and MAP
2419  // rr0
2420 #endif
2421  }
2422 
2423  } // KAUST
2424 
2425  queso_require_equal_to_msg(currChain.subSequenceSize(), currOptions.m_rawChainSize, "currChain (first one) has been generated with invalid size");
2426 
2427  double levelRunTime = MiscGetEllapsedSeconds(&timevalLevel);
2428  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2429  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2430  << ": ending level " << m_currLevel+LEVEL_REF_ID
2431  << ", total level time = " << levelRunTime << " seconds"
2432  << std::endl;
2433  }
2434 
2435  return;
2436 }
2437 //---------------------------------------------------
2438 template <class P_V,class P_M>
2439 void
2440 MLSampling<P_V,P_M>::generateSequence_Step01_inter0(
2441  const MLSamplingLevelOptions* currOptions, // input
2442  unsigned int& unifiedRequestedNumSamples) // output
2443 {
2444  int iRC = UQ_OK_RC;
2445  struct timeval timevalStep;
2446  iRC = gettimeofday(&timevalStep, NULL);
2447  if (iRC) {}; // just to remove compiler warning
2448 
2449  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2450  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2451  << ", level " << m_currLevel+LEVEL_REF_ID
2452  << ", step " << m_currStep
2453  << ": beginning step 1 of 11"
2454  << std::endl;
2455  }
2456 
2457  unsigned int tmpSize = currOptions->m_rawChainSize;
2458  // This computed 'unifiedRequestedNumSamples' needs to be recomputed only at the last
2459  // level, when 'currOptions' is replaced by 'lastLevelOptions' (see step 3 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");
2463 
2464  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2465  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateSequence()"
2466  << ", level " << m_currLevel+LEVEL_REF_ID
2467  << ", step " << m_currStep
2468  << ", currOptions->m_rawChainSize = " << currOptions->m_rawChainSize // Ok to use rawChainSize
2469  << std::endl;
2470  }
2471 
2472  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
2473  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2474  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2475  << ", level " << m_currLevel+LEVEL_REF_ID
2476  << ", step " << m_currStep
2477  << ", after " << stepRunTime << " seconds"
2478  << std::endl;
2479  }
2480 
2481  return;
2482 }
2483 //---------------------------------------------------
2484 template <class P_V,class P_M>
2485 void
2486 MLSampling<P_V,P_M>::generateSequence_Step02_inter0(
2487  const MLSamplingLevelOptions* currOptions, // input
2488  SequenceOfVectors<P_V,P_M>& currChain, // input/output
2489  ScalarSequence<double>& currLogLikelihoodValues, // input/output
2490  ScalarSequence<double>& currLogTargetValues, // input/output
2491  SequenceOfVectors<P_V,P_M>& prevChain, // output
2492  ScalarSequence<double>& prevLogLikelihoodValues, // output
2493  ScalarSequence<double>& prevLogTargetValues, // output
2494  unsigned int& indexOfFirstWeight, // output
2495  unsigned int& indexOfLastWeight) // output
2496 {
2497  int iRC = UQ_OK_RC;
2498  struct timeval timevalStep;
2499  iRC = gettimeofday(&timevalStep, NULL);
2500  if (iRC) {}; // just to remove compiler warning
2501 
2502  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2503  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2504  << ", level " << m_currLevel+LEVEL_REF_ID
2505  << ", step " << m_currStep
2506  << ": beginning step 2 of 11"
2507  << std::endl;
2508  }
2509 
2510  prevChain = currChain;
2511  currChain.clear();
2512  currChain.setName(currOptions->m_prefix + "rawChain");
2513 
2514  prevLogLikelihoodValues = currLogLikelihoodValues; // likelihood is important
2515  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2516  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
2517  << ", level " << m_currLevel+LEVEL_REF_ID
2518  << ", step " << m_currStep
2519  << ", prevLogLikelihoodValues[0] = " << prevLogLikelihoodValues[0]
2520  << std::endl;
2521  }
2522  prevLogTargetValues = currLogTargetValues;
2523 
2524  currLogLikelihoodValues.clear();
2525  currLogLikelihoodValues.setName(currOptions->m_prefix + "rawLogLikelihood");
2526 
2527  currLogTargetValues.clear();
2528  currLogTargetValues.setName(currOptions->m_prefix + "rawLogTarget");
2529 
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()"
2534  << ", level " << m_currLevel+LEVEL_REF_ID
2535  << ", step " << m_currStep
2536  << ":"
2537  << std::endl;
2538  for (unsigned int i = 0; i < prevChain.subSequenceSize(); ++i) {
2539  prevChain.getPositionValues(i,prevPosition);
2540  *m_env.subDisplayFile() << " prevChain[" << i
2541  << "] = " << prevPosition
2542  << ", prevLogLikelihoodValues[" << i
2543  << "] = " << prevLogLikelihoodValues[i]
2544  << ", prevLogTargetValues[" << i
2545  << "] = " << prevLogTargetValues[i]
2546  << std::endl;
2547  }
2548  }
2549 #endif
2550 
2551  unsigned int quantity1 = prevChain.unifiedSequenceSize();
2552  unsigned int quantity2 = currChain.unifiedSequenceSize();
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()"
2559  << ", level " << m_currLevel+LEVEL_REF_ID
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
2567  << std::endl;
2568  }
2569 
2570  queso_require_equal_to_msg(prevChain.subSequenceSize(), prevLogLikelihoodValues.subSequenceSize(), "different sizes between previous chain and previous sequence of likelihood values");
2571 
2572  queso_require_equal_to_msg(prevChain.subSequenceSize(), prevLogTargetValues.subSequenceSize(), "different sizes between previous chain and previous sequence of target values");
2573 
2574  // Set 'indexOfFirstWeight' and 'indexOfLastWeight' // KAUST
2575  indexOfFirstWeight = 0;
2576  indexOfLastWeight = indexOfFirstWeight + prevChain.subSequenceSize()-1;
2577  {
2578  //std::cout << "m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc() << std::endl;
2579  int r = m_env.inter0Rank();
2580  //std::cout << "r = " << r << std::endl;
2581  m_env.inter0Comm().Barrier();
2582  unsigned int auxUint = 0;
2583  if (r > 0) {
2584  RawType_MPI_Status status;
2585  //std::cout << "Rank " << r << " is entering MPI_Recv()" << std::endl;
2586  m_env.inter0Comm().Recv((void*) &auxUint, 1, RawValue_MPI_UNSIGNED, r-1, r-1, &status,
2587  "MLSampling<P_V,P_M>::generateSequence()",
2588  "failed MPI.Recv()");
2589  //std::cout << "Rank " << r << " received auxUint = " << auxUint << std::endl;
2590  indexOfFirstWeight = auxUint;
2591  indexOfLastWeight = indexOfFirstWeight + prevChain.subSequenceSize()-1;
2592  }
2593  if (r < (m_env.inter0Comm().NumProc()-1)) {
2594  auxUint = indexOfLastWeight + 1;
2595  //std::cout << "Rank " << r << " is sending auxUint = " << auxUint << std::endl;
2596  m_env.inter0Comm().Send((void*) &auxUint, 1, RawValue_MPI_UNSIGNED, r+1, r,
2597  "MLSampling<P_V,P_M>::generateSequence()",
2598  "failed MPI.Send()");
2599  //std::cout << "Rank " << r << " sent auxUint = " << auxUint << std::endl;
2600  }
2601  m_env.inter0Comm().Barrier();
2602  }
2603 
2604  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
2605  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2606  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2607  << ", level " << m_currLevel+LEVEL_REF_ID
2608  << ", step " << m_currStep
2609  << ", after " << stepRunTime << " seconds"
2610  << std::endl;
2611  }
2612 
2613  return;
2614 }
2615 //---------------------------------------------------
2616 template <class P_V,class P_M>
2617 void
2618 MLSampling<P_V,P_M>::generateSequence_Step03_inter0(
2619  const MLSamplingLevelOptions* currOptions, // input
2620  const ScalarSequence<double>& prevLogLikelihoodValues, // input
2621  double prevExponent, // input
2622  double failedExponent, // input // gpmsa1
2623  double& currExponent, // output
2624  ScalarSequence<double>& weightSequence) // output
2625 {
2626  int iRC = UQ_OK_RC;
2627  struct timeval timevalStep;
2628  iRC = gettimeofday(&timevalStep, NULL);
2629  if (iRC) {}; // just to remove compiler warning
2630 
2631  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2632  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2633  << ", level " << m_currLevel+LEVEL_REF_ID
2634  << ", step " << m_currStep
2635  << ", failedExponent = " << failedExponent // gpmsa1
2636  << ": beginning step 3 of 11"
2637  << std::endl;
2638  }
2639 
2640  std::vector<double> exponents(2,0.);
2641  exponents[0] = prevExponent;
2642  exponents[1] = 1.;
2643 
2644  double nowExponent = 1.; // Try '1.' right away
2645  double nowEffectiveSizeRatio = 0.; // To be computed
2646 
2647  unsigned int nowAttempt = 0;
2648  bool testResult = false;
2649  double meanEffectiveSizeRatio = .5*(currOptions->m_minEffectiveSizeRatio + currOptions->m_maxEffectiveSizeRatio);
2650  ScalarSequence<double> omegaLnDiffSequence(m_env,prevLogLikelihoodValues.subSequenceSize(),"");
2651 
2652  double nowUnifiedEvidenceLnFactor = 0.;
2653  do {
2654  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2655  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2656  << ", level " << m_currLevel+LEVEL_REF_ID
2657  << ", step " << m_currStep
2658  << ", failedExponent = " << failedExponent // gpmsa1
2659  << ": entering loop for computing next exponent"
2660  << ", with nowAttempt = " << nowAttempt
2661  << std::endl;
2662  }
2663 
2664  if (failedExponent > 0.) { // gpmsa1
2665  nowExponent = .5*(prevExponent+failedExponent);
2666  }
2667  else {
2668  if (nowAttempt > 0) {
2669  if (nowEffectiveSizeRatio > meanEffectiveSizeRatio) {
2670  exponents[0] = nowExponent;
2671  }
2672  else {
2673  exponents[1] = nowExponent;
2674  }
2675  nowExponent = .5*(exponents[0] + exponents[1]);
2676  }
2677  }
2678  double auxExponent = nowExponent;
2679  if (prevExponent != 0.) {
2680  auxExponent /= prevExponent;
2681  auxExponent -= 1.;
2682  }
2683  double subWeightRatioSum = 0.;
2684  double unifiedWeightRatioSum = 0.;
2685 
2686  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2687  omegaLnDiffSequence[i] = prevLogLikelihoodValues[i]*auxExponent; // likelihood is important
2688  }
2689 
2690 #if 1 // prudenci-2012-07-06
2691  //double unifiedOmegaLnMin = omegaLnDiffSequence.unifiedMinPlain(m_vectorSpace.numOfProcsForStorage() == 1);
2692  double unifiedOmegaLnMax = omegaLnDiffSequence.unifiedMaxPlain(m_vectorSpace.numOfProcsForStorage() == 1);
2693 #else
2694  double unifiedOmegaLnMin = 0.;
2695  double unifiedOmegaLnMax = 0.;
2696  omegaLnDiffSequence.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1, // KAUST3
2697  0,
2698  omegaLnDiffSequence.subSequenceSize(),
2699  unifiedOmegaLnMin,
2700  unifiedOmegaLnMax);
2701 #endif
2702  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
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()"
2710  << ", level " << m_currLevel+LEVEL_REF_ID
2711  << ", step " << m_currStep
2712  << ", i = " << i
2713  << ", prevLogLikelihoodValues[i] = " << prevLogLikelihoodValues[i]
2714  << ", omegaLnDiffSequence[i] = " << omegaLnDiffSequence[i]
2715  << ", weightSequence[i] = " << weightSequence[i]
2716  //<< ", subWeightRatioSum = " << subWeightRatioSum
2717  << std::endl;
2718  }
2719  }
2720 #endif
2721  }
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");
2725 
2726  unsigned int auxQuantity = weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2727  nowUnifiedEvidenceLnFactor = log(unifiedWeightRatioSum) + unifiedOmegaLnMax - log(auxQuantity);
2728 
2729  double effectiveSampleSize = 0.;
2730  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2731  weightSequence[i] /= unifiedWeightRatioSum;
2732  effectiveSampleSize += weightSequence[i]*weightSequence[i];
2733  //if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2734  // *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2735  // << ", level " << m_currLevel+LEVEL_REF_ID
2736  // << ", step " << m_currStep
2737  // << ": i = " << i
2738  // << ", effectiveSampleSize = " << effectiveSampleSize
2739  // << std::endl;
2740  //}
2741  }
2742 
2743  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2744  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2745  << ", level " << m_currLevel+LEVEL_REF_ID
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
2758  << std::endl;
2759  }
2760 
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()"
2764  << ", level " << m_currLevel+LEVEL_REF_ID
2765  << ", step " << m_currStep
2766  << ":"
2767  << std::endl;
2768  }
2769  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2770  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2771  *m_env.subDisplayFile() << " weightSequence[" << i
2772  << "] = " << weightSequence[i]
2773  << std::endl;
2774  }
2775  }
2776 #endif
2777 
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");
2783 
2784  effectiveSampleSize = 1./effectiveSampleSize;
2785  nowEffectiveSizeRatio = effectiveSampleSize/((double) weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1));
2786  queso_require_less_equal_msg(nowEffectiveSizeRatio, (1.+1.e-8), "effective sample size ratio cannot be > 1");
2787 
2788  // m_env.worldRank(),
2789  // "MLSampling<P_V,P_M>::generateSequence()",
2790  // "effective sample size ratio cannot be < 1");
2791 
2792  if (failedExponent > 0.) { // gpmsa1
2793  testResult = true;
2794  }
2795  else {
2796  //bool aux1 = (nowEffectiveSizeRatio == meanEffectiveSizeRatio);
2797  bool aux2 = (nowExponent == 1. )
2798  &&
2799  (nowEffectiveSizeRatio > meanEffectiveSizeRatio);
2800  bool aux3 = (nowEffectiveSizeRatio >= currOptions->m_minEffectiveSizeRatio)
2801  &&
2802  (nowEffectiveSizeRatio <= currOptions->m_maxEffectiveSizeRatio);
2803  testResult = aux2 || aux3;
2804  }
2805 
2806  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2807  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2808  << ", level " << m_currLevel+LEVEL_REF_ID
2809  << ", step " << m_currStep
2810  << ": nowAttempt = " << nowAttempt
2811  << ", prevExponent = " << prevExponent
2812  << ", failedExponent = " << failedExponent // gpmsa1
2813  << ", exponents[0] = " << exponents[0]
2814  << ", nowExponent = " << nowExponent
2815  << ", exponents[1] = " << exponents[1]
2816  << ", effectiveSampleSize = " << effectiveSampleSize
2817  << ", weightSequenceSize = " << weightSequence.subSequenceSize()
2818  << ", minEffectiveSizeRatio = " << currOptions->m_minEffectiveSizeRatio
2819  << ", nowEffectiveSizeRatio = " << nowEffectiveSizeRatio
2820  << ", maxEffectiveSizeRatio = " << currOptions->m_maxEffectiveSizeRatio
2821  //<< ", aux2 = " << aux2
2822  //<< ", aux3 = " << aux3
2823  << ", testResult = " << testResult
2824  << std::endl;
2825  }
2826  nowAttempt++;
2827 
2828  // Make sure all nodes in 'inter0Comm' have the same value of 'nowExponent'
2829  if (MiscCheckForSameValueInAllNodes(nowExponent,
2830  0., // kept 'zero' on 2010/03/05
2831  m_env.inter0Comm(),
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()"
2835  << ", level " << m_currLevel+LEVEL_REF_ID
2836  << ", step " << m_currStep
2837  << ": nowAttempt = " << nowAttempt
2838  << ", MiscCheck for 'nowExponent' detected a problem"
2839  << std::endl;
2840  }
2841  }
2842 
2843  // Make sure all nodes in 'inter0Comm' have the same value of 'testResult'
2844  if (MiscCheckForSameValueInAllNodes(testResult,
2845  0., // kept 'zero' on 2010/03/05
2846  m_env.inter0Comm(),
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()"
2850  << ", level " << m_currLevel+LEVEL_REF_ID
2851  << ", step " << m_currStep
2852  << ": nowAttempt = " << nowAttempt
2853  << ", MiscCheck for 'testResult' detected a problem"
2854  << std::endl;
2855  }
2856  }
2857  } while (testResult == false);
2858  currExponent = nowExponent;
2859  if (failedExponent > 0.) { // gpmsa1
2860  m_logEvidenceFactors[m_logEvidenceFactors.size()-1] = nowUnifiedEvidenceLnFactor;
2861  }
2862  else {
2863  m_logEvidenceFactors.push_back(nowUnifiedEvidenceLnFactor); // restart
2864  }
2865 
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()"
2869  << ", level " << m_currLevel+LEVEL_REF_ID
2870  << ", step " << m_currStep
2871  << ": weightSequence.subSequenceSize() = " << weightSequence.subSequenceSize()
2872  << ", weightSequence.unifiedSequenceSize() = " << quantity1
2873  << ", failedExponent = " << failedExponent // gpmsa1
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])
2878  << std::endl;
2879 
2880  //unsigned int numZeros = 0;
2881  //for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2882  // *m_env.subDisplayFile() << "weightSequence[" << i
2883  // << "] = " << weightSequence[i]
2884  // << std::endl;
2885  // if (weightSequence[i] == 0.) numZeros++;
2886  //}
2887  //*m_env.subDisplayFile() << "Number of zeros in weightSequence = " << numZeros
2888  // << std::endl;
2889  }
2890 
2891  // Make sure all nodes in 'inter0Comm' have the same value of 'logEvidenceFactor'
2892  if (MiscCheckForSameValueInAllNodes(m_logEvidenceFactors[m_logEvidenceFactors.size()-1],
2893  3.0e-16, // changed from 'zero' on 2010/03/03
2894  m_env.inter0Comm(),
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()"
2898  << ", level " << m_currLevel+LEVEL_REF_ID
2899  << ", step " << m_currStep
2900  << ", failedExponent = " << failedExponent // gpmsa1
2901  << ": nowAttempt = " << nowAttempt
2902  << ", MiscCheck for 'logEvidenceFactor' detected a problem"
2903  << std::endl;
2904  }
2905  }
2906 
2907  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
2908  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2909  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2910  << ", level " << m_currLevel+LEVEL_REF_ID
2911  << ", step " << m_currStep
2912  << ", failedExponent = " << failedExponent // gpmsa
2913  << ", after " << stepRunTime << " seconds"
2914  << std::endl;
2915  }
2916 
2917  return;
2918 }
2919 //---------------------------------------------------
2920 template <class P_V,class P_M>
2921 void
2922 MLSampling<P_V,P_M>::generateSequence_Step04_inter0(
2923  const SequenceOfVectors<P_V,P_M>& prevChain, // input
2924  const ScalarSequence<double>& weightSequence, // input
2925  P_M& unifiedCovMatrix) // output
2926 {
2927  int iRC = UQ_OK_RC;
2928  struct timeval timevalStep;
2929  iRC = gettimeofday(&timevalStep, NULL);
2930  if (iRC) {}; // just to remove compiler warning
2931 
2932  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2933  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2934  << ", level " << m_currLevel+LEVEL_REF_ID
2935  << ", step " << m_currStep
2936  << ": beginning step 4 of 11"
2937  << std::endl;
2938  }
2939 
2940  P_V auxVec(m_vectorSpace.zeroVector());
2941  P_V subWeightedMeanVec(m_vectorSpace.zeroVector());
2942  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2943  prevChain.getPositionValues(i,auxVec);
2944  subWeightedMeanVec += weightSequence[i]*auxVec;
2945  }
2946 
2947  // Todd Oliver 2010-09-07: compute weighted mean over all processors
2948  P_V unifiedWeightedMeanVec(m_vectorSpace.zeroVector());
2949  if (m_env.inter0Rank() >= 0) {
2950  subWeightedMeanVec.mpiAllReduce(RawValue_MPI_SUM,m_env.inter0Comm(),unifiedWeightedMeanVec);
2951  }
2952  else {
2953  unifiedWeightedMeanVec = subWeightedMeanVec;
2954  }
2955 
2956  P_V diffVec(m_vectorSpace.zeroVector());
2957  P_M subCovMatrix(m_vectorSpace.zeroVector());
2958  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2959  prevChain.getPositionValues(i,auxVec);
2960  diffVec = auxVec - unifiedWeightedMeanVec;
2961  subCovMatrix += weightSequence[i]*matrixProduct(diffVec,diffVec);
2962  }
2963 
2964  for (unsigned int i = 0; i < unifiedCovMatrix.numRowsLocal(); ++i) { // KAUST5
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");
2972  }
2973  else {
2974  sumValue = localValue;
2975  }
2976  unifiedCovMatrix(i,j) = sumValue;
2977  }
2978  }
2979 
2980  if (m_numDisabledParameters > 0) { // gpmsa2
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.;
2985  }
2986  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2987  unifiedCovMatrix(paramId,j) = 0.;
2988  }
2989  unifiedCovMatrix(paramId,paramId) = 1.;
2990  }
2991  }
2992  }
2993 
2994  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2995  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2996  << ", level " << m_currLevel+LEVEL_REF_ID
2997  << ", step " << m_currStep
2998  << ": unifiedCovMatrix = " << unifiedCovMatrix
2999  << std::endl;
3000  }
3001 
3002  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3003  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3004  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3005  << ", level " << m_currLevel+LEVEL_REF_ID
3006  << ", step " << m_currStep
3007  << ", after " << stepRunTime << " seconds"
3008  << std::endl;
3009  }
3010 
3011  return;
3012 }
3013 //---------------------------------------------------
3014 template <class P_V,class P_M>
3015 void
3016 MLSampling<P_V,P_M>::generateSequence_Step05_inter0(
3017  unsigned int unifiedRequestedNumSamples, // input
3018  const ScalarSequence<double>& weightSequence, // input
3019  std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // output
3020  std::vector<double>& unifiedWeightStdVectorAtProc0Only) // output
3021 {
3022  int iRC = UQ_OK_RC;
3023  struct timeval timevalStep;
3024  iRC = gettimeofday(&timevalStep, NULL);
3025  if (iRC) {}; // just to remove compiler warning
3026 
3027  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3028  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3029  << ", level " << m_currLevel+LEVEL_REF_ID
3030  << ", step " << m_currStep
3031  << ": beginning step 5 of 11"
3032  << std::endl;
3033  }
3034 
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()"
3038  << ", level " << m_currLevel+LEVEL_REF_ID
3039  << ", step " << m_currStep
3040  << ", before weightSequence.getUnifiedContentsAtProc0Only()"
3041  << ":"
3042  << std::endl;
3043  }
3044  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
3045  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3046  *m_env.subDisplayFile() << ", weightSequence[" << i
3047  << "] = " << weightSequence[i]
3048  << std::endl;
3049  }
3050  }
3051 #endif
3052 
3053  weightSequence.getUnifiedContentsAtProc0Only(m_vectorSpace.numOfProcsForStorage() == 1,
3054  unifiedWeightStdVectorAtProc0Only);
3055 
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()"
3059  << ", level " << m_currLevel+LEVEL_REF_ID
3060  << ", step " << m_currStep
3061  << ", after weightSequence.getUnifiedContentsAtProc0Only()"
3062  << ":"
3063  << std::endl;
3064  }
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]
3069  << std::endl;
3070  }
3071  }
3072 #endif
3073  sampleIndexes_proc0(unifiedRequestedNumSamples, // input
3074  unifiedWeightStdVectorAtProc0Only, // input
3075  unifiedIndexCountersAtProc0Only); // output
3076 
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");
3080  }
3081 
3082  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3083  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3084  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3085  << ", level " << m_currLevel+LEVEL_REF_ID
3086  << ", step " << m_currStep
3087  << ", after " << stepRunTime << " seconds"
3088  << std::endl;
3089  }
3090 
3091  return;
3092 }
3093 //---------------------------------------------------
3094 template <class P_V,class P_M>
3095 void
3096 MLSampling<P_V,P_M>::generateSequence_Step06_all(
3097  const MLSamplingLevelOptions* currOptions, // input
3098  unsigned int indexOfFirstWeight, // input
3099  unsigned int indexOfLastWeight, // input
3100  const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // input
3101  bool& useBalancedChains, // output
3102  std::vector<ExchangeInfoStruct>& exchangeStdVec) // output
3103 {
3104  int iRC = UQ_OK_RC;
3105  struct timeval timevalStep;
3106  iRC = gettimeofday(&timevalStep, NULL);
3107  if (iRC) {}; // just to remove compiler warning
3108 
3109  useBalancedChains = decideOnBalancedChains_all(currOptions, // input
3110  indexOfFirstWeight, // input
3111  indexOfLastWeight, // input
3112  unifiedIndexCountersAtProc0Only, // input
3113  exchangeStdVec); // output
3114 
3115  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3116  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3117  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3118  << ", level " << m_currLevel+LEVEL_REF_ID
3119  << ", step " << m_currStep
3120  << ", after " << stepRunTime << " seconds"
3121  << std::endl;
3122  }
3123 
3124  return;
3125 }
3126 //---------------------------------------------------
3127 template <class P_V,class P_M>
3128 void
3129 MLSampling<P_V,P_M>::generateSequence_Step07_inter0(
3130  bool useBalancedChains, // input
3131  unsigned int indexOfFirstWeight, // input
3132  unsigned int indexOfLastWeight, // input
3133  const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // input
3134  UnbalancedLinkedChainsPerNodeStruct& unbalancedLinkControl, // (possible) output
3135  const MLSamplingLevelOptions* currOptions, // input
3136  const SequenceOfVectors<P_V,P_M>& prevChain, // input
3137  double prevExponent, // input
3138  double currExponent, // input
3139  const ScalarSequence<double>& prevLogLikelihoodValues, // input
3140  const ScalarSequence<double>& prevLogTargetValues, // input
3141  std::vector<ExchangeInfoStruct>& exchangeStdVec, // (possible) input/output
3142  BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl) // (possible) output
3143 {
3144  int iRC = UQ_OK_RC;
3145  struct timeval timevalStep;
3146  iRC = gettimeofday(&timevalStep, NULL);
3147  if (iRC) {}; // just to remove compiler warning
3148 
3149  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3150  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3151  << ", level " << m_currLevel+LEVEL_REF_ID
3152  << ", step " << m_currStep
3153  << ": beginning step 7 of 11"
3154  << std::endl;
3155  }
3156 
3157  if (useBalancedChains) {
3158  prepareBalLinkedChains_inter0(currOptions, // input
3159  prevChain, // input
3160  prevExponent, // input
3161  currExponent, // input
3162  prevLogLikelihoodValues, // input
3163  prevLogTargetValues, // input
3164  exchangeStdVec, // input/output
3165  balancedLinkControl); // output
3166  }
3167  else {
3168  prepareUnbLinkedChains_inter0(indexOfFirstWeight, // input
3169  indexOfLastWeight, // input
3170  unifiedIndexCountersAtProc0Only, // input
3171  unbalancedLinkControl); // output
3172  }
3173 
3174  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3175  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3176  << ", level " << m_currLevel+LEVEL_REF_ID
3177  << ", step " << m_currStep
3178  << ": balancedLinkControl.balLinkedChains.size() = " << balancedLinkControl.balLinkedChains.size()
3179  << ", unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.unbLinkedChains.size()
3180  << std::endl;
3181  }
3182 
3183  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3184  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3185  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3186  << ", level " << m_currLevel+LEVEL_REF_ID
3187  << ", step " << m_currStep
3188  << ", after " << stepRunTime << " seconds"
3189  << std::endl;
3190  }
3191 
3192  return;
3193 }
3194 //---------------------------------------------------
3195 template <class P_V,class P_M>
3196 void
3197 MLSampling<P_V,P_M>::generateSequence_Step08_all(
3198  BayesianJointPdf<P_V,P_M>& currPdf, // input/output
3199  GenericVectorRV<P_V,P_M>& currRv) // output
3200 {
3201  int iRC = UQ_OK_RC;
3202  struct timeval timevalStep;
3203  iRC = gettimeofday(&timevalStep, NULL);
3204  if (iRC) {}; // just to remove compiler warning
3205 
3206  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3207  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3208  << ", level " << m_currLevel+LEVEL_REF_ID
3209  << ", step " << m_currStep
3210  << ": beginning step 8 of 11"
3211  << std::endl;
3212  }
3213 
3214  currRv.setPdf(currPdf);
3215 
3216  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3217  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3218  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3219  << ", level " << m_currLevel+LEVEL_REF_ID
3220  << ", step " << m_currStep
3221  << ", after " << stepRunTime << " seconds"
3222  << std::endl;
3223  }
3224 
3225  return;
3226 }
3227 //---------------------------------------------------
3228 template <class P_V,class P_M>
3229 void
3230 MLSampling<P_V,P_M>::generateSequence_Step09_all(
3231  const SequenceOfVectors<P_V,P_M>& prevChain, // input
3232  double prevExponent, // input
3233  double currExponent, // input
3234  const ScalarSequence<double>& prevLogLikelihoodValues, // input
3235  const ScalarSequence<double>& prevLogTargetValues, // input
3236  unsigned int indexOfFirstWeight, // input
3237  unsigned int indexOfLastWeight, // input
3238  const std::vector<double>& unifiedWeightStdVectorAtProc0Only, // input
3239  const ScalarSequence<double>& weightSequence, // input
3240  double prevEta, // input
3241  const GenericVectorRV<P_V,P_M>& currRv, // input
3242  MLSamplingLevelOptions* currOptions, // input (changed temporarily internally)
3243  P_M& unifiedCovMatrix, // input/output
3244  double& currEta) // output
3245 {
3246  int iRC = UQ_OK_RC;
3247  struct timeval timevalStep;
3248  iRC = gettimeofday(&timevalStep, NULL);
3249  if (iRC) {}; // just to remove compiler warning
3250 
3251  if (currOptions->m_scaleCovMatrix == false) {
3252  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3253  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3254  << ", level " << m_currLevel+LEVEL_REF_ID
3255  << ", step " << m_currStep
3256  << ": skipping step 9 of 11"
3257  << std::endl;
3258  }
3259  }
3260  else {
3261  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3262  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3263  << ", level " << m_currLevel+LEVEL_REF_ID
3264  << ", step " << m_currStep
3265  << ": beginning step 9 of 11"
3266  << std::endl;
3267  }
3268 
3269  double beforeEta = prevEta;
3270  double beforeRejectionRate = 0.; // To be updated
3271  bool beforeRejectionRateIsBelowRange = true; // To be updated
3272 
3273  double nowEta = prevEta;
3274  double nowRejectionRate = 0.; // To be computed
3275  bool nowRejectionRateIsBelowRange = true; // To be computed
3276 
3277  std::vector<double> etas(2,0.);
3278  etas[0] = beforeEta;
3279  etas[1] = 1.;
3280 
3281  std::vector<double> rejs(2,0.);
3282  rejs[0] = 0.; // To be computed
3283  rejs[1] = 0.; // To be computed
3284 
3285  unsigned int nowAttempt = 0;
3286  bool testResult = false;
3287  double meanRejectionRate = .5*(currOptions->m_minRejectionRate + currOptions->m_maxRejectionRate);
3288  bool useMiddlePointLogicForEta = false;
3289  P_M nowCovMatrix(unifiedCovMatrix);
3290 #if 0 // KAUST, to check
3291  std::vector<double> unifiedWeightStdVectorAtProc0Only(0);
3292  weightSequence.getUnifiedContentsAtProc0Only(m_vectorSpace.numOfProcsForStorage() == 1,
3293  unifiedWeightStdVectorAtProc0Only);
3294 #endif
3295  do {
3296  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3297  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3298  << ", level " << m_currLevel+LEVEL_REF_ID
3299  << ", step " << m_currStep
3300  << ": entering loop for assessing rejection rate"
3301  << ", with nowAttempt = " << nowAttempt
3302  << ", nowRejectionRate = " << nowRejectionRate
3303  << std::endl;
3304  }
3305  nowCovMatrix = unifiedCovMatrix;
3306 
3307  if (nowRejectionRate < currOptions->m_minRejectionRate) {
3308  nowRejectionRateIsBelowRange = true;
3309  }
3310  else if (nowRejectionRate > currOptions->m_maxRejectionRate) {
3311  nowRejectionRateIsBelowRange = false;
3312  }
3313  else {
3314  queso_error_msg("nowRejectionRate should be out of the requested range at this point of the logic");
3315  }
3316 
3317  if (m_env.inter0Rank() >= 0) { // KAUST
3318  if (nowAttempt > 0) {
3319  if (useMiddlePointLogicForEta == false) {
3320  if (nowAttempt == 1) {
3321  // Ok, keep useMiddlePointLogicForEta = false
3322  }
3323  else if ((beforeRejectionRateIsBelowRange == true) &&
3324  (nowRejectionRateIsBelowRange == true)) {
3325  // Ok
3326  }
3327  else if ((beforeRejectionRateIsBelowRange == false) &&
3328  (nowRejectionRateIsBelowRange == false)) {
3329  // Ok
3330  }
3331  else if ((beforeRejectionRateIsBelowRange == true ) &&
3332  (nowRejectionRateIsBelowRange == false)) {
3333  useMiddlePointLogicForEta = true;
3334 
3335  // This is the first time the middle point logic will be used below
3336  etas[0] = std::min(beforeEta,nowEta);
3337  etas[1] = std::max(beforeEta,nowEta);
3338 
3339  if (etas[0] == beforeEta) {
3340  rejs[0] = beforeRejectionRate;
3341  rejs[1] = nowRejectionRate;
3342  }
3343  else {
3344  rejs[0] = nowRejectionRate;
3345  rejs[1] = beforeRejectionRate;
3346  }
3347  }
3348  else if ((beforeRejectionRateIsBelowRange == false) &&
3349  (nowRejectionRateIsBelowRange == true )) {
3350  useMiddlePointLogicForEta = true;
3351 
3352  // This is the first time the middle point logic will be used below
3353  etas[0] = std::min(beforeEta,nowEta);
3354  etas[1] = std::max(beforeEta,nowEta);
3355  }
3356  else {
3357  queso_error_msg("before and now range flags are inconsistent");
3358  }
3359  } // if (useMiddlePointLogicForEta == false)
3360 
3361  beforeEta = nowEta;
3362  beforeRejectionRate = nowRejectionRate;
3363  beforeRejectionRateIsBelowRange = nowRejectionRateIsBelowRange;
3364  if (useMiddlePointLogicForEta == false) {
3365  if (beforeRejectionRateIsBelowRange) nowEta *= 4.;
3366  else nowEta /= 4.;
3367  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3368  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3369  << ", level " << m_currLevel+LEVEL_REF_ID
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
3375  << std::endl;
3376  }
3377  }
3378  else {
3379  if (nowRejectionRate > meanRejectionRate) {
3380  if (rejs[0] > meanRejectionRate) {
3381  etas[0] = nowEta;
3382  etas[1] = etas[1];
3383  }
3384  else {
3385  etas[0] = etas[0];
3386  etas[1] = nowEta;
3387  }
3388  }
3389  else {
3390  if (rejs[0] < meanRejectionRate) {
3391  etas[0] = nowEta;
3392  etas[1] = etas[1];
3393  }
3394  else {
3395  etas[0] = etas[0];
3396  etas[1] = nowEta;
3397  }
3398  }
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()"
3402  << ", level " << m_currLevel+LEVEL_REF_ID
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]
3410  << std::endl;
3411  }
3412  }
3413  } // if (nowAttempt > 0)
3414  } // if (m_env.inter0Rank() >= 0) // KAUST
3415 
3416  nowCovMatrix *= nowEta;
3417 
3418  // prudencio 2010-12-09: logic 'originalSubNumSamples += 1' added because of the difference of results between GNU and INTEL compiled codes
3419  double doubSubNumSamples = (1.-meanRejectionRate)/meanRejectionRate/currOptions->m_covRejectionRate/currOptions->m_covRejectionRate; // e.g. 19.99...; or 20.0; or 20.1; or 20.9
3420  unsigned int originalSubNumSamples = 1 + (unsigned int) (doubSubNumSamples); // e.g. 20; or 21; or 21; or 21
3421  double auxDouble = (double) originalSubNumSamples; // e.g. 20.0; or 21.0; or 21.0; or 21.0
3422  if ((auxDouble - doubSubNumSamples) < 1.e-8) { // e.g. 0.00...01; or 1.0; or 0.9; or 0.1
3423  originalSubNumSamples += 1;
3424  }
3425 
3426  if (m_env.inter0Rank() >= 0) { // KAUST
3427  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3428  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3429  << ", level " << m_currLevel+LEVEL_REF_ID
3430  << ", step " << m_currStep
3431  << ": in loop for assessing rejection rate"
3432  << ", about to sample " << originalSubNumSamples << " indexes"
3433  << ", meanRejectionRate = " << meanRejectionRate
3434  << ", covRejectionRate = " << currOptions->m_covRejectionRate
3435  << std::endl;
3436  }
3437  } // KAUST
3438 
3439  std::vector<unsigned int> nowUnifiedIndexCountersAtProc0Only(0); // It will be resized by 'sampleIndexes_proc0()' below
3440  if (m_env.inter0Rank() >= 0) { // KAUST
3441  unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3442  sampleIndexes_proc0(tmpUnifiedNumSamples, // input
3443  unifiedWeightStdVectorAtProc0Only, // input
3444  nowUnifiedIndexCountersAtProc0Only); // output
3445 
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");
3449  }
3450 
3451  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3452  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3453  << ", level " << m_currLevel+LEVEL_REF_ID
3454  << ", step " << m_currStep
3455  << ": in loop for assessing rejection rate"
3456  << ", about to distribute sampled assessment indexes"
3457  << std::endl;
3458  }
3459  } // KAUST
3460 
3461  std::vector<ExchangeInfoStruct> exchangeStdVec(0);
3462  BalancedLinkedChainsPerNodeStruct<P_V> nowBalLinkControl;
3463  UnbalancedLinkedChainsPerNodeStruct nowUnbLinkControl; // KAUST
3464 
3465  // All processors should call this routine in order to have the same decision value
3466  bool useBalancedChains = decideOnBalancedChains_all(currOptions, // input
3467  indexOfFirstWeight, // input
3468  indexOfLastWeight, // input
3469  nowUnifiedIndexCountersAtProc0Only, // input
3470  exchangeStdVec); // output
3471 
3472  if (m_env.inter0Rank() >= 0) { // KAUST
3473  if (useBalancedChains) {
3474  prepareBalLinkedChains_inter0(currOptions, // input
3475  prevChain, // input
3476  prevExponent, // input
3477  currExponent, // input
3478  prevLogLikelihoodValues, // input
3479  prevLogTargetValues, // input
3480  exchangeStdVec, // input/output
3481  nowBalLinkControl); // output
3482  }
3483  else {
3484  prepareUnbLinkedChains_inter0(indexOfFirstWeight, // input
3485  indexOfLastWeight, // input
3486  nowUnifiedIndexCountersAtProc0Only, // input
3487  nowUnbLinkControl); // output
3488  }
3489  } // KAUST
3490 
3491  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3492  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3493  << ", level " << m_currLevel+LEVEL_REF_ID
3494  << ", step " << m_currStep
3495  << ": in loop for assessing rejection rate"
3496  << ", about to generate assessment chain"
3497  << std::endl;
3498  }
3499 
3500  SequenceOfVectors<P_V,P_M> nowChain(m_vectorSpace,
3501  0,
3502  m_options.m_prefix+"now_chain");
3503  double nowRunTime = 0.;
3504  unsigned int nowRejections = 0;
3505 
3506  // KAUST: all nodes should call here
3507  bool savedTotallyMute = currOptions->m_totallyMute; // HERE - ENHANCEMENT
3508  unsigned int savedRawChainSize = currOptions->m_rawChainSize; // Ok to use rawChainSize
3509 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3510  bool savedRawChainComputeStats = currOptions->m_rawChainComputeStats;
3511 #endif
3512  bool savedFilteredChainGenerate = currOptions->m_filteredChainGenerate;
3513  unsigned int savedDrMaxNumExtraStages = currOptions->m_drMaxNumExtraStages;
3514  unsigned int savedAmAdaptInterval = currOptions->m_amAdaptInterval;
3515 
3516  currOptions->m_totallyMute = true;
3517  if (m_env.displayVerbosity() >= 999999) {
3518  currOptions->m_totallyMute = false;
3519  }
3520  currOptions->m_rawChainSize = 0; // will be set inside generateXYZLinkedChains()
3521 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3522  currOptions->m_rawChainComputeStats = false;
3523 #endif
3524  currOptions->m_filteredChainGenerate = false;
3525  currOptions->m_drMaxNumExtraStages = 0;
3526  currOptions->m_amAdaptInterval = 0;
3527 
3528  // KAUST: all nodes in 'subComm' should call here, important
3529  if (useBalancedChains) {
3530  generateBalLinkedChains_all(*currOptions, // input, only m_rawChainSize changes
3531  nowCovMatrix, // input
3532  currRv, // input
3533  nowBalLinkControl, // input // Round Rock
3534  nowChain, // output
3535  nowRunTime, // output
3536  nowRejections, // output
3537  NULL, // output
3538  NULL); // output
3539  }
3540  else {
3541  generateUnbLinkedChains_all(*currOptions, // input, only m_rawChainSize changes
3542  nowCovMatrix, // input
3543  currRv, // input
3544  nowUnbLinkControl, // input // Round Rock
3545  indexOfFirstWeight, // input // Round Rock
3546  prevChain, // input // Round Rock
3547  prevExponent, // input
3548  currExponent, // input
3549  prevLogLikelihoodValues, // input
3550  prevLogTargetValues, // input
3551  nowChain, // output
3552  nowRunTime, // output
3553  nowRejections, // output
3554  NULL, // output
3555  NULL); // output
3556  }
3557 
3558  // KAUST: all nodes should call here
3559  currOptions->m_totallyMute = savedTotallyMute;
3560  currOptions->m_rawChainSize = savedRawChainSize;
3561 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3562  currOptions->m_rawChainComputeStats = savedRawChainComputeStats;
3563 #endif
3564  currOptions->m_filteredChainGenerate = savedFilteredChainGenerate; // FIX ME
3565  currOptions->m_drMaxNumExtraStages = savedDrMaxNumExtraStages;
3566  currOptions->m_amAdaptInterval = savedAmAdaptInterval;
3567 
3568  for (unsigned int i = 0; i < nowBalLinkControl.balLinkedChains.size(); ++i) {
3569  queso_require_msg(nowBalLinkControl.balLinkedChains[i].initialPosition, "Initial position pointer in step 9 should not be NULL");
3570  delete nowBalLinkControl.balLinkedChains[i].initialPosition;
3571  nowBalLinkControl.balLinkedChains[i].initialPosition = NULL;
3572  }
3573  nowBalLinkControl.balLinkedChains.clear();
3574 
3575  if (m_env.inter0Rank() >= 0) { // KAUST
3576  // If only one cov matrix is used, then the rejection should be assessed among all inter0Comm nodes // KAUST3
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");
3581 
3582 #if 0 // Round Rock 2009 12 29
3583  unsigned int tmpUnifiedNumSamples = 0;
3584  m_env.inter0Comm().Allreduce((void *) &tmpSubNumSamples, (void *) &tmpUnifiedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
3585  "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3586  "failed MPI.Allreduce() for num samples in step 9");
3587 #endif
3588 
3589  unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3590  nowRejectionRate = ((double) nowUnifiedRejections) / ((double) tmpUnifiedNumSamples);
3591 
3592  //bool aux1 = (nowRejectionRate == meanRejectionRate);
3593  bool aux2 = (nowRejectionRate >= currOptions->m_minRejectionRate)
3594  &&
3595  (nowRejectionRate <= currOptions->m_maxRejectionRate);
3596  testResult = aux2;
3597 
3598  // Make sure all nodes in 'inter0Comm' have the same value of 'testResult'
3599  if (MiscCheckForSameValueInAllNodes(testResult,
3600  0., // kept 'zero' on 2010/03/03
3601  m_env.inter0Comm(),
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()"
3605  << ", level " << m_currLevel+LEVEL_REF_ID
3606  << ", step " << m_currStep
3607  << ": nowAttempt = " << nowAttempt
3608  << ", MiscCheck for 'testResult' detected a problem"
3609  << std::endl;
3610  }
3611  }
3612  } // if (m_env.inter0Rank() >= 0) { // KAUST
3613 
3614  // KAUST: all nodes in 'subComm' should have the same 'testResult'
3615  unsigned int tmpUint = (unsigned int) testResult;
3616  m_env.subComm().Bcast((void *) &tmpUint, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'subComm', important
3617  "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3618  "failed MPI.Bcast() for testResult");
3619  testResult = (bool) tmpUint;
3620 
3621  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3622  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3623  << ", level " << m_currLevel+LEVEL_REF_ID
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]
3631  << ", minRejectionRate = " << currOptions->m_minRejectionRate
3632  << ", nowRejectionRate = " << nowRejectionRate
3633  << ", maxRejectionRate = " << currOptions->m_maxRejectionRate
3634  << std::endl;
3635  }
3636  nowAttempt++;
3637 
3638  if (m_env.inter0Rank() >= 0) { // KAUST
3639  // Make sure all nodes in 'inter0Comm' have the same value of 'nowEta'
3640  if (MiscCheckForSameValueInAllNodes(nowEta,
3641  1.0e-16, // changed from 'zero' on 2009/11/dd
3642  m_env.inter0Comm(),
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()"
3646  << ", level " << m_currLevel+LEVEL_REF_ID
3647  << ", step " << m_currStep
3648  << ": nowAttempt = " << nowAttempt
3649  << ", MiscCheck for 'nowEta' detected a problem"
3650  << std::endl;
3651  }
3652  }
3653  }
3654  } while (testResult == false);
3655  currEta = nowEta;
3656  if (currEta != 1.) {
3657  unifiedCovMatrix *= currEta;
3658  if (m_numDisabledParameters > 0) { // gpmsa2
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.;
3663  }
3664  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
3665  unifiedCovMatrix(paramId,j) = 0.;
3666  }
3667  unifiedCovMatrix(paramId,paramId) = 1.;
3668  }
3669  }
3670  }
3671  }
3672 
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()"
3676  << ", level " << m_currLevel+LEVEL_REF_ID
3677  << ", step " << m_currStep
3678  << ": weightSequence.subSequenceSize() = " << weightSequence.subSequenceSize()
3679  << ", weightSequence.unifiedSequenceSize() = " << quantity1
3680  << ", currEta = " << currEta
3681  << ", assessed rejection rate = " << nowRejectionRate
3682  << std::endl;
3683  }
3684  }
3685 
3686  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3687  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3688  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3689  << ", level " << m_currLevel+LEVEL_REF_ID
3690  << ", step " << m_currStep
3691  << ", after " << stepRunTime << " seconds"
3692  << std::endl;
3693  }
3694 
3695  return;
3696 }
3697 //---------------------------------------------------
3698 template <class P_V,class P_M>
3699 void
3700 MLSampling<P_V,P_M>::generateSequence_Step10_all(
3701  MLSamplingLevelOptions& currOptions, // input (changed temporarily internally)
3702  const P_M& unifiedCovMatrix, // input
3703  const GenericVectorRV <P_V,P_M>& currRv, // input
3704  bool useBalancedChains, // input
3705  const UnbalancedLinkedChainsPerNodeStruct& unbalancedLinkControl, // input // Round Rock
3706  unsigned int indexOfFirstWeight, // input // Round Rock
3707  const SequenceOfVectors<P_V,P_M>& prevChain, // input // Round Rock
3708  double prevExponent, // input
3709  double currExponent, // input
3710  const ScalarSequence<double>& prevLogLikelihoodValues, // input
3711  const ScalarSequence<double>& prevLogTargetValues, // input
3712  const BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl, // input // Round Rock
3713  SequenceOfVectors <P_V,P_M>& currChain, // output
3714  double& cumulativeRawChainRunTime, // output
3715  unsigned int& cumulativeRawChainRejections, // output
3716  ScalarSequence <double>* currLogLikelihoodValues, // output
3717  ScalarSequence <double>* currLogTargetValues) // output
3718 {
3719  int iRC = UQ_OK_RC;
3720  struct timeval timevalStep;
3721  iRC = gettimeofday(&timevalStep, NULL);
3722  if (iRC) {}; // just to remove compiler warning
3723 
3724  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3725  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3726  << ", level " << m_currLevel+LEVEL_REF_ID
3727  << ", step " << m_currStep
3728  << ": beginning step 10 of 11"
3729  << ", currLogLikelihoodValues = " << currLogLikelihoodValues
3730  << std::endl;
3731  }
3732 
3733  // All nodes should call here
3734  bool savedTotallyMute = currOptions.m_totallyMute; // HERE - ENHANCEMENT
3735  unsigned int savedRawChainSize = currOptions.m_rawChainSize; // Ok to use rawChainSize
3736 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3737  bool savedRawChainComputeStats = currOptions.m_rawChainComputeStats;
3738 #endif
3739  bool savedFilteredChainGenerate = currOptions.m_filteredChainGenerate;
3740 
3741  currOptions.m_totallyMute = true;
3742  if (m_env.displayVerbosity() >= 999999) {
3743  currOptions.m_totallyMute = false;
3744  }
3745  currOptions.m_rawChainSize = 0; // will be set inside generateXYZLinkedChains()
3746 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3747  currOptions.m_rawChainComputeStats = false;
3748 #endif
3749  currOptions.m_filteredChainGenerate = false;
3750 
3751  // All nodes should call here
3752  if (useBalancedChains) {
3753  generateBalLinkedChains_all(currOptions, // input, only m_rawChainSize changes
3754  unifiedCovMatrix, // input
3755  currRv, // input
3756  balancedLinkControl, // input // Round Rock
3757  currChain, // output
3758  cumulativeRawChainRunTime, // output
3759  cumulativeRawChainRejections, // output
3760  currLogLikelihoodValues, // output // likelihood is important
3761  currLogTargetValues); // output
3762  }
3763  else {
3764  generateUnbLinkedChains_all(currOptions, // input, only m_rawChainSize changes
3765  unifiedCovMatrix, // input
3766  currRv, // input
3767  unbalancedLinkControl, // input // Round Rock
3768  indexOfFirstWeight, // input // Round Rock
3769  prevChain, // input // Round Rock
3770  prevExponent, // input
3771  currExponent, // input
3772  prevLogLikelihoodValues, // input
3773  prevLogTargetValues, // input
3774  currChain, // output
3775  cumulativeRawChainRunTime, // output
3776  cumulativeRawChainRejections, // output
3777  currLogLikelihoodValues, // output // likelihood is important
3778  currLogTargetValues); // output
3779  }
3780 
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()"
3785  << ", level " << m_currLevel+LEVEL_REF_ID
3786  << ", step " << m_currStep
3787  << ", after chain generatrion"
3788  << ", currLogLikelihoodValues[0] = " << tmpValue
3789  << std::endl;
3790  }
3791 
3792  // All nodes should call here
3793  currOptions.m_totallyMute = savedTotallyMute;
3794  currOptions.m_rawChainSize = savedRawChainSize;
3795 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3796  currOptions.m_rawChainComputeStats = savedRawChainComputeStats;
3797 #endif
3798  currOptions.m_filteredChainGenerate = savedFilteredChainGenerate; // FIX ME
3799 
3800  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3801  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3802  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3803  << ", level " << m_currLevel+LEVEL_REF_ID
3804  << ", step " << m_currStep
3805  << ", after " << stepRunTime << " seconds"
3806  << std::endl;
3807  }
3808 
3809  return;
3810 }
3811 //---------------------------------------------------
3812 template <class P_V,class P_M>
3813 void
3814 MLSampling<P_V,P_M>::generateSequence_Step11_inter0(
3815  const MLSamplingLevelOptions* currOptions, // input
3816  unsigned int unifiedRequestedNumSamples, // input
3817  unsigned int cumulativeRawChainRejections, // input
3818  SequenceOfVectors<P_V,P_M>& currChain, // input/output
3819  ScalarSequence<double>& currLogLikelihoodValues, // input/output
3820  ScalarSequence<double>& currLogTargetValues, // input/output
3821  unsigned int& unifiedNumberOfRejections) // output
3822 {
3823  int iRC = UQ_OK_RC;
3824  struct timeval timevalStep;
3825  iRC = gettimeofday(&timevalStep, NULL);
3826  if (iRC) {}; // just to remove compiler warning
3827 
3828  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3829  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3830  << ", level " << m_currLevel+LEVEL_REF_ID
3831  << ", step " << m_currStep
3832  << ": beginning step 11 of 11"
3833  << std::endl;
3834  }
3835 
3836  //if (m_env.subComm().MyPID() == 0) std::cout << "Aqui 000" << std::endl;
3837 
3838 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3839  if (currOptions->m_rawChainComputeStats) {
3840  FilePtrSetStruct filePtrSet;
3841  m_env.openOutputFile(currOptions->m_dataOutputFileName,
3842  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
3843  currOptions->m_dataOutputAllowedSet,
3844  false,
3845  filePtrSet);
3846 
3847  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) { // output debug
3848  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3849  << ", level " << m_currLevel+LEVEL_REF_ID
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
3855  << std::endl;
3856  }
3857  //m_env.syncPrintDebugMsg("At step 11, calling computeStatistics for raw chain",1,10,m_env.inter0Comm()); // output debug
3858  currChain.computeStatistics(*currOptions->m_rawChainStatisticalOptionsObj,
3859  filePtrSet.ofsVar);
3860 
3861  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
3862  }
3863  // Compute MLE and MAP
3864  // rr0
3865 #endif
3866  if (currOptions->m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
3867  currChain.unifiedWriteContents(currOptions->m_rawChainDataOutputFileName,
3868  currOptions->m_rawChainDataOutputFileType); // KAUST5
3869  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3870  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3871  << ", level " << m_currLevel+LEVEL_REF_ID
3872  << ", step " << m_currStep
3873  << ", before calling currLogLikelihoodValues.unifiedWriteContents()"
3874  << ", currLogLikelihoodValues[0] = " << currLogLikelihoodValues[0]
3875  << std::endl;
3876  }
3877  currLogLikelihoodValues.unifiedWriteContents(currOptions->m_rawChainDataOutputFileName,
3878  currOptions->m_rawChainDataOutputFileType);
3879  currLogTargetValues.unifiedWriteContents (currOptions->m_rawChainDataOutputFileName,
3880  currOptions->m_rawChainDataOutputFileType);
3881  }
3882 
3883  if (currOptions->m_filteredChainGenerate) {
3884  FilePtrSetStruct filePtrSet;
3885  m_env.openOutputFile(currOptions->m_dataOutputFileName,
3886  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
3887  currOptions->m_dataOutputAllowedSet,
3888  false,
3889  filePtrSet);
3890 
3891  unsigned int filterInitialPos = (unsigned int) (currOptions->m_filteredChainDiscardedPortion * (double) currChain.subSequenceSize());
3892  unsigned int filterSpacing = currOptions->m_filteredChainLag;
3893  if (filterSpacing == 0) {
3894  currChain.computeFilterParams(filePtrSet.ofsVar,
3895  filterInitialPos,
3896  filterSpacing);
3897  }
3898 
3899  // Filter positions from the converged portion of the chain
3900  currChain.filter(filterInitialPos,
3901  filterSpacing);
3902  currChain.setName(currOptions->m_prefix + "filtChain");
3903 
3904  currLogLikelihoodValues.filter(filterInitialPos,
3905  filterSpacing);
3906  currLogLikelihoodValues.setName(currOptions->m_prefix + "filtLogLikelihood");
3907 
3908  currLogTargetValues.filter(filterInitialPos,
3909  filterSpacing);
3910  currLogTargetValues.setName(currOptions->m_prefix + "filtLogTarget");
3911 
3912 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3913  if (currOptions->m_filteredChainComputeStats) {
3914  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) { // output debug
3915  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3916  << ", level " << m_currLevel+LEVEL_REF_ID
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
3922  << std::endl;
3923  }
3924 
3925  //m_env.syncPrintDebugMsg("At step 11, calling computeStatistics for filtered chain",1,10,m_env.inter0Comm()); // output debug
3926  currChain.computeStatistics(*currOptions->m_filteredChainStatisticalOptionsObj,
3927  filePtrSet.ofsVar);
3928 
3929  }
3930 #endif
3931  // Compute MLE and MAP
3932  // rr0
3933  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
3934 
3935  if (currOptions->m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
3936  currChain.unifiedWriteContents (currOptions->m_filteredChainDataOutputFileName,
3937  currOptions->m_filteredChainDataOutputFileType);
3938  currLogLikelihoodValues.unifiedWriteContents(currOptions->m_filteredChainDataOutputFileName,
3939  currOptions->m_filteredChainDataOutputFileType);
3940  currLogTargetValues.unifiedWriteContents (currOptions->m_filteredChainDataOutputFileName,
3941  currOptions->m_filteredChainDataOutputFileType);
3942  }
3943  } // if (currOptions->m_filteredChainGenerate)
3944 
3945  if (currOptions->m_filteredChainGenerate) {
3946  // Do not check
3947  }
3948  else {
3949  // Check if unified size of generated chain matches the unified requested size // KAUST
3950  unsigned int tmpSize = currChain.subSequenceSize();
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");
3955  //std::cout << "unifiedGeneratedNumSamples = " << unifiedGeneratedNumSamples
3956  // << ", unifiedRequestedNumSamples = " << unifiedRequestedNumSamples
3957  // << std::endl;
3958  queso_require_equal_to_msg(unifiedGeneratedNumSamples, unifiedRequestedNumSamples, "currChain (linked one) has been generated with invalid size");
3959  }
3960 
3961  // Compute unified number of rejections
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");
3965 
3966  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3967  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3968  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3969  << ", level " << m_currLevel+LEVEL_REF_ID
3970  << ", step " << m_currStep
3971  << ", after " << stepRunTime << " seconds"
3972  << std::endl;
3973  }
3974 
3975  return;
3976 }
3977 
3978 // Default constructor -----------------------------
3979 template<class P_V,class P_M>
3980 MLSampling<P_V,P_M>::MLSampling(
3981  const char* prefix,
3982  const BaseVectorRV <P_V,P_M>& priorRv,
3983  const BaseScalarFunction<P_V,P_M>& likelihoodFunction)
3984  :
3985  m_env (priorRv.env()),
3986  m_priorRv (priorRv),
3987  m_likelihoodFunction(likelihoodFunction),
3988  m_vectorSpace (m_priorRv.imageSet().vectorSpace()),
3989  m_targetDomain (InstantiateIntersection(m_priorRv.pdf().domainSet(),m_likelihoodFunction.domainSet())),
3990  m_numDisabledParameters (0), // gpmsa2
3991  m_parameterEnabledStatus(m_vectorSpace.dimLocal(),true), // gpmsa2
3992  m_options (m_env,prefix),
3993  m_currLevel (0),
3994  m_currStep (0),
3995  m_debugExponent (0.),
3996  m_logEvidenceFactors(0),
3997  m_logEvidence (0.),
3998  m_meanLogLikelihood (0.),
3999  m_eig (0.)
4000 {
4001  if (m_env.subDisplayFile()) {
4002  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::constructor()"
4003  << std::endl;
4004  }
4005 
4006  if (m_env.subDisplayFile()) {
4007  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::constructor()"
4008  << std::endl;
4009  }
4010 }
4011 // Destructor ---------------------------------------
4012 template<class P_V,class P_M>
4013 MLSampling<P_V,P_M>::~MLSampling()
4014 {
4015  m_numDisabledParameters = 0; // gpmsa2
4016  m_parameterEnabledStatus.clear(); // gpmsa2
4017  if (m_targetDomain) delete m_targetDomain;
4018 }
4019 // Statistical methods-------------------------------
4020 /* This operation currently implements the PAMSSA algorithm (S. H. Cheung and E. E. Prudencio. Parallel adaptive multilevel
4021  * sampling algorithms for the Bayesian analysis of mathematical models. International Journal
4022  * for Uncertainty Quantification, 2(3):215237, 2012.)*/
4023 template <class P_V,class P_M>
4024 void
4025 MLSampling<P_V,P_M>::generateSequence(
4026  BaseVectorSequence<P_V,P_M>& workingChain,
4027  ScalarSequence<double>* workingLogLikelihoodValues,
4028  ScalarSequence<double>* workingLogTargetValues)
4029 {
4030  struct timeval timevalRoutineBegin;
4031  int iRC = 0;
4032  iRC = gettimeofday(&timevalRoutineBegin, NULL);
4033  if (iRC) {}; // just to remove compiler warning
4034 
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..."
4040  << std::endl;
4041  }
4042 
4043  //***********************************************************
4044  // Declaration of Variables
4045  //***********************************************************
4046  double currExponent = 0.; // restate
4047  double currEta = 1.; // restate
4048  unsigned int currUnifiedRequestedNumSamples = 0;
4049  SequenceOfVectors<P_V,P_M> currChain (m_vectorSpace, // restate
4050  0,
4051  m_options.m_prefix+"curr_chain");
4052  ScalarSequence<double> currLogLikelihoodValues(m_env, // restate
4053  0,
4054  "");
4055  ScalarSequence<double> currLogTargetValues (m_env, // restate
4056  0,
4057  "");
4058 
4059  bool stopAtEndOfLevel = false;
4060  char levelPrefix[256];
4061 
4062  //***********************************************************
4063  // Take care of first level (level '0')
4064  //***********************************************************
4065  MLSamplingLevelOptions defaultLevelOptions(m_env,(m_options.m_prefix + "default_").c_str());
4066  defaultLevelOptions.scanOptionsValues(NULL);
4067 
4068  MLSamplingLevelOptions lastLevelOptions(m_env,(m_options.m_prefix + "last_").c_str());
4069  lastLevelOptions.scanOptionsValues(&defaultLevelOptions);
4070 
4071  if (m_options.m_restartInput_baseNameForFiles != ".") {
4072  restartML(currExponent, // output
4073  currEta, // output
4074  currChain, // output
4075  currLogLikelihoodValues, // output
4076  currLogTargetValues); // output
4077 
4078  if (currExponent == 1.) {
4079  if (lastLevelOptions.m_parameterDisabledSet.size() > 0) { // gpmsa2
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;
4085  }
4086  }
4087  }
4088 
4089 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4090  if (lastLevelOptions.m_rawChainComputeStats) {
4091  FilePtrSetStruct filePtrSet;
4092  m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4093  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
4094  lastLevelOptions.m_dataOutputAllowedSet,
4095  false,
4096  filePtrSet);
4097 
4098  currChain.computeStatistics(*lastLevelOptions.m_rawChainStatisticalOptionsObj,
4099  filePtrSet.ofsVar);
4100 
4101  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
4102  }
4103 #endif
4104  // Compute MLE and MAP
4105  // rr0
4106 
4107  if (lastLevelOptions.m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
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);
4111  }
4112 
4113  if (lastLevelOptions.m_filteredChainGenerate) {
4114  FilePtrSetStruct filePtrSet;
4115  m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4116  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
4117  lastLevelOptions.m_dataOutputAllowedSet,
4118  false,
4119  filePtrSet);
4120 
4121  unsigned int filterInitialPos = (unsigned int) (lastLevelOptions.m_filteredChainDiscardedPortion * (double) currChain.subSequenceSize());
4122  unsigned int filterSpacing = lastLevelOptions.m_filteredChainLag;
4123  if (filterSpacing == 0) {
4124  currChain.computeFilterParams(filePtrSet.ofsVar,
4125  filterInitialPos,
4126  filterSpacing);
4127  }
4128 
4129  // Filter positions from the converged portion of the chain
4130  currChain.filter(filterInitialPos,
4131  filterSpacing);
4132  currChain.setName(lastLevelOptions.m_prefix + "filtChain");
4133 
4134  currLogLikelihoodValues.filter(filterInitialPos,
4135  filterSpacing);
4136  currLogLikelihoodValues.setName(lastLevelOptions.m_prefix + "filtLogLikelihood");
4137 
4138  currLogTargetValues.filter(filterInitialPos,
4139  filterSpacing);
4140  currLogTargetValues.setName(lastLevelOptions.m_prefix + "filtLogTarget");
4141 
4142 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4143  if (lastLevelOptions.m_filteredChainComputeStats) {
4144  currChain.computeStatistics(*lastLevelOptions.m_filteredChainStatisticalOptionsObj,
4145  filePtrSet.ofsVar);
4146  }
4147 #endif
4148  // Compute MLE and MAP
4149  // rr0
4150  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
4151 
4152  if (lastLevelOptions.m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
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);
4156  }
4157  } // if (lastLevelOptions.m_filteredChainGenerate)
4158  }
4159  }
4160  else {
4161  sprintf(levelPrefix,"%d_",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
4162  MLSamplingLevelOptions currOptions(m_env,(m_options.m_prefix + levelPrefix).c_str());
4163  currOptions.scanOptionsValues(&defaultLevelOptions);
4164 
4165  if (currOptions.m_parameterDisabledSet.size() > 0) { // gpmsa2
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;
4171  }
4172  }
4173  }
4174 
4175  generateSequence_Level0_all(currOptions, // input
4176  currUnifiedRequestedNumSamples, // output
4177  currChain, // output
4178  currLogLikelihoodValues, // output
4179  currLogTargetValues); // output
4180 
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 );
4185  }
4186  if (performCheckpoint) {
4187  checkpointML(currExponent, // input
4188  currEta, // input
4189  currChain, // input
4190  currLogLikelihoodValues, // input
4191  currLogTargetValues); // input
4192  }
4193  }
4194  //std::cout << "In QUESO: end of level 0. Exiting on purpose" << std::endl;
4195  //exit(1);
4196 
4197  double minLogLike = -INFINITY;
4198  double maxLogLike = INFINITY;
4199  currLogLikelihoodValues.subMinMaxExtra(0,
4200  currLogLikelihoodValues.subSequenceSize(),
4201  minLogLike,
4202  maxLogLike);
4203  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4204  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4205  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4206  << ", sub minLogLike = " << minLogLike
4207  << ", sub maxLogLike = " << maxLogLike
4208  << std::endl;
4209  }
4210 
4211  m_env.fullComm().Barrier();
4212 
4213  minLogLike = -INFINITY;
4214  maxLogLike = INFINITY;
4215  currLogLikelihoodValues.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4216  0,
4217  currLogLikelihoodValues.subSequenceSize(),
4218  minLogLike,
4219  maxLogLike);
4220  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4221  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4222  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4223  << ", unified minLogLike = " << minLogLike
4224  << ", unified maxLogLike = " << maxLogLike
4225  << std::endl;
4226  }
4227 
4228  //***********************************************************
4229  // Take care of next levels
4230  //***********************************************************
4231  while ((currExponent < 1. ) && // begin level while
4232  (stopAtEndOfLevel == false)) {
4233  m_currLevel++; // restate
4234 
4235  struct timeval timevalLevelBegin;
4236  iRC = 0;
4237  iRC = gettimeofday(&timevalLevelBegin, NULL);
4238 
4239  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4240  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4241  << ": beginning level " << m_currLevel+LEVEL_REF_ID
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"
4247  << std::endl;
4248  }
4249 
4250  iRC = UQ_OK_RC;
4251  struct timeval timevalLevel;
4252  iRC = gettimeofday(&timevalLevel, NULL);
4253  double cumulativeRawChainRunTime = 0.;
4254  unsigned int cumulativeRawChainRejections = 0;
4255 
4256  bool tryExponentEta = true; // gpmsa1
4257  double failedExponent = 0.; // gpmsa1
4258  double failedEta = 0.; // gpmsa1
4259 
4260  MLSamplingLevelOptions* currOptions = NULL; // step 1
4261  SequenceOfVectors<P_V,P_M>* prevChain = NULL; // step 2
4262  ScalarSequence<double> prevLogLikelihoodValues(m_env, 0, "");
4263  ScalarSequence<double> prevLogTargetValues (m_env, 0, "");
4264  double prevExponent = 0.; // step 3
4265  unsigned int indexOfFirstWeight = 0; // step 2
4266  unsigned int indexOfLastWeight = 0; // step 2
4267  P_M* unifiedCovMatrix = NULL; // step 4
4268  bool useBalancedChains = false; // step 6
4269  BalancedLinkedChainsPerNodeStruct<P_V>* balancedLinkControl = NULL; // step 7
4270  UnbalancedLinkedChainsPerNodeStruct* unbalancedLinkControl = NULL; // step 7
4271  BayesianJointPdf<P_V,P_M>* currPdf = NULL; // step 8
4272  GenericVectorRV<P_V,P_M>* currRv = NULL; // step 8
4273 
4274  unsigned int exponentEtaTriedAmount = 0;
4275  while (tryExponentEta) { // gpmsa1
4276  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4277  *m_env.subDisplayFile() << "In IMLSampling<P_V,P_M>::generateSequence()"
4278  << ", level " << m_currLevel+LEVEL_REF_ID
4279  << ", beginning 'do-while(tryExponentEta)"
4280  << ": failedExponent = " << failedExponent
4281  << ", failedEta = " << failedEta
4282  << std::endl;
4283  }
4284 
4285  //***********************************************************
4286  // Step 1 of 11: read options
4287  //***********************************************************
4288  m_currStep = 1;
4289  sprintf(levelPrefix,"%d_",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
4290  currOptions = new MLSamplingLevelOptions(m_env,(m_options.m_prefix + levelPrefix).c_str());
4291  currOptions->scanOptionsValues(&defaultLevelOptions);
4292 
4293  if (currOptions->m_parameterDisabledSet.size() > 0) { // gpmsa2
4294  for (std::set<unsigned int>::iterator setIt = currOptions->m_parameterDisabledSet.begin(); setIt != currOptions->m_parameterDisabledSet.end(); ++setIt) {
4295  unsigned int paramId = *setIt;
4296  if (paramId < m_vectorSpace.dimLocal()) {
4297  m_numDisabledParameters++;
4298  m_parameterEnabledStatus[paramId] = false;
4299  }
4300  }
4301  }
4302 
4303  if (m_env.inter0Rank() >= 0) {
4304  generateSequence_Step01_inter0(currOptions, // input
4305  currUnifiedRequestedNumSamples); // output
4306  }
4307 
4308  //***********************************************************
4309  // Step 2 of 11: save [chain and corresponding target pdf values] from previous level
4310  //***********************************************************
4311  m_currStep = 2;
4312  prevExponent = currExponent;
4313  double prevEta = currEta;
4314  unsigned int prevUnifiedRequestedNumSamples = currUnifiedRequestedNumSamples;
4315  prevChain = new SequenceOfVectors<P_V,P_M>(m_vectorSpace,
4316  0,
4317  m_options.m_prefix+"prev_chain");
4318 
4319  indexOfFirstWeight = 0;
4320  indexOfLastWeight = 0;
4321 
4322  //std::cout << "m_env.inter0Rank() = " << m_env.inter0Rank() << std::endl;
4323  if (m_env.inter0Rank() >= 0) {
4324  generateSequence_Step02_inter0(currOptions, // input
4325  currChain, // input/output // restate
4326  currLogLikelihoodValues, // input/output // restate
4327 
4328  currLogTargetValues, // input/output // restate
4329  *prevChain, // output
4330  prevLogLikelihoodValues, // output
4331  prevLogTargetValues, // output
4332  indexOfFirstWeight, // output
4333  indexOfLastWeight); // output
4334  }
4335 
4336  //***********************************************************
4337  // Step 3 of 11: compute [currExponent and sequence of weights] for current level
4338  // update 'm_logEvidenceFactors'
4339  //***********************************************************
4340  m_currStep = 3;
4341  ScalarSequence<double> weightSequence(m_env,prevLogLikelihoodValues.subSequenceSize(),"");
4342  if (m_env.inter0Rank() >= 0) {
4343  generateSequence_Step03_inter0(currOptions, // input
4344  prevLogLikelihoodValues, // input
4345  prevExponent, // input
4346  failedExponent, // input // gpmsa1
4347  currExponent, // output
4348  weightSequence); // output
4349  }
4350 
4351  // All nodes in 'subComm' should have the same 'currExponent'
4352  m_env.subComm().Bcast((void *) &currExponent, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm', important
4353  "MLSampling<P_V,P_M>::generateSequence()",
4354  "failed MPI.Bcast() for currExponent");
4355  m_debugExponent = currExponent;
4356 
4357  if (currExponent == 1.) {
4358  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4359  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4360  << ", level " << m_currLevel+LEVEL_REF_ID
4361  << ", step " << m_currStep
4362  << ": copying 'last' level options to current options" // In all nodes of 'subComm', important
4363  << std::endl;
4364  }
4365  delete currOptions;
4366  currOptions = &lastLevelOptions;
4367 
4368  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4369  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4370  << ", level " << m_currLevel+LEVEL_REF_ID
4371  << ", step " << m_currStep
4372  << ": after copying 'last' level options to current options, the current options are"
4373  << "\n" << *currOptions
4374  << std::endl;
4375  }
4376 
4377  if (m_env.inter0Rank() >= 0) {
4378  // It is necessary to recompute 'currUnifiedRequestedNumSamples' because
4379  // 'currOptions' has just been replaced by 'lastLevelOptions'
4380  unsigned int tmpSize = currOptions->m_rawChainSize;
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");
4384  }
4385  }
4386 
4387  //***********************************************************
4388  // Step 4 of 11: create covariance matrix for current level
4389  //***********************************************************
4390  m_currStep = 4;
4391  P_V oneVec(m_vectorSpace.zeroVector());
4392  oneVec.cwSet(1.);
4393 
4394  unifiedCovMatrix = NULL;
4395  if (m_env.inter0Rank() >= 0) {
4396  unifiedCovMatrix = m_vectorSpace.newMatrix();
4397  }
4398  else {
4399  unifiedCovMatrix = new P_M(oneVec);
4400  }
4401 
4402  if (m_env.inter0Rank() >= 0) {
4403  generateSequence_Step04_inter0(*prevChain, // input
4404  weightSequence, // input
4405  *unifiedCovMatrix); // output
4406  }
4407 
4408  //***********************************************************
4409  // Step 5 of 11: create *unified* finite distribution for current level
4410  //***********************************************************
4411  m_currStep = 5;
4412  std::vector<unsigned int> unifiedIndexCountersAtProc0Only(0);
4413  std::vector<double> unifiedWeightStdVectorAtProc0Only(0); // KAUST, to check
4414  if (m_env.inter0Rank() >= 0) {
4415  generateSequence_Step05_inter0(currUnifiedRequestedNumSamples, // input
4416  weightSequence, // input
4417  unifiedIndexCountersAtProc0Only, // output
4418  unifiedWeightStdVectorAtProc0Only); // output
4419  }
4420 
4421  //***********************************************************
4422  // Step 6 of 11: decide on using balanced chains or not
4423  //***********************************************************
4424  m_currStep = 6;
4425  useBalancedChains = false;
4426  std::vector<ExchangeInfoStruct> exchangeStdVec(0);
4427  // All processors should call this routine in order to have the same decision value
4428  generateSequence_Step06_all(currOptions, // input
4429  indexOfFirstWeight, // input
4430  indexOfLastWeight, // input
4431  unifiedIndexCountersAtProc0Only, // input
4432  useBalancedChains, // output
4433  exchangeStdVec); // output
4434 
4435  //***********************************************************
4436  // Step 7 of 11: plan for number of linked chains for each node so that all
4437  // nodes generate the closest possible to the same number of positions
4438  //***********************************************************
4439  m_currStep = 7;
4440  balancedLinkControl = new BalancedLinkedChainsPerNodeStruct<P_V>();
4441  unbalancedLinkControl = new UnbalancedLinkedChainsPerNodeStruct ();
4442  if (m_env.inter0Rank() >= 0) {
4443  generateSequence_Step07_inter0(useBalancedChains, // input
4444  indexOfFirstWeight, // input
4445  indexOfLastWeight, // input
4446  unifiedIndexCountersAtProc0Only, // input
4447  *unbalancedLinkControl, // (possible) output
4448  currOptions, // input
4449  *prevChain, // input
4450  prevExponent, // input
4451  currExponent, // input
4452  prevLogLikelihoodValues, // input
4453  prevLogTargetValues, // input
4454  exchangeStdVec, // (possible) input/output
4455  *balancedLinkControl); // (possible) output
4456  }
4457 
4458  // aqui: prevChain might not be needed anymore. Delete it to save memory
4459 
4460  //***********************************************************
4461  // Step 8 of 11: create vector RV for current level
4462  //***********************************************************
4463  m_currStep = 8;
4464  currPdf = new BayesianJointPdf<P_V,P_M> (m_options.m_prefix.c_str(),
4465  m_priorRv.pdf(),
4466  m_likelihoodFunction,
4467  currExponent,
4468  *m_targetDomain);
4469 
4470  currRv = new GenericVectorRV<P_V,P_M> (m_options.m_prefix.c_str(),
4471  *m_targetDomain);
4472 
4473  // All nodes should set 'currRv'
4474  generateSequence_Step08_all(*currPdf,
4475  *currRv);
4476 
4477  //***********************************************************
4478  // Step 9 of 11: scale unified covariance matrix until min <= rejection rate <= max
4479  //***********************************************************
4480  m_currStep = 9;
4481  generateSequence_Step09_all(*prevChain, // input
4482  prevExponent, // input
4483  currExponent, // input
4484  prevLogLikelihoodValues, // input
4485  prevLogTargetValues, // input
4486  indexOfFirstWeight, // input
4487  indexOfLastWeight, // input
4488  unifiedWeightStdVectorAtProc0Only, // input
4489  weightSequence, // input
4490  prevEta, // input
4491  *currRv, // input
4492  currOptions, // input (changed temporarily internally)
4493  *unifiedCovMatrix, // input/output
4494  currEta); // output
4495 
4496  tryExponentEta = false; // gpmsa1
4497  if ((currOptions->m_minAcceptableEta > 0. ) && // gpmsa1
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()"
4501  << ", level " << m_currLevel+LEVEL_REF_ID
4502  << ", preparing to retry ExponentEta"
4503  << ": currExponent = " << currExponent
4504  << ", currEta = " << currEta
4505  << std::endl;
4506  }
4507  exponentEtaTriedAmount++;
4508  tryExponentEta = true;
4509  failedExponent = currExponent;
4510  failedEta = currEta;
4511 
4512  // "Return" to previous level
4513  delete currRv; // Step 8
4514  currRv = NULL; // Step 8
4515  delete currPdf; // Step 8
4516  currPdf = NULL; // Step 8
4517 
4518  delete balancedLinkControl; // Step 7
4519  balancedLinkControl = NULL; // Step 7
4520  delete unbalancedLinkControl; // Step 7
4521  unbalancedLinkControl = NULL; // Step 7
4522 
4523  delete unifiedCovMatrix; // Step 4
4524  unifiedCovMatrix = NULL; // Step 4
4525 
4526  currExponent = prevExponent;
4527  currEta = 1.; // prevEta;
4528  currUnifiedRequestedNumSamples = prevUnifiedRequestedNumSamples;
4529 
4530  currChain.clear(); // Step 2
4531  currChain = (*prevChain); // Step 2
4532  delete prevChain; // Step 2
4533  prevChain = NULL; // Step 2
4534 
4535  currLogLikelihoodValues = prevLogLikelihoodValues;
4536  currLogTargetValues = prevLogTargetValues;
4537  }
4538  } // while (tryExponentEta) // gpmsa1
4539 
4540  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) { // gpmsa1
4541  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4542  << ", level " << m_currLevel+LEVEL_REF_ID
4543  << ", exited 'do-while(tryExponentEta)"
4544  << ", failedExponent = " << failedExponent
4545  << ", failedEta = " << failedEta
4546  << std::endl;
4547  }
4548 
4549  //***********************************************************
4550  // Step 10 of 11: sample vector RV of current level
4551  //***********************************************************
4552  m_currStep = 10;
4553  // All nodes should call here
4554  generateSequence_Step10_all(*currOptions, // input (changed temporarily internally)
4555  *unifiedCovMatrix, // input
4556  *currRv, // input
4557  useBalancedChains, // input
4558  *unbalancedLinkControl, // input // Round Rock
4559  indexOfFirstWeight, // input // Round Rock
4560  *prevChain, // input // Round Rock
4561  prevExponent, // input
4562  currExponent, // input
4563  prevLogLikelihoodValues, // input
4564  prevLogTargetValues, // input
4565  *balancedLinkControl, // input // Round Rock
4566  currChain, // output
4567  cumulativeRawChainRunTime, // output
4568  cumulativeRawChainRejections, // output
4569  &currLogLikelihoodValues, // output // likelihood is important
4570  &currLogTargetValues); // output);
4571 
4572  //***********************************************************
4573  // Perform checkpoint if necessary
4574  //***********************************************************
4575  stopAtEndOfLevel = currOptions->m_stopAtEnd;
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;
4581  }
4582  }
4583  if (performCheckpoint) {
4584  checkpointML(currExponent, // input
4585  currEta, // input
4586  currChain, // input
4587  currLogLikelihoodValues, // input
4588  currLogTargetValues); // input
4589  }
4590 
4591  //***********************************************************
4592  // Just free some memory
4593  //***********************************************************
4594  {
4595  delete unifiedCovMatrix;
4596 
4597  for (unsigned int i = 0; i < balancedLinkControl->balLinkedChains.size(); ++i) {
4598  queso_require_msg(balancedLinkControl->balLinkedChains[i].initialPosition, "Initial position pointer in step 9 should not be NULL");
4599  delete balancedLinkControl->balLinkedChains[i].initialPosition;
4600  balancedLinkControl->balLinkedChains[i].initialPosition = NULL;
4601  }
4602  balancedLinkControl->balLinkedChains.clear();
4603  }
4604 
4605  //***********************************************************
4606  // Step 11 of 11: filter chain if requested
4607  //***********************************************************
4608  m_currStep = 11;
4609  unsigned int unifiedNumberOfRejections = 0;
4610  if (m_env.inter0Rank() >= 0) {
4611  generateSequence_Step11_inter0(currOptions, // input
4612  currUnifiedRequestedNumSamples, // input
4613  cumulativeRawChainRejections, // input
4614  currChain, // input/output
4615  currLogLikelihoodValues, // input/output
4616  currLogTargetValues, // input/output
4617  unifiedNumberOfRejections); // output
4618  }
4619 
4620  minLogLike = -INFINITY;
4621  maxLogLike = INFINITY;
4622  currLogLikelihoodValues.subMinMaxExtra(0,
4623  currLogLikelihoodValues.subSequenceSize(),
4624  minLogLike,
4625  maxLogLike);
4626  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4627  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4628  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4629  << ", sub minLogLike = " << minLogLike
4630  << ", sub maxLogLike = " << maxLogLike
4631  << std::endl;
4632  }
4633 
4634  m_env.fullComm().Barrier();
4635 
4636  minLogLike = -INFINITY;
4637  maxLogLike = INFINITY;
4638  currLogLikelihoodValues.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4639  0,
4640  currLogLikelihoodValues.subSequenceSize(),
4641  minLogLike,
4642  maxLogLike);
4643  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4644  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4645  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4646  << ", unified minLogLike = " << minLogLike
4647  << ", unified maxLogLike = " << maxLogLike
4648  << std::endl;
4649  }
4650 
4651  //***********************************************************
4652  // Prepare to end current level
4653  //***********************************************************
4654  double levelRunTime = MiscGetEllapsedSeconds(&timevalLevel);
4655  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4656  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4657  << ": ending level " << m_currLevel+LEVEL_REF_ID
4658  << ", having generated " << currChain.subSequenceSize()
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
4668  << std::endl;
4669  }
4670 
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");
4676 
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");
4681 
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());
4687 
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");
4692 
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");
4697 
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());
4703 
4704  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4705  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4706  << ", level " << m_currLevel+LEVEL_REF_ID
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
4713  << std::endl;
4714  }
4715  }
4716 
4717  if (currExponent != 1.) delete currOptions;
4718 
4719  struct timeval timevalLevelEnd;
4720  iRC = 0;
4721  iRC = gettimeofday(&timevalLevelEnd, NULL);
4722 
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"
4731  << std::endl;
4732  }
4733  } // end of level while
4734 
4735 
4736  // m_env.worldRank(),
4737  // "MLSampling<P_V,P_M>::generateSequence()",
4738  // "exponent has not achieved value '1' even after exiting level loop");
4739 
4740  //***********************************************************
4741  // Compute information gain
4742  // ln( \pi(D|M) ) = E[ln( \pi(D|\theta,M) )] - E[ln( \pi(\theta|D,M) / \pi(\theta|M) )]
4743  //***********************************************************
4744  if (m_env.inter0Rank() >= 0) { // KAUST
4745  queso_require_equal_to_msg(m_currLevel, m_logEvidenceFactors.size(), "invalid m_currLevel at the exit of the level loop");
4746  m_logEvidence = 0.;
4747  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
4748  m_logEvidence += m_logEvidenceFactors[i];
4749  }
4750 
4751 #if 1 // prudenci-2012-07-06
4752  m_meanLogLikelihood = currLogLikelihoodValues.unifiedMeanPlain(m_vectorSpace.numOfProcsForStorage() == 1);
4753 #else
4754  m_meanLogLikelihood = currLogLikelihoodValues.unifiedMeanExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4755  0,
4756  currLogLikelihoodValues.subSequenceSize());
4757 #endif
4758 
4759  m_eig = m_meanLogLikelihood - m_logEvidence;
4760 
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
4767  << std::endl;
4768  }
4769  }
4770 
4771  m_env.subComm().Bcast((void *) &m_logEvidence, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm'
4772  "MLSampling<P_V,P_M>::generateSequence()",
4773  "failed MPI.Bcast() for m_logEvidence");
4774 
4775  m_env.subComm().Bcast((void *) &m_meanLogLikelihood, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm'
4776  "MLSampling<P_V,P_M>::generateSequence()",
4777  "failed MPI.Bcast() for m_meanLogLikelihood");
4778 
4779  m_env.subComm().Bcast((void *) &m_eig, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm'
4780  "MLSampling<P_V,P_M>::generateSequence()",
4781  "failed MPI.Bcast() for m_eig");
4782 
4783  //***********************************************************
4784  // Prepare to return
4785  //***********************************************************
4786  workingChain.clear();
4787  workingChain.resizeSequence(currChain.subSequenceSize());
4788  P_V auxVec(m_vectorSpace.zeroVector());
4789  for (unsigned int i = 0; i < workingChain.subSequenceSize(); ++i) {
4790  if (m_env.inter0Rank() >= 0) {
4791  currChain.getPositionValues(i,auxVec);
4792  }
4793  workingChain.setPositionValues(i,auxVec);
4794  }
4795 
4796  if (workingLogLikelihoodValues) *workingLogLikelihoodValues = currLogLikelihoodValues;
4797  if (workingLogTargetValues ) *workingLogTargetValues = currLogTargetValues;
4798 
4799  struct timeval timevalRoutineEnd;
4800  iRC = 0;
4801  iRC = gettimeofday(&timevalRoutineEnd, NULL);
4802 
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"
4810  << std::endl;
4811  }
4812 
4813  return;
4814 }
4815 
4816 template<class P_V,class P_M>
4817 std::ostream& operator<<(std::ostream& os, const MLSampling<P_V,P_M>& obj)
4818 {
4819  obj.print(os);
4820 
4821  return os;
4822 }
4823 
4824 template <class P_V,class P_M>
4825 double MLSampling<P_V,P_M>::logEvidence() const
4826 {
4827  return m_logEvidence;
4828 }
4829 
4830 template <class P_V,class P_M>
4831 double MLSampling<P_V,P_M>::meanLogLikelihood() const
4832 {
4833  return m_meanLogLikelihood;
4834 }
4835 
4836 template <class P_V,class P_M>
4837 double MLSampling<P_V,P_M>::eig() const
4838 {
4839  return m_eig;
4840 }
4841 
4842 } // End namespace QUESO
4843 
4844 template class QUESO::MLSampling<QUESO::GslVector, QUESO::GslMatrix>;
UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
QUESO::MLSampling::generateSequence_Step02_inter0
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)...
Definition: MLSampling.C:2486
QUESO::MLSamplingLevelOptions::m_stopAtEnd
bool m_stopAtEnd
Stop at end of such level.
Definition: MLSamplingLevelOptions.h:158
QUESO::SequenceOfVectors::unifiedReadContents
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
Definition: SequenceOfVectors.C:1763
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:81
QUESO::MLSampling::sampleIndexes_proc0
void sampleIndexes_proc0(unsigned int unifiedRequestedNumSamples, const std::vector< double > &unifiedWeightStdVectorAtProc0Only, std::vector< unsigned int > &unifiedIndexCountersAtProc0Only)
Definition: MLSampling.C:90
QUESO::MLSampling::generateSequence_Step06_all
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).
Definition: MLSampling.C:3096
RawValue_MPI_UNSIGNED
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:68
QUESO::MLSamplingLevelOptions::m_filteredChainDataOutputFileType
std::string m_filteredChainDataOutputFileType
Type of output file for filtered chain.
Definition: MLSamplingLevelOptions.h:281
QUESO::BaseVectorSequence::append
void append(const BaseVectorSequence< V, M > &src, unsigned int initialPos, unsigned int numPos)
Appends the vector src to this vector.
Definition: VectorSequence.C:354
QUESO::ScalarSequence::getUnifiedContentsAtProc0Only
void getUnifiedContentsAtProc0Only(bool useOnlyInter0Comm, std::vector< T > &outputVec) const
Gets the unified contents of processor of rank equals to 0.
Definition: ScalarSequence.C:235
QUESO::MLSampling::MLSampling
MLSampling(const char *prefix, const BaseVectorRV< P_V, P_M > &priorRv, const BaseScalarFunction< P_V, P_M > &likelihoodFunction)
Constructor.
Definition: MLSampling.C:3980
queso_error_msg
#define queso_error_msg(msg)
Definition: asserts.h:47
QUESO::MLSamplingLevelOptions::m_minEffectiveSizeRatio
double m_minEffectiveSizeRatio
Minimum allowed effective size ratio wrt previous level.
Definition: MLSamplingLevelOptions.h:179
QUESO::MLSamplingLevelOptions::m_initialPositionUsePreviousLevelLikelihood
bool m_initialPositionUsePreviousLevelLikelihood
Use previous level likelihood for initial chain position instead of re-computing it from target pdf...
Definition: MLSamplingLevelOptions.h:216
QUESO::MLSamplingLevelOptions::m_filteredChainLag
unsigned int m_filteredChainLag
Spacing for chain filtering.
Definition: MLSamplingLevelOptions.h:275
QUESO::GenericVectorRV::setPdf
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
Definition: GenericVectorRV.C:90
QUESO::RawType_MPI_Status
int RawType_MPI_Status
Definition: MpiComm.h:62
QUESO::BaseVectorSequence::computeFilterParams
void computeFilterParams(std::ofstream *passedOfs, unsigned int &initialPos, unsigned int &spacing)
Computes the filtering parameters spacing for the sequence of vectors.
Definition: VectorSequence.C:615
QUESO::SequenceOfVectors::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Definition: SequenceOfVectors.C:1493
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
QUESO::MLSampling::meanLogLikelihood
double meanLogLikelihood() const
Method to calculate the mean of the logarithm of the likelihood.
Definition: MLSampling.C:4831
QUESO::MLSamplingLevelOptions::m_initialValuesOfDisabledParameters
std::vector< double > m_initialValuesOfDisabledParameters
Definition: MLSamplingLevelOptions.h:221
QUESO::ScalarSequence::clear
void clear()
Clears the sequence of scalars.
Definition: ScalarSequence.C:123
QUESO::MLSampling::generateSequence_Level0_all
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.
Definition: MLSampling.C:2293
QUESO::MetropolisHastingsSG::getRawChainInfo
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
Definition: MetropolisHastingsSG.C:1430
queso_require_less_msg
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:87
QUESO::MLSamplingLevelOptions::m_loadBalanceTreshold
double m_loadBalanceTreshold
Perform load balancing if load unbalancing ratio &gt; threshold.
Definition: MLSamplingLevelOptions.h:176
QUESO::MLSamplingLevelOptions::m_filteredChainDataOutputFileName
std::string m_filteredChainDataOutputFileName
Name of output file for filtered chain.
Definition: MLSamplingLevelOptions.h:278
QUESO::MLSamplingLevelOptions::m_dataOutputAllowedSet
std::set< unsigned int > m_dataOutputAllowedSet
subEnvs that will write to generic output file.
Definition: MLSamplingLevelOptions.h:167
QUESO::MLSampling::justBalance_proc0
void justBalance_proc0(const MLSamplingLevelOptions *currOptions, std::vector< ExchangeInfoStruct > &exchangeStdVec)
Definition: MLSampling.C:1495
QUESO::MLSamplingLevelOptions::m_filteredChainDiscardedPortion
double m_filteredChainDiscardedPortion
Initial discarded portion for chain filtering.
Definition: MLSamplingLevelOptions.h:272
QUESO::MLSamplingLevelOptions::m_covRejectionRate
double m_covRejectionRate
c.o.v. for judging attempted rejection rate at current level.
Definition: MLSamplingLevelOptions.h:194
QUESO::MLSampling::restartML
void restartML(double &currExponent, double &currEta, SequenceOfVectors< P_V, P_M > &currChain, ScalarSequence< double > &currLogLikelihoodValues, ScalarSequence< double > &currLogTargetValues)
Restarts ML algorithm.
Definition: MLSampling.C:2113
QUESO::MLSampling::generateSequence_Step05_inter0
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).
Definition: MLSampling.C:3016
QUESO::MLSamplingLevelOptions::m_rawChainDataOutputFileType
std::string m_rawChainDataOutputFileType
Type of output file for raw chain.
Definition: MLSamplingLevelOptions.h:253
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2001
ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
#define ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
Definition: MLSampling.h:46
QUESO::SequenceOfVectors::setPositionValues
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
Definition: SequenceOfVectors.C:206
QUESO::MLSampling::generateSequence_Step08_all
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).
Definition: MLSampling.C:3197
QUESO::FiniteDistribution
A templated class for a finite distribution.
Definition: FiniteDistribution.h:45
QUESO::MLSampling::logEvidence
double logEvidence() const
Method to calculate the logarithm of the evidence.
Definition: MLSampling.C:4825
QUESO::MLSamplingLevelOptions::m_dataOutputFileName
std::string m_dataOutputFileName
Name of generic output file.
Definition: MLSamplingLevelOptions.h:161
QUESO::MLSampling::generateUnbLinkedChains_all
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)
Definition: MLSampling.C:886
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:282
QUESO::ExchangeInfoStruct::numberOfPositions
unsigned int numberOfPositions
Definition: MLSampling.h:71
QUESO::ScalarSequence::unifiedMaxPlain
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
Definition: ScalarSequence.C:382
queso_require_less_equal_msg
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:89
queso_require_msg
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
QUESO::SequenceOfVectors::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: SequenceOfVectors.C:2042
QUESO::MLSampling::eig
double eig() const
Calculates the expected information gain value, EIG.
Definition: MLSampling.C:4837
queso_require_greater_equal_msg
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
QUESO::ScalarSequence::unifiedMinMaxExtra
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...
Definition: ScalarSequence.C:1541
QUESO::MLSamplingLevelOptions::m_maxEffectiveSizeRatio
double m_maxEffectiveSizeRatio
Maximum allowed effective size ratio wrt previous level.
Definition: MLSamplingLevelOptions.h:182
QUESO::BaseVectorSequence::clear
void clear()
Reset the values and the size of the sequence of vectors.
Definition: VectorSequence.C:133
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:89
QUESO::MLSampling::m_env
const BaseEnvironment & m_env
Queso enviroment.
Definition: MLSampling.h:470
QUESO::MLSampling::checkpointML
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.
Definition: MLSampling.C:2006
QUESO::UnbalancedLinkedChainsPerNodeStruct::unbLinkedChains
std::vector< UnbalancedLinkedChainControlStruct > unbLinkedChains
Definition: MLSampling.h:101
RawValue_MPI_MIN
#define RawValue_MPI_MIN
Definition: MpiComm.h:69
QUESO::MLSamplingLevelOptions::m_drMaxNumExtraStages
unsigned int m_drMaxNumExtraStages
&#39;dr&#39; maximum number of extra stages.
Definition: MLSamplingLevelOptions.h:309
QUESO::MLSampling::~MLSampling
~MLSampling()
Destructor.
Definition: MLSampling.C:4013
QUESO::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2683
QUESO::MLSamplingLevelOptions::m_scaleCovMatrix
bool m_scaleCovMatrix
Whether or not scale proposal covariance matrix.
Definition: MLSamplingLevelOptions.h:185
QUESO::ScalarSequence::unifiedReadContents
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
Definition: ScalarSequence.C:2908
QUESO::SequenceOfVectors::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
Definition: SequenceOfVectors.C:103
QUESO::ScalarSequence::unifiedMeanExtra
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...
Definition: ScalarSequence.C:876
QUESO::MLSampling::generateSequence_Step07_inter0
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...
Definition: MLSampling.C:3129
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:72
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:278
QUESO::InstantiateIntersection
VectorSet< V, M > * InstantiateIntersection(const VectorSet< V, M > &domain1, const VectorSet< V, M > &domain2)
This method calculates the intersection of domain1 and domain2.
Definition: InstantiateIntersection.C:37
QUESO::MLSampling::generateSequence_Step10_all
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).
Definition: MLSampling.C:3700
QUESO::ScalarSequence::append
void append(const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
Appends the scalar sequence src to this sequence.
Definition: ScalarSequence.C:2499
QUESO::SequenceOfVectors::getPositionValues
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
Definition: SequenceOfVectors.C:188
QUESO::MLSampling::generateSequence_Step03_inter0
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 &#39;m_logEvidenceFactors&#39; (St...
Definition: MLSampling.C:2618
QUESO::MLSampling::decideOnBalancedChains_all
bool decideOnBalancedChains_all(const MLSamplingLevelOptions *currOptions, unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, std::vector< ExchangeInfoStruct > &exchangeStdVec)
Definition: MLSampling.C:144
QUESO::MLSampling::prepareBalLinkedChains_inter0
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)
Definition: MLSampling.C:325
QUESO::FiniteDistribution::sample
unsigned int sample() const
Samples.
Definition: FiniteDistribution.C:140
QUESO::ScalarSequence::setName
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
Definition: ScalarSequence.C:115
QUESO::MiscCheckForSameValueInAllNodes
bool MiscCheckForSameValueInAllNodes(T &inputValue, double acceptableTreshold, const MpiComm &comm, const char *whereString)
Definition: Miscellaneous.C:479
QUESO::MLSampling::generateSequence_Step09_all
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 &lt;= rejection rate &lt;= max (Step 09 from ML algorithm)...
Definition: MLSampling.C:3230
RawValue_MPI_MAX
#define RawValue_MPI_MAX
Definition: MpiComm.h:70
QUESO::MLSamplingLevelOptions::m_parameterDisabledSet
std::set< unsigned int > m_parameterDisabledSet
Definition: MLSamplingLevelOptions.h:219
QUESO::MLSampling::generateSequence_Step01_inter0
void generateSequence_Step01_inter0(const MLSamplingLevelOptions *currOptions, unsigned int &unifiedRequestedNumSamples)
Reads options for the ML algorithm (Step 01 from ML algorithm).
Definition: MLSampling.C:2440
QUESO::MetropolisHastingsSG
A templated class that represents a Metropolis-Hastings generator of samples.
Definition: MetropolisHastingsSG.h:123
reason
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)
QUESO::ScalarSequence::unifiedMeanPlain
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Definition: ScalarSequence.C:407
QUESO::MLSampling::generateSequence_Step11_inter0
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).
Definition: MLSampling.C:3814
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
QUESO::MLSampling::generateSequence
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
Definition: MLSampling.C:4025
LEVEL_REF_ID
#define LEVEL_REF_ID
Definition: MLSamplingLevelOptions.h:28
QUESO::MHRawChainInfoStruct
A struct that represents a Metropolis-Hastings sample.
Definition: MetropolisHastingsSG.h:58
QUESO::ScalarSequence::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:136
QUESO::ScalarSequence::subMinMaxExtra
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...
Definition: ScalarSequence.C:1512
QUESO::MLSamplingLevelOptions::m_filteredChainGenerate
bool m_filteredChainGenerate
Whether or not to generate filtered chain.
Definition: MLSamplingLevelOptions.h:269
QUESO::MLSamplingLevelOptions::scanOptionsValues
void scanOptionsValues(const MLSamplingLevelOptions *defaultOptions)
Definition: MLSamplingLevelOptions.C:438
QUESO::ExchangeInfoStruct::originalIndexOfInitialPosition
unsigned int originalIndexOfInitialPosition
Definition: MLSampling.h:69
QUESO::BaseVectorSequence::resizeSequence
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2431
QUESO::MetropolisHastingsSG::generateSequence
void generateSequence(BaseVectorSequence< P_V, P_M > &workingChain, ScalarSequence< double > *workingLogLikelihoodValues, ScalarSequence< double > *workingLogTargetValues)
Method to generate the chain.
Definition: MetropolisHastingsSG.C:977
QUESO::MLSamplingLevelOptions
This class provides options for each level of the Multilevel sequence generator if no input file is a...
Definition: MLSamplingLevelOptions.h:117
QUESO::MLSampling::generateSequence_Step04_inter0
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).
Definition: MLSampling.C:2922
QUESO::MLSampling::generateBalLinkedChains_all
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)
Definition: MLSampling.C:628
QUESO::MLSamplingLevelOptions::m_minRejectionRate
double m_minRejectionRate
minimum allowed attempted rejection rate at current level
Definition: MLSamplingLevelOptions.h:188
QUESO::MLSamplingLevelOptions::m_totallyMute
bool m_totallyMute
Whether or not to be totally mute (no printout message).
Definition: MLSamplingLevelOptions.h:201
QUESO::MLSampling::prepareUnbLinkedChains_inter0
void prepareUnbLinkedChains_inter0(unsigned int indexOfFirstWeight, unsigned int indexOfLastWeight, const std::vector< unsigned int > &unifiedIndexCountersAtProc0Only, UnbalancedLinkedChainsPerNodeStruct &unbalancedLinkControl)
Definition: MLSampling.C:465
QUESO::MLSamplingLevelOptions::m_rawChainSize
unsigned int m_rawChainSize
Size of raw chain.
Definition: MLSamplingLevelOptions.h:232
QUESO::ScalarSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
Definition: ScalarSequence.C:143
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
QUESO::MLSamplingLevelOptions::m_minAcceptableEta
double m_minAcceptableEta
minimum acceptable eta,
Definition: MLSamplingLevelOptions.h:198
QUESO::MLSamplingLevelOptions::m_loadBalanceAlgorithmId
unsigned int m_loadBalanceAlgorithmId
Perform load balancing with chosen algorithm (0 = no balancing).
Definition: MLSamplingLevelOptions.h:173
QUESO::MLSamplingLevelOptions::m_rawChainDataOutputFileName
std::string m_rawChainDataOutputFileName
Name of output file for raw chain.
Definition: MLSamplingLevelOptions.h:247
UQ_MH_SG_FILENAME_FOR_NO_FILE
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
Definition: MetropolisHastingsSGOptions.h:36
RawValue_MPI_CHAR
#define RawValue_MPI_CHAR
Definition: MpiComm.h:65
QUESO::MLSamplingLevelOptions::m_maxRejectionRate
double m_maxRejectionRate
maximum allowed attempted rejection rate at current level.
Definition: MLSamplingLevelOptions.h:191
QUESO::BalancedLinkedChainsPerNodeStruct::balLinkedChains
std::vector< BalancedLinkedChainControlStruct< P_V > > balLinkedChains
Definition: MLSampling.h:88
QUESO::BaseVectorSequence::setName
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
Definition: VectorSequence.C:125
RawValue_MPI_DOUBLE
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
QUESO::BaseVectorSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize() const
Calculates the size of the unified sequence of vectors.
Definition: VectorSequence.C:70
QUESO::BaseVectorSequence::setPositionValues
virtual void setPositionValues(unsigned int posId, const V &vec)=0
Set the values in vec at position posId of the sequence. See template specialization.
QUESO::SequenceOfVectors::resizeSequence
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
Definition: SequenceOfVectors.C:110
QUESO::MLSampling::mpiExchangePositions_inter0
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)
Definition: MLSampling.C:1783
QUESO::MLSamplingLevelOptions::m_amAdaptInterval
unsigned int m_amAdaptInterval
&#39;am&#39; adaptation interval.
Definition: MLSamplingLevelOptions.h:325
QUESO::MLSampling
A templated class that represents a Multilevel generator of samples.
Definition: MLSampling.h:122
Generated on Fri Jun 17 2016 14:17:42 for queso-0.55.0 by  doxygen 1.8.5