queso-0.53.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().Gather((void *) &auxUInt, 1, RawValue_MPI_UNSIGNED, (void *) &allFirstIndexes[0], (int) 1, RawValue_MPI_UNSIGNED, 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().Gather((void *) &auxUInt, 1, RawValue_MPI_UNSIGNED, (void *) &allLastIndexes[0], (int) 1, RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &minRatio, (int) auxBuf.size(), RawValue_MPI_DOUBLE, 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().Allreduce((void *) &auxBuf[0], (void *) &maxRatio, (int) auxBuf.size(), RawValue_MPI_DOUBLE, 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().Allreduce((void *) &auxBuf[0], (void *) &minModifiedSubNumSamples, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &maxModifiedSubNumSamples, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &sumModifiedSubNumSamples, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &minNumberOfPositions, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &maxNumberOfPositions, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &sumNumberOfPositions, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &minNumberOfPositions, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &maxNumberOfPositions, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Allreduce((void *) &auxBuf[0], (void *) &sumNumberOfPositions, (int) auxBuf.size(), RawValue_MPI_UNSIGNED, 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().Gatherv((void *) &sendbuf[0], (int) sendcnt, RawValue_MPI_DOUBLE, (void *) &recvbuf[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, r, // LOAD BALANCE
1944  "MLSampling<P_V,P_M>::mpiExchangePositions_inter0()",
1945  "failed MPI.Gatherv()");
1946 #endif
1947 
1949  // Make sanity checks
1951 
1953  // Transfer data from 'recvbuf' to 'balancedLinkControl'
1954  // Remember that finalNumChainsPerNode[r] = (totalNumberOfInitialPositionsNodeRHasToReceive + numberOfInitialPositionsNodeRAlreadyHas)
1955  // Remember that totalNumberOfInitialPositionsNodeRHasToReceive = totalNumberOfChainLenghtsNodeRHasToInherit
1957  if (m_env.inter0Rank() == (int) r) {
1958  balancedLinkControl.balLinkedChains.resize(finalNumChainsPerNode[r]);
1959  unsigned int auxIndex = 0;
1960 
1961  for (unsigned int i = 0; i < Nc; ++i) {
1962  if ((exchangeStdVec[i].finalNodeOfInitialPosition == (int) r) &&
1963  (exchangeStdVec[i].originalNodeOfInitialPosition == (int) r)) {
1964  unsigned int originalIndex = exchangeStdVec[i].originalIndexOfInitialPosition;
1965  prevChain.getPositionValues(originalIndex, auxInitialPosition);
1966  balancedLinkControl.balLinkedChains[auxIndex].initialPosition = new P_V(auxInitialPosition);
1967  balancedLinkControl.balLinkedChains[auxIndex].initialLogPrior = prevLogTargetValues[originalIndex] - prevLogLikelihoodValues[originalIndex];
1968  balancedLinkControl.balLinkedChains[auxIndex].initialLogLikelihood = expRatio*prevLogLikelihoodValues[originalIndex];
1969  balancedLinkControl.balLinkedChains[auxIndex].numberOfPositions = exchangeStdVec[i].numberOfPositions;
1970  auxIndex++;
1971  }
1972  }
1973 
1974  for (unsigned int i = 0; i < totalNumberOfInitialPositionsNodeRHasToReceive; ++i) {
1975  for (unsigned int j = 0; j < dimSize; ++j) {
1976  auxInitialPosition[j] = recvbuf[i*nValuesPerInitialPosition + j];
1977  }
1978  balancedLinkControl.balLinkedChains[auxIndex].initialPosition = new P_V(auxInitialPosition);
1979  double prevLogLikelihood = recvbuf[i*nValuesPerInitialPosition + dimSize];
1980  double prevLogTarget = recvbuf[i*nValuesPerInitialPosition + dimSize + 1];
1981  balancedLinkControl.balLinkedChains[auxIndex].initialLogPrior = prevLogTarget - prevLogLikelihood;
1982  balancedLinkControl.balLinkedChains[auxIndex].initialLogLikelihood = expRatio*prevLogLikelihood;
1983  balancedLinkControl.balLinkedChains[auxIndex].numberOfPositions = chainLenghtsNodeRHasToInherit[i]; // aqui 3
1984  auxIndex++;
1985  }
1986  }
1987 
1988  m_env.inter0Comm().Barrier();
1989  } // for 'r'
1990 
1991  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1992  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::mpiExchangePositions_inter0()"
1993  << ", level " << m_currLevel+LEVEL_REF_ID
1994  << ", step " << m_currStep
1995  << std::endl;
1996  }
1997 
1998  return;
1999 }
2000 
2001 // Statistical/private methods-----------------------
2002 template <class P_V,class P_M>
2003 void
2004 MLSampling<P_V,P_M>::checkpointML(
2005  double currExponent, // input
2006  double currEta, // input
2007  const SequenceOfVectors<P_V,P_M>& currChain, // input
2008  const ScalarSequence<double>& currLogLikelihoodValues, // input
2009  const ScalarSequence<double>& currLogTargetValues) // input
2010 {
2011  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2012  *m_env.subDisplayFile() << "\n CHECKPOINTING initiating at level " << m_currLevel
2013  << "\n" << std::endl;
2014  }
2015 
2016  //******************************************************************************
2017  // Write 'control' file without 'level' spefication in name
2018  //******************************************************************************
2019  unsigned int quantity1 = currChain.unifiedSequenceSize();
2020  unsigned int quantity2 = currLogLikelihoodValues.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2021  unsigned int quantity3 = currLogTargetValues.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2022  if (m_env.inter0Rank() >= 0) {
2023  queso_require_equal_to_msg(m_logEvidenceFactors.size(), m_currLevel, "number of evidence factors is not consistent");
2024  queso_require_equal_to_msg(quantity1, quantity2, "quantity2 is not consistent");
2025  queso_require_equal_to_msg(quantity1, quantity3, "quantity3 is not consistent");
2026  }
2027 
2028  if (m_env.fullRank() == 0) {
2029  std::ofstream* ofsVar = new std::ofstream((m_options.m_restartOutput_baseNameForFiles + "Control.txt").c_str(),
2030  std::ofstream::out | std::ofstream::trunc);
2031  *ofsVar << m_currLevel << std::endl // 1
2032  << m_vectorSpace.dimGlobal() << std::endl // 2
2033  << currExponent << std::endl // 3
2034  << currEta << std::endl // 4
2035  << quantity1 << std::endl; // 5
2036  unsigned int savedPrecision = ofsVar->precision();
2037  ofsVar->precision(16);
2038  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2039  *ofsVar << m_logEvidenceFactors[i] << std::endl;
2040  }
2041  ofsVar->precision(savedPrecision);
2042  *ofsVar << "COMPLETE" << std::endl; // 6 = ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
2043 
2044  delete ofsVar;
2045  }
2046  m_env.fullComm().Barrier();
2047 
2048  //******************************************************************************
2049  // Write three 'data' files
2050  //******************************************************************************
2051  char levelSufix[256];
2052  sprintf(levelSufix,"%d",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
2053 
2054  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2055  *m_env.subDisplayFile() << "\n CHECKPOINTING chain at level " << m_currLevel
2056  << "\n" << std::endl;
2057  }
2058  currChain.unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles + "Chain_l" + levelSufix,
2059  m_options.m_restartOutput_fileType);
2060  m_env.fullComm().Barrier();
2061 
2062  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2063  *m_env.subDisplayFile() << "\n CHECKPOINTING like at level " << m_currLevel
2064  << "\n" << std::endl;
2065  }
2066  currLogLikelihoodValues.unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles + "LogLike_l" + levelSufix,
2067  m_options.m_restartOutput_fileType);
2068  m_env.fullComm().Barrier();
2069 
2070  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2071  *m_env.subDisplayFile() << "\n CHECKPOINTING target at level " << m_currLevel
2072  << "\n" << std::endl;
2073  }
2074  currLogTargetValues.unifiedWriteContents(m_options.m_restartOutput_baseNameForFiles + "LogTarget_l" + levelSufix,
2075  m_options.m_restartOutput_fileType);
2076  m_env.fullComm().Barrier();
2077 
2078  //******************************************************************************
2079  // Write 'control' file *with* 'level' spefication in name
2080  //******************************************************************************
2081  if (m_env.fullRank() == 0) {
2082  std::ofstream* ofsVar = new std::ofstream((m_options.m_restartOutput_baseNameForFiles + "Control_l" + levelSufix + ".txt").c_str(),
2083  std::ofstream::out | std::ofstream::trunc);
2084  *ofsVar << m_currLevel << std::endl // 1
2085  << m_vectorSpace.dimGlobal() << std::endl // 2
2086  << currExponent << std::endl // 3
2087  << currEta << std::endl // 4
2088  << quantity1 << std::endl; // 5
2089  unsigned int savedPrecision = ofsVar->precision();
2090  ofsVar->precision(16);
2091  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2092  *ofsVar << m_logEvidenceFactors[i] << std::endl;
2093  }
2094  ofsVar->precision(savedPrecision);
2095  *ofsVar << "COMPLETE" << std::endl; // 6 = ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
2096 
2097  delete ofsVar;
2098  }
2099  m_env.fullComm().Barrier();
2100 
2101  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2102  *m_env.subDisplayFile() << "\n CHECKPOINTING done at level " << m_currLevel
2103  << "\n" << std::endl;
2104  }
2105 
2106  return;
2107 }
2108 //---------------------------------------------------
2109 template <class P_V,class P_M>
2110 void
2111 MLSampling<P_V,P_M>::restartML(
2112  double& currExponent, // output
2113  double& currEta, // output
2114  SequenceOfVectors<P_V,P_M>& currChain, // output
2115  ScalarSequence<double>& currLogLikelihoodValues, // output
2116  ScalarSequence<double>& currLogTargetValues) // output
2117 {
2118  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2119  *m_env.subDisplayFile() << "\n RESTARTING initiating at level " << m_currLevel
2120  << "\n" << std::endl;
2121  }
2122 
2123  //******************************************************************************
2124  // Read 'control' file
2125  //******************************************************************************
2126  unsigned int vectorSpaceDim = 0;
2127  unsigned int quantity1 = 0;
2128  std::string checkingString("");
2129  if (m_env.fullRank() == 0) {
2130  std::ifstream* ifsVar = new std::ifstream((m_options.m_restartInput_baseNameForFiles + "Control.txt").c_str(),
2131  std::ifstream::in);
2132 
2133  //******************************************************************************
2134  // Determine number of lines
2135  //******************************************************************************
2136  unsigned int numLines = std::count(std::istreambuf_iterator<char>(*ifsVar),
2137  std::istreambuf_iterator<char>(),
2138  '\n');
2139  ifsVar->seekg(0,std::ios_base::beg);
2140  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2141  *m_env.subDisplayFile() << "Restart input file has " << numLines
2142  << " lines"
2143  << std::endl;
2144  }
2145 
2146  //******************************************************************************
2147  // Read all values
2148  //******************************************************************************
2149  *ifsVar >> m_currLevel; // 1
2150  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");
2151 
2152  m_logEvidenceFactors.clear();
2153  m_logEvidenceFactors.resize(m_currLevel,0.);
2154  *ifsVar >> vectorSpaceDim // 2
2155  >> currExponent // 3
2156  >> currEta // 4
2157  >> quantity1; // 5
2158  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2159  *ifsVar >> m_logEvidenceFactors[i];
2160  }
2161  *ifsVar >> checkingString; // 6 = ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
2162  queso_require_equal_to_msg(checkingString, "COMPLETE", "control txt input file is not complete");
2163 
2164  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2165  *m_env.subDisplayFile() << "Restart input file has the following information:"
2166  << "\n m_currLevel = " << m_currLevel
2167  << "\n vectorSpaceDim = " << vectorSpaceDim
2168  << "\n currExponent = " << currExponent
2169  << "\n currEta = " << currEta
2170  << "\n quantity1 = " << quantity1;
2171  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2172  *m_env.subDisplayFile() << "\n [" << i << "] = " << m_logEvidenceFactors[i];
2173  }
2174  *m_env.subDisplayFile() << std::endl;
2175  }
2176 
2177 #if 0 // For debug only
2178  std::string tmpString;
2179  for (unsigned int i = 0; i < 2; ++i) {
2180  *ifsVar >> tmpString;
2181  std::cout << "Just read '" << tmpString << "'" << std::endl;
2182  }
2183  while ((lineId < numLines) && (ifsVar->eof() == false)) {
2184  }
2185  ifsVar->ignore(maxCharsPerLine,'\n');
2186 #endif
2187 
2188  delete ifsVar;
2189  } // if (m_env.fullRank() == 0)
2190  m_env.fullComm().Barrier();
2191 
2192  //******************************************************************************
2193  // MPI_Bcast 'm_currLevel'
2194  //******************************************************************************
2195  unsigned int tmpUint = (unsigned int) m_currLevel;
2196  m_env.fullComm().Bcast((void *) &tmpUint, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'fullComm'
2197  "MLSampling<P_V,P_M>::restartML()",
2198  "failed MPI.Bcast() for m_currLevel");
2199  if (m_env.fullRank() != 0) {
2200  m_currLevel = tmpUint;
2201  }
2202 
2203  //******************************************************************************
2204  // MPI_Bcast the rest of the information just read
2205  //******************************************************************************
2206  std::vector<double> tmpData(ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA-1+m_currLevel,0.);
2207  if (m_env.fullRank() == 0) {
2208  tmpData[0] = vectorSpaceDim;
2209  tmpData[1] = currExponent;
2210  tmpData[2] = currEta;
2211  tmpData[3] = quantity1;
2212  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2213  tmpData[4+i] = m_logEvidenceFactors[i];
2214  }
2215  }
2216  else {
2217  m_logEvidenceFactors.clear();
2218  m_logEvidenceFactors.resize(m_currLevel,0.);
2219  }
2220  m_env.fullComm().Bcast((void *) &tmpData[0], (int) tmpData.size(), RawValue_MPI_DOUBLE, 0, // Yes, 'fullComm'
2221  "MLSampling<P_V,P_M>::restartML()",
2222  "failed MPI.Bcast() for rest of information read from input file");
2223  if (m_env.fullRank() != 0) {
2224  vectorSpaceDim = tmpData[0];
2225  currExponent = tmpData[1];
2226  currEta = tmpData[2];
2227  quantity1 = tmpData[3];
2228  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
2229  m_logEvidenceFactors[i] = tmpData[4+i];
2230  }
2231  }
2232 
2233  //******************************************************************************
2234  // Process read data in all MPI nodes now
2235  //******************************************************************************
2236  queso_require_equal_to_msg(vectorSpaceDim, m_vectorSpace.dimGlobal(), "read vector space dimension is not consistent");
2237  queso_require_msg(!((currExponent < 0.) || (currExponent > 1.)), "read currExponent is not consistent");
2238  queso_require_equal_to_msg((quantity1 % m_env.numSubEnvironments()), 0, "read size of chain should be a multiple of the number of subenvironments");
2239  unsigned int subSequenceSize = 0;
2240  subSequenceSize = ((double) quantity1) / ((double) m_env.numSubEnvironments());
2241 
2242  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2243  *m_env.subDisplayFile() << "Restart input file has the following information"
2244  << ": subSequenceSize = " << subSequenceSize
2245  << std::endl;
2246  }
2247 
2248  //******************************************************************************
2249  // Read three 'data' files
2250  //******************************************************************************
2251  char levelSufix[256];
2252  sprintf(levelSufix,"%d",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
2253 
2254  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2255  *m_env.subDisplayFile() << "\n RESTARTING chain at level " << m_currLevel
2256  << "\n" << std::endl;
2257  }
2258  currChain.unifiedReadContents(m_options.m_restartInput_baseNameForFiles + "Chain_l" + levelSufix,
2259  m_options.m_restartInput_fileType,
2260  subSequenceSize);
2261  m_env.fullComm().Barrier();
2262 
2263  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2264  *m_env.subDisplayFile() << "\n RESTARTING like at level " << m_currLevel
2265  << "\n" << std::endl;
2266  }
2267  currLogLikelihoodValues.unifiedReadContents(m_options.m_restartInput_baseNameForFiles + "LogLike_l" + levelSufix,
2268  m_options.m_restartInput_fileType,
2269  subSequenceSize);
2270  m_env.fullComm().Barrier();
2271 
2272  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2273  *m_env.subDisplayFile() << "\n RESTARTING target at level " << m_currLevel
2274  << "\n" << std::endl;
2275  }
2276  currLogTargetValues.unifiedReadContents(m_options.m_restartInput_baseNameForFiles + "LogTarget_l" + levelSufix,
2277  m_options.m_restartInput_fileType,
2278  subSequenceSize);
2279  m_env.fullComm().Barrier();
2280 
2281  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2282  *m_env.subDisplayFile() << "\n RESTARTING done at level " << m_currLevel
2283  << "\n" << std::endl;
2284  }
2285 
2286  return;
2287 }
2288 //---------------------------------------------------
2289 template <class P_V,class P_M>
2290 void
2291 MLSampling<P_V,P_M>::generateSequence_Level0_all(
2292  const MLSamplingLevelOptions& currOptions, // input
2293  unsigned int& unifiedRequestedNumSamples, // output
2294  SequenceOfVectors<P_V,P_M>& currChain, // output
2295  ScalarSequence<double>& currLogLikelihoodValues, // output
2296  ScalarSequence<double>& currLogTargetValues) // output
2297 {
2298  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2299  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateSequence()"
2300  << ": beginning level " << m_currLevel+LEVEL_REF_ID
2301  << ", currOptions.m_rawChainSize = " << currOptions.m_rawChainSize // Ok to use rawChainSize
2302  << std::endl;
2303  }
2304 
2305  int iRC = UQ_OK_RC;
2306  struct timeval timevalLevel;
2307  iRC = gettimeofday(&timevalLevel, NULL);
2308  if (iRC) {}; // just to remove compiler warning
2309 
2310  if (m_env.inter0Rank() >= 0) {
2311  unsigned int tmpSize = currOptions.m_rawChainSize;
2312  m_env.inter0Comm().Allreduce((void *) &tmpSize, (void *) &unifiedRequestedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
2313  "MLSampling<P_V,P_M>::generateSequence()",
2314  "failed MPI.Allreduce() for requested num samples in level 0");
2315  }
2316  else {
2317  unifiedRequestedNumSamples = currOptions.m_rawChainSize;
2318  }
2319 
2320  currChain.setName (currOptions.m_prefix + "rawChain" );
2321  currLogLikelihoodValues.setName(currOptions.m_prefix + "rawLogLikelihood");
2322  currLogTargetValues.setName (currOptions.m_prefix + "rawLogTarget" );
2323 
2324  currChain.resizeSequence (currOptions.m_rawChainSize); // Ok to use rawChainSize
2325  currLogLikelihoodValues.resizeSequence(currOptions.m_rawChainSize); // Ok to use rawChainSize
2326  currLogTargetValues.resizeSequence (currOptions.m_rawChainSize); // Ok to use rawChainSize
2327 
2328  P_V auxVec(m_vectorSpace.zeroVector());
2329  ScalarFunctionSynchronizer<P_V,P_M> likelihoodSynchronizer(m_likelihoodFunction,auxVec); // prudencio 2010-08-01
2330  for (unsigned int i = 0; i < currChain.subSequenceSize(); ++i) {
2331  //std::cout << "In QUESO: before prior realizer with i = " << i << std::endl;
2332  bool outOfSupport = true;
2333  do {
2334 
2335  m_priorRv.realizer().realization(auxVec); // gpmsa2
2336  if (m_numDisabledParameters > 0) { // gpmsa2
2337  unsigned int disabledCounter = 0;
2338  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2339  if (m_parameterEnabledStatus[paramId] == false) {
2340  auxVec[paramId] = currOptions.m_initialValuesOfDisabledParameters[disabledCounter];
2341  disabledCounter++;
2342  }
2343  }
2344  }
2345  auxVec.mpiBcast(0, m_env.subComm()); // prudencio 2010-08-01
2346 
2347  outOfSupport = !(m_targetDomain->contains(auxVec));
2348  } while (outOfSupport); // prudenci 2011-Oct-04
2349 
2350  currChain.setPositionValues(i,auxVec);
2351  // KAUST: all nodes should call likelihood
2352 #if 1 // prudencio 2010-08-01
2353  currLogLikelihoodValues[i] = likelihoodSynchronizer.callFunction(&auxVec,NULL,NULL,NULL,NULL,NULL,NULL); // likelihood is important
2354 #else
2355  currLogLikelihoodValues[i] = m_likelihoodFunction.lnValue(auxVec,NULL,NULL,NULL,NULL); // likelihood is important
2356 #endif
2357  currLogTargetValues[i] = m_priorRv.pdf().lnValue(auxVec,NULL,NULL,NULL,NULL) + currLogLikelihoodValues[i];
2358  //std::cout << "In QUESO: currLogTargetValues[" << i << "] = " << currLogTargetValues[i] << std::endl;
2359  }
2360 
2361  if (m_env.inter0Rank() >= 0) { // KAUST
2362 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2363  if (currOptions.m_rawChainComputeStats) {
2364  FilePtrSetStruct filePtrSet;
2365  m_env.openOutputFile(currOptions.m_dataOutputFileName,
2366  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
2367  currOptions.m_dataOutputAllowedSet,
2368  false,
2369  filePtrSet);
2370 
2371  //m_env.syncPrintDebugMsg("At level 0, calling computeStatistics for chain",1,10,m_env.inter0Comm()); // output debug
2372  currChain.computeStatistics(*currOptions.m_rawChainStatisticalOptionsObj,
2373  filePtrSet.ofsVar);
2374 
2375  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
2376  }
2377  // Compute MLE and MAP
2378  // rr0
2379 #endif
2380  if (currOptions.m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
2381  currChain.unifiedWriteContents (currOptions.m_rawChainDataOutputFileName,
2382  currOptions.m_rawChainDataOutputFileType);
2383  currLogLikelihoodValues.unifiedWriteContents(currOptions.m_rawChainDataOutputFileName,
2384  currOptions.m_rawChainDataOutputFileType);
2385  currLogTargetValues.unifiedWriteContents (currOptions.m_rawChainDataOutputFileName,
2386  currOptions.m_rawChainDataOutputFileType);
2387  }
2388 
2389  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2390  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2391  << ", level " << m_currLevel+LEVEL_REF_ID
2392  << ": finished generating " << currChain.subSequenceSize()
2393  << " chain positions"
2394  << std::endl;
2395 
2396  //unsigned int numZeros = 0;
2397  //for (unsigned int i = 0; i < currTargetValues.subSequenceSize(); ++i) {
2398  // *m_env.subDisplayFile() << "currTargetValues[" << i
2399  // << "] = " << currTargetValues[i]
2400  // << std::endl;
2401  // if (currTargetValues[i] == 0.) numZeros++;
2402  //}
2403  //*m_env.subDisplayFile() << "Number of zeros in currTargetValues = " << numZeros
2404  // << std::endl;
2405  }
2406 
2407  if (currOptions.m_filteredChainGenerate) {
2408  // todo
2409 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
2410  if (currOptions.m_filteredChainComputeStats) {
2411  // todo
2412 
2413  //currChain.computeStatistics(*currOptions.m_filteredChainStatisticalOptionsObj,
2414  // filePtrSet.ofsVar);
2415  }
2416  // Compute MLE and MAP
2417  // rr0
2418 #endif
2419  }
2420 
2421  } // KAUST
2422 
2423  queso_require_equal_to_msg(currChain.subSequenceSize(), currOptions.m_rawChainSize, "currChain (first one) has been generated with invalid size");
2424 
2425  double levelRunTime = MiscGetEllapsedSeconds(&timevalLevel);
2426  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2427  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2428  << ": ending level " << m_currLevel+LEVEL_REF_ID
2429  << ", total level time = " << levelRunTime << " seconds"
2430  << std::endl;
2431  }
2432 
2433  return;
2434 }
2435 //---------------------------------------------------
2436 template <class P_V,class P_M>
2437 void
2438 MLSampling<P_V,P_M>::generateSequence_Step01_inter0(
2439  const MLSamplingLevelOptions* currOptions, // input
2440  unsigned int& unifiedRequestedNumSamples) // output
2441 {
2442  int iRC = UQ_OK_RC;
2443  struct timeval timevalStep;
2444  iRC = gettimeofday(&timevalStep, NULL);
2445  if (iRC) {}; // just to remove compiler warning
2446 
2447  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2448  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2449  << ", level " << m_currLevel+LEVEL_REF_ID
2450  << ", step " << m_currStep
2451  << ": beginning step 1 of 11"
2452  << std::endl;
2453  }
2454 
2455  unsigned int tmpSize = currOptions->m_rawChainSize;
2456  // This computed 'unifiedRequestedNumSamples' needs to be recomputed only at the last
2457  // level, when 'currOptions' is replaced by 'lastLevelOptions' (see step 3 of 11)
2458  m_env.inter0Comm().Allreduce((void *) &tmpSize, (void *) &unifiedRequestedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
2459  "MLSampling<P_V,P_M>::generateSequence()",
2460  "failed MPI.Allreduce() for requested num samples in step 1");
2461 
2462  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2463  *m_env.subDisplayFile() << "KEY In MLSampling<P_V,P_M>::generateSequence()"
2464  << ", level " << m_currLevel+LEVEL_REF_ID
2465  << ", step " << m_currStep
2466  << ", currOptions->m_rawChainSize = " << currOptions->m_rawChainSize // Ok to use rawChainSize
2467  << std::endl;
2468  }
2469 
2470  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
2471  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2472  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2473  << ", level " << m_currLevel+LEVEL_REF_ID
2474  << ", step " << m_currStep
2475  << ", after " << stepRunTime << " seconds"
2476  << std::endl;
2477  }
2478 
2479  return;
2480 }
2481 //---------------------------------------------------
2482 template <class P_V,class P_M>
2483 void
2484 MLSampling<P_V,P_M>::generateSequence_Step02_inter0(
2485  const MLSamplingLevelOptions* currOptions, // input
2486  SequenceOfVectors<P_V,P_M>& currChain, // input/output
2487  ScalarSequence<double>& currLogLikelihoodValues, // input/output
2488  ScalarSequence<double>& currLogTargetValues, // input/output
2489  SequenceOfVectors<P_V,P_M>& prevChain, // output
2490  ScalarSequence<double>& prevLogLikelihoodValues, // output
2491  ScalarSequence<double>& prevLogTargetValues, // output
2492  unsigned int& indexOfFirstWeight, // output
2493  unsigned int& indexOfLastWeight) // output
2494 {
2495  int iRC = UQ_OK_RC;
2496  struct timeval timevalStep;
2497  iRC = gettimeofday(&timevalStep, NULL);
2498  if (iRC) {}; // just to remove compiler warning
2499 
2500  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2501  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2502  << ", level " << m_currLevel+LEVEL_REF_ID
2503  << ", step " << m_currStep
2504  << ": beginning step 2 of 11"
2505  << std::endl;
2506  }
2507 
2508  prevChain = currChain;
2509  currChain.clear();
2510  currChain.setName(currOptions->m_prefix + "rawChain");
2511 
2512  prevLogLikelihoodValues = currLogLikelihoodValues; // likelihood is important
2513  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2514  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
2515  << ", level " << m_currLevel+LEVEL_REF_ID
2516  << ", step " << m_currStep
2517  << ", prevLogLikelihoodValues[0] = " << prevLogLikelihoodValues[0]
2518  << std::endl;
2519  }
2520  prevLogTargetValues = currLogTargetValues;
2521 
2522  currLogLikelihoodValues.clear();
2523  currLogLikelihoodValues.setName(currOptions->m_prefix + "rawLogLikelihood");
2524 
2525  currLogTargetValues.clear();
2526  currLogTargetValues.setName(currOptions->m_prefix + "rawLogTarget");
2527 
2528 #if 0 // For debug only
2529  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2530  P_V prevPosition(m_vectorSpace.zeroVector());
2531  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2532  << ", level " << m_currLevel+LEVEL_REF_ID
2533  << ", step " << m_currStep
2534  << ":"
2535  << std::endl;
2536  for (unsigned int i = 0; i < prevChain.subSequenceSize(); ++i) {
2537  prevChain.getPositionValues(i,prevPosition);
2538  *m_env.subDisplayFile() << " prevChain[" << i
2539  << "] = " << prevPosition
2540  << ", prevLogLikelihoodValues[" << i
2541  << "] = " << prevLogLikelihoodValues[i]
2542  << ", prevLogTargetValues[" << i
2543  << "] = " << prevLogTargetValues[i]
2544  << std::endl;
2545  }
2546  }
2547 #endif
2548 
2549  unsigned int quantity1 = prevChain.unifiedSequenceSize();
2550  unsigned int quantity2 = currChain.unifiedSequenceSize();
2551  unsigned int quantity3 = prevLogLikelihoodValues.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2552  unsigned int quantity4 = currLogLikelihoodValues.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2553  unsigned int quantity5 = prevLogTargetValues.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2554  unsigned int quantity6 = currLogTargetValues.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2555  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2556  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2557  << ", level " << m_currLevel+LEVEL_REF_ID
2558  << ", step " << m_currStep
2559  << ": prevChain.unifiedSequenceSize() = " << quantity1
2560  << ", currChain.unifiedSequenceSize() = " << quantity2
2561  << ", prevLogLikelihoodValues.unifiedSequenceSize() = " << quantity3
2562  << ", currLogLikelihoodValues.unifiedSequenceSize() = " << quantity4
2563  << ", prevLogTargetValues.unifiedSequenceSize() = " << quantity5
2564  << ", currLogTargetValues.unifiedSequenceSize() = " << quantity6
2565  << std::endl;
2566  }
2567 
2568  queso_require_equal_to_msg(prevChain.subSequenceSize(), prevLogLikelihoodValues.subSequenceSize(), "different sizes between previous chain and previous sequence of likelihood values");
2569 
2570  queso_require_equal_to_msg(prevChain.subSequenceSize(), prevLogTargetValues.subSequenceSize(), "different sizes between previous chain and previous sequence of target values");
2571 
2572  // Set 'indexOfFirstWeight' and 'indexOfLastWeight' // KAUST
2573  indexOfFirstWeight = 0;
2574  indexOfLastWeight = indexOfFirstWeight + prevChain.subSequenceSize()-1;
2575  {
2576  //std::cout << "m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc() << std::endl;
2577  int r = m_env.inter0Rank();
2578  //std::cout << "r = " << r << std::endl;
2579  m_env.inter0Comm().Barrier();
2580  unsigned int auxUint = 0;
2581  if (r > 0) {
2582  RawType_MPI_Status status;
2583  //std::cout << "Rank " << r << " is entering MPI_Recv()" << std::endl;
2584  m_env.inter0Comm().Recv((void*) &auxUint, 1, RawValue_MPI_UNSIGNED, r-1, r-1, &status,
2585  "MLSampling<P_V,P_M>::generateSequence()",
2586  "failed MPI.Recv()");
2587  //std::cout << "Rank " << r << " received auxUint = " << auxUint << std::endl;
2588  indexOfFirstWeight = auxUint;
2589  indexOfLastWeight = indexOfFirstWeight + prevChain.subSequenceSize()-1;
2590  }
2591  if (r < (m_env.inter0Comm().NumProc()-1)) {
2592  auxUint = indexOfLastWeight + 1;
2593  //std::cout << "Rank " << r << " is sending auxUint = " << auxUint << std::endl;
2594  m_env.inter0Comm().Send((void*) &auxUint, 1, RawValue_MPI_UNSIGNED, r+1, r,
2595  "MLSampling<P_V,P_M>::generateSequence()",
2596  "failed MPI.Send()");
2597  //std::cout << "Rank " << r << " sent auxUint = " << auxUint << std::endl;
2598  }
2599  m_env.inter0Comm().Barrier();
2600  }
2601 
2602  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
2603  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2604  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2605  << ", level " << m_currLevel+LEVEL_REF_ID
2606  << ", step " << m_currStep
2607  << ", after " << stepRunTime << " seconds"
2608  << std::endl;
2609  }
2610 
2611  return;
2612 }
2613 //---------------------------------------------------
2614 template <class P_V,class P_M>
2615 void
2616 MLSampling<P_V,P_M>::generateSequence_Step03_inter0(
2617  const MLSamplingLevelOptions* currOptions, // input
2618  const ScalarSequence<double>& prevLogLikelihoodValues, // input
2619  double prevExponent, // input
2620  double failedExponent, // input // gpmsa1
2621  double& currExponent, // output
2622  ScalarSequence<double>& weightSequence) // output
2623 {
2624  int iRC = UQ_OK_RC;
2625  struct timeval timevalStep;
2626  iRC = gettimeofday(&timevalStep, NULL);
2627  if (iRC) {}; // just to remove compiler warning
2628 
2629  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2630  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2631  << ", level " << m_currLevel+LEVEL_REF_ID
2632  << ", step " << m_currStep
2633  << ", failedExponent = " << failedExponent // gpmsa1
2634  << ": beginning step 3 of 11"
2635  << std::endl;
2636  }
2637 
2638  std::vector<double> exponents(2,0.);
2639  exponents[0] = prevExponent;
2640  exponents[1] = 1.;
2641 
2642  double nowExponent = 1.; // Try '1.' right away
2643  double nowEffectiveSizeRatio = 0.; // To be computed
2644 
2645  unsigned int nowAttempt = 0;
2646  bool testResult = false;
2647  double meanEffectiveSizeRatio = .5*(currOptions->m_minEffectiveSizeRatio + currOptions->m_maxEffectiveSizeRatio);
2648  ScalarSequence<double> omegaLnDiffSequence(m_env,prevLogLikelihoodValues.subSequenceSize(),"");
2649 
2650  double nowUnifiedEvidenceLnFactor = 0.;
2651  do {
2652  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2653  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2654  << ", level " << m_currLevel+LEVEL_REF_ID
2655  << ", step " << m_currStep
2656  << ", failedExponent = " << failedExponent // gpmsa1
2657  << ": entering loop for computing next exponent"
2658  << ", with nowAttempt = " << nowAttempt
2659  << std::endl;
2660  }
2661 
2662  if (failedExponent > 0.) { // gpmsa1
2663  nowExponent = .5*(prevExponent+failedExponent);
2664  }
2665  else {
2666  if (nowAttempt > 0) {
2667  if (nowEffectiveSizeRatio > meanEffectiveSizeRatio) {
2668  exponents[0] = nowExponent;
2669  }
2670  else {
2671  exponents[1] = nowExponent;
2672  }
2673  nowExponent = .5*(exponents[0] + exponents[1]);
2674  }
2675  }
2676  double auxExponent = nowExponent;
2677  if (prevExponent != 0.) {
2678  auxExponent /= prevExponent;
2679  auxExponent -= 1.;
2680  }
2681  double subWeightRatioSum = 0.;
2682  double unifiedWeightRatioSum = 0.;
2683 
2684  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2685  omegaLnDiffSequence[i] = prevLogLikelihoodValues[i]*auxExponent; // likelihood is important
2686  }
2687 
2688 #if 1 // prudenci-2012-07-06
2689  //double unifiedOmegaLnMin = omegaLnDiffSequence.unifiedMinPlain(m_vectorSpace.numOfProcsForStorage() == 1);
2690  double unifiedOmegaLnMax = omegaLnDiffSequence.unifiedMaxPlain(m_vectorSpace.numOfProcsForStorage() == 1);
2691 #else
2692  double unifiedOmegaLnMin = 0.;
2693  double unifiedOmegaLnMax = 0.;
2694  omegaLnDiffSequence.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1, // KAUST3
2695  0,
2696  omegaLnDiffSequence.subSequenceSize(),
2697  unifiedOmegaLnMin,
2698  unifiedOmegaLnMax);
2699 #endif
2700  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2701  omegaLnDiffSequence[i] -= unifiedOmegaLnMax;
2702  weightSequence[i] = exp(omegaLnDiffSequence[i]);
2703  subWeightRatioSum += weightSequence[i];
2704 #if 0 // For debug only
2705  if ((m_currLevel == 1) && (nowAttempt == 6)) {
2706  if (m_env.subDisplayFile() && (m_env.displayVerbosity() >= 99)) {
2707  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2708  << ", level " << m_currLevel+LEVEL_REF_ID
2709  << ", step " << m_currStep
2710  << ", i = " << i
2711  << ", prevLogLikelihoodValues[i] = " << prevLogLikelihoodValues[i]
2712  << ", omegaLnDiffSequence[i] = " << omegaLnDiffSequence[i]
2713  << ", weightSequence[i] = " << weightSequence[i]
2714  //<< ", subWeightRatioSum = " << subWeightRatioSum
2715  << std::endl;
2716  }
2717  }
2718 #endif
2719  }
2720  m_env.inter0Comm().Allreduce((void *) &subWeightRatioSum, (void *) &unifiedWeightRatioSum, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
2721  "MLSampling<P_V,P_M>::generateSequence()",
2722  "failed MPI.Allreduce() for weight ratio sum");
2723 
2724  unsigned int auxQuantity = weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2725  nowUnifiedEvidenceLnFactor = log(unifiedWeightRatioSum) + unifiedOmegaLnMax - log(auxQuantity);
2726 
2727  double effectiveSampleSize = 0.;
2728  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2729  weightSequence[i] /= unifiedWeightRatioSum;
2730  effectiveSampleSize += weightSequence[i]*weightSequence[i];
2731  //if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2732  // *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2733  // << ", level " << m_currLevel+LEVEL_REF_ID
2734  // << ", step " << m_currStep
2735  // << ": i = " << i
2736  // << ", effectiveSampleSize = " << effectiveSampleSize
2737  // << std::endl;
2738  //}
2739  }
2740 
2741  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2742  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2743  << ", level " << m_currLevel+LEVEL_REF_ID
2744  << ", step " << m_currStep
2745  << ": nowAttempt = " << nowAttempt
2746  << ", prevExponent = " << prevExponent
2747  << ", exponents[0] = " << exponents[0]
2748  << ", nowExponent = " << nowExponent
2749  << ", exponents[1] = " << exponents[1]
2750  << ", subWeightRatioSum = " << subWeightRatioSum
2751  << ", unifiedWeightRatioSum = " << unifiedWeightRatioSum
2752  << ", unifiedOmegaLnMax = " << unifiedOmegaLnMax
2753  << ", weightSequence.unifiedSequenceSize() = " << auxQuantity
2754  << ", nowUnifiedEvidenceLnFactor = " << nowUnifiedEvidenceLnFactor
2755  << ", effectiveSampleSize = " << effectiveSampleSize
2756  << std::endl;
2757  }
2758 
2759 #if 0 // For debug only
2760  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2761  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2762  << ", level " << m_currLevel+LEVEL_REF_ID
2763  << ", step " << m_currStep
2764  << ":"
2765  << std::endl;
2766  }
2767  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2768  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2769  *m_env.subDisplayFile() << " weightSequence[" << i
2770  << "] = " << weightSequence[i]
2771  << std::endl;
2772  }
2773  }
2774 #endif
2775 
2776  double subQuantity = effectiveSampleSize;
2777  effectiveSampleSize = 0.;
2778  m_env.inter0Comm().Allreduce((void *) &subQuantity, (void *) &effectiveSampleSize, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
2779  "MLSampling<P_V,P_M>::generateSequence()",
2780  "failed MPI.Allreduce() for effective sample size");
2781 
2782  effectiveSampleSize = 1./effectiveSampleSize;
2783  nowEffectiveSizeRatio = effectiveSampleSize/((double) weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1));
2784  queso_require_less_equal_msg(nowEffectiveSizeRatio, (1.+1.e-8), "effective sample size ratio cannot be > 1");
2785 
2786  // m_env.worldRank(),
2787  // "MLSampling<P_V,P_M>::generateSequence()",
2788  // "effective sample size ratio cannot be < 1");
2789 
2790  if (failedExponent > 0.) { // gpmsa1
2791  testResult = true;
2792  }
2793  else {
2794  //bool aux1 = (nowEffectiveSizeRatio == meanEffectiveSizeRatio);
2795  bool aux2 = (nowExponent == 1. )
2796  &&
2797  (nowEffectiveSizeRatio > meanEffectiveSizeRatio);
2798  bool aux3 = (nowEffectiveSizeRatio >= currOptions->m_minEffectiveSizeRatio)
2799  &&
2800  (nowEffectiveSizeRatio <= currOptions->m_maxEffectiveSizeRatio);
2801  testResult = aux2 || aux3;
2802  }
2803 
2804  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2805  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2806  << ", level " << m_currLevel+LEVEL_REF_ID
2807  << ", step " << m_currStep
2808  << ": nowAttempt = " << nowAttempt
2809  << ", prevExponent = " << prevExponent
2810  << ", failedExponent = " << failedExponent // gpmsa1
2811  << ", exponents[0] = " << exponents[0]
2812  << ", nowExponent = " << nowExponent
2813  << ", exponents[1] = " << exponents[1]
2814  << ", effectiveSampleSize = " << effectiveSampleSize
2815  << ", weightSequenceSize = " << weightSequence.subSequenceSize()
2816  << ", minEffectiveSizeRatio = " << currOptions->m_minEffectiveSizeRatio
2817  << ", nowEffectiveSizeRatio = " << nowEffectiveSizeRatio
2818  << ", maxEffectiveSizeRatio = " << currOptions->m_maxEffectiveSizeRatio
2819  //<< ", aux2 = " << aux2
2820  //<< ", aux3 = " << aux3
2821  << ", testResult = " << testResult
2822  << std::endl;
2823  }
2824  nowAttempt++;
2825 
2826  // Make sure all nodes in 'inter0Comm' have the same value of 'nowExponent'
2827  if (MiscCheckForSameValueInAllNodes(nowExponent,
2828  0., // kept 'zero' on 2010/03/05
2829  m_env.inter0Comm(),
2830  "MLSampling<P_V,P_M>::generateSequence(), step 3, nowExponent") == false) {
2831  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2832  *m_env.subDisplayFile() << "WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2833  << ", level " << m_currLevel+LEVEL_REF_ID
2834  << ", step " << m_currStep
2835  << ": nowAttempt = " << nowAttempt
2836  << ", MiscCheck for 'nowExponent' detected a problem"
2837  << std::endl;
2838  }
2839  }
2840 
2841  // Make sure all nodes in 'inter0Comm' have the same value of 'testResult'
2842  if (MiscCheckForSameValueInAllNodes(testResult,
2843  0., // kept 'zero' on 2010/03/05
2844  m_env.inter0Comm(),
2845  "MLSampling<P_V,P_M>::generateSequence(), step 3, testResult") == false) {
2846  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2847  *m_env.subDisplayFile() << "WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2848  << ", level " << m_currLevel+LEVEL_REF_ID
2849  << ", step " << m_currStep
2850  << ": nowAttempt = " << nowAttempt
2851  << ", MiscCheck for 'testResult' detected a problem"
2852  << std::endl;
2853  }
2854  }
2855  } while (testResult == false);
2856  currExponent = nowExponent;
2857  if (failedExponent > 0.) { // gpmsa1
2858  m_logEvidenceFactors[m_logEvidenceFactors.size()-1] = nowUnifiedEvidenceLnFactor;
2859  }
2860  else {
2861  m_logEvidenceFactors.push_back(nowUnifiedEvidenceLnFactor); // restart
2862  }
2863 
2864  unsigned int quantity1 = weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
2865  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2866  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2867  << ", level " << m_currLevel+LEVEL_REF_ID
2868  << ", step " << m_currStep
2869  << ": weightSequence.subSequenceSize() = " << weightSequence.subSequenceSize()
2870  << ", weightSequence.unifiedSequenceSize() = " << quantity1
2871  << ", failedExponent = " << failedExponent // gpmsa1
2872  << ", currExponent = " << currExponent
2873  << ", effective ratio = " << nowEffectiveSizeRatio
2874  << ", log(evidence factor) = " << m_logEvidenceFactors[m_logEvidenceFactors.size()-1]
2875  << ", evidence factor = " << exp(m_logEvidenceFactors[m_logEvidenceFactors.size()-1])
2876  << std::endl;
2877 
2878  //unsigned int numZeros = 0;
2879  //for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2880  // *m_env.subDisplayFile() << "weightSequence[" << i
2881  // << "] = " << weightSequence[i]
2882  // << std::endl;
2883  // if (weightSequence[i] == 0.) numZeros++;
2884  //}
2885  //*m_env.subDisplayFile() << "Number of zeros in weightSequence = " << numZeros
2886  // << std::endl;
2887  }
2888 
2889  // Make sure all nodes in 'inter0Comm' have the same value of 'logEvidenceFactor'
2890  if (MiscCheckForSameValueInAllNodes(m_logEvidenceFactors[m_logEvidenceFactors.size()-1],
2891  3.0e-16, // changed from 'zero' on 2010/03/03
2892  m_env.inter0Comm(),
2893  "MLSampling<P_V,P_M>::generateSequence(), step 3, logEvidenceFactor") == false) {
2894  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2895  *m_env.subDisplayFile() << "WARNING, In MLSampling<P_V,P_M>::generateSequence()"
2896  << ", level " << m_currLevel+LEVEL_REF_ID
2897  << ", step " << m_currStep
2898  << ", failedExponent = " << failedExponent // gpmsa1
2899  << ": nowAttempt = " << nowAttempt
2900  << ", MiscCheck for 'logEvidenceFactor' detected a problem"
2901  << std::endl;
2902  }
2903  }
2904 
2905  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
2906  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2907  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
2908  << ", level " << m_currLevel+LEVEL_REF_ID
2909  << ", step " << m_currStep
2910  << ", failedExponent = " << failedExponent // gpmsa
2911  << ", after " << stepRunTime << " seconds"
2912  << std::endl;
2913  }
2914 
2915  return;
2916 }
2917 //---------------------------------------------------
2918 template <class P_V,class P_M>
2919 void
2920 MLSampling<P_V,P_M>::generateSequence_Step04_inter0(
2921  const SequenceOfVectors<P_V,P_M>& prevChain, // input
2922  const ScalarSequence<double>& weightSequence, // input
2923  P_M& unifiedCovMatrix) // output
2924 {
2925  int iRC = UQ_OK_RC;
2926  struct timeval timevalStep;
2927  iRC = gettimeofday(&timevalStep, NULL);
2928  if (iRC) {}; // just to remove compiler warning
2929 
2930  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2931  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2932  << ", level " << m_currLevel+LEVEL_REF_ID
2933  << ", step " << m_currStep
2934  << ": beginning step 4 of 11"
2935  << std::endl;
2936  }
2937 
2938  P_V auxVec(m_vectorSpace.zeroVector());
2939  P_V subWeightedMeanVec(m_vectorSpace.zeroVector());
2940  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2941  prevChain.getPositionValues(i,auxVec);
2942  subWeightedMeanVec += weightSequence[i]*auxVec;
2943  }
2944 
2945  // Todd Oliver 2010-09-07: compute weighted mean over all processors
2946  P_V unifiedWeightedMeanVec(m_vectorSpace.zeroVector());
2947  if (m_env.inter0Rank() >= 0) {
2948  subWeightedMeanVec.mpiAllReduce(RawValue_MPI_SUM,m_env.inter0Comm(),unifiedWeightedMeanVec);
2949  }
2950  else {
2951  unifiedWeightedMeanVec = subWeightedMeanVec;
2952  }
2953 
2954  P_V diffVec(m_vectorSpace.zeroVector());
2955  P_M subCovMatrix(m_vectorSpace.zeroVector());
2956  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
2957  prevChain.getPositionValues(i,auxVec);
2958  diffVec = auxVec - unifiedWeightedMeanVec;
2959  subCovMatrix += weightSequence[i]*matrixProduct(diffVec,diffVec);
2960  }
2961 
2962  for (unsigned int i = 0; i < unifiedCovMatrix.numRowsLocal(); ++i) { // KAUST5
2963  for (unsigned int j = 0; j < unifiedCovMatrix.numCols(); ++j) {
2964  double localValue = subCovMatrix(i,j);
2965  double sumValue = 0.;
2966  if (m_env.inter0Rank() >= 0) {
2967  m_env.inter0Comm().Allreduce((void *) &localValue, (void *) &sumValue, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
2968  "MLSampling<P_V,P_M>::generateSequence()",
2969  "failed MPI.Allreduce() for cov matrix");
2970  }
2971  else {
2972  sumValue = localValue;
2973  }
2974  unifiedCovMatrix(i,j) = sumValue;
2975  }
2976  }
2977 
2978  if (m_numDisabledParameters > 0) { // gpmsa2
2979  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
2980  if (m_parameterEnabledStatus[paramId] == false) {
2981  for (unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
2982  unifiedCovMatrix(i,paramId) = 0.;
2983  }
2984  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
2985  unifiedCovMatrix(paramId,j) = 0.;
2986  }
2987  unifiedCovMatrix(paramId,paramId) = 1.;
2988  }
2989  }
2990  }
2991 
2992  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2993  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
2994  << ", level " << m_currLevel+LEVEL_REF_ID
2995  << ", step " << m_currStep
2996  << ": unifiedCovMatrix = " << unifiedCovMatrix
2997  << std::endl;
2998  }
2999 
3000  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3001  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3002  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3003  << ", level " << m_currLevel+LEVEL_REF_ID
3004  << ", step " << m_currStep
3005  << ", after " << stepRunTime << " seconds"
3006  << std::endl;
3007  }
3008 
3009  return;
3010 }
3011 //---------------------------------------------------
3012 template <class P_V,class P_M>
3013 void
3014 MLSampling<P_V,P_M>::generateSequence_Step05_inter0(
3015  unsigned int unifiedRequestedNumSamples, // input
3016  const ScalarSequence<double>& weightSequence, // input
3017  std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // output
3018  std::vector<double>& unifiedWeightStdVectorAtProc0Only) // output
3019 {
3020  int iRC = UQ_OK_RC;
3021  struct timeval timevalStep;
3022  iRC = gettimeofday(&timevalStep, NULL);
3023  if (iRC) {}; // just to remove compiler warning
3024 
3025  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3026  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3027  << ", level " << m_currLevel+LEVEL_REF_ID
3028  << ", step " << m_currStep
3029  << ": beginning step 5 of 11"
3030  << std::endl;
3031  }
3032 
3033 #if 0 // For debug only
3034  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3035  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3036  << ", level " << m_currLevel+LEVEL_REF_ID
3037  << ", step " << m_currStep
3038  << ", before weightSequence.getUnifiedContentsAtProc0Only()"
3039  << ":"
3040  << std::endl;
3041  }
3042  for (unsigned int i = 0; i < weightSequence.subSequenceSize(); ++i) {
3043  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3044  *m_env.subDisplayFile() << ", weightSequence[" << i
3045  << "] = " << weightSequence[i]
3046  << std::endl;
3047  }
3048  }
3049 #endif
3050 
3051  weightSequence.getUnifiedContentsAtProc0Only(m_vectorSpace.numOfProcsForStorage() == 1,
3052  unifiedWeightStdVectorAtProc0Only);
3053 
3054 #if 0 // For debug only
3055  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3056  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3057  << ", level " << m_currLevel+LEVEL_REF_ID
3058  << ", step " << m_currStep
3059  << ", after weightSequence.getUnifiedContentsAtProc0Only()"
3060  << ":"
3061  << std::endl;
3062  }
3063  for (unsigned int i = 0; i < unifiedWeightStdVectorAtProc0Only.size(); ++i) {
3064  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3065  *m_env.subDisplayFile() << " unifiedWeightStdVectorAtProc0Only[" << i
3066  << "] = " << unifiedWeightStdVectorAtProc0Only[i]
3067  << std::endl;
3068  }
3069  }
3070 #endif
3071  sampleIndexes_proc0(unifiedRequestedNumSamples, // input
3072  unifiedWeightStdVectorAtProc0Only, // input
3073  unifiedIndexCountersAtProc0Only); // output
3074 
3075  unsigned int auxUnifiedSize = weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3076  if (m_env.inter0Rank() == 0) {
3077  queso_require_equal_to_msg(unifiedIndexCountersAtProc0Only.size(), auxUnifiedSize, "wrong output from sampleIndexesAtProc0() in step 5");
3078  }
3079 
3080  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3081  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3082  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3083  << ", level " << m_currLevel+LEVEL_REF_ID
3084  << ", step " << m_currStep
3085  << ", after " << stepRunTime << " seconds"
3086  << std::endl;
3087  }
3088 
3089  return;
3090 }
3091 //---------------------------------------------------
3092 template <class P_V,class P_M>
3093 void
3094 MLSampling<P_V,P_M>::generateSequence_Step06_all(
3095  const MLSamplingLevelOptions* currOptions, // input
3096  unsigned int indexOfFirstWeight, // input
3097  unsigned int indexOfLastWeight, // input
3098  const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // input
3099  bool& useBalancedChains, // output
3100  std::vector<ExchangeInfoStruct>& exchangeStdVec) // output
3101 {
3102  int iRC = UQ_OK_RC;
3103  struct timeval timevalStep;
3104  iRC = gettimeofday(&timevalStep, NULL);
3105  if (iRC) {}; // just to remove compiler warning
3106 
3107  useBalancedChains = decideOnBalancedChains_all(currOptions, // input
3108  indexOfFirstWeight, // input
3109  indexOfLastWeight, // input
3110  unifiedIndexCountersAtProc0Only, // input
3111  exchangeStdVec); // output
3112 
3113  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3114  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3115  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3116  << ", level " << m_currLevel+LEVEL_REF_ID
3117  << ", step " << m_currStep
3118  << ", after " << stepRunTime << " seconds"
3119  << std::endl;
3120  }
3121 
3122  return;
3123 }
3124 //---------------------------------------------------
3125 template <class P_V,class P_M>
3126 void
3127 MLSampling<P_V,P_M>::generateSequence_Step07_inter0(
3128  bool useBalancedChains, // input
3129  unsigned int indexOfFirstWeight, // input
3130  unsigned int indexOfLastWeight, // input
3131  const std::vector<unsigned int>& unifiedIndexCountersAtProc0Only, // input
3132  UnbalancedLinkedChainsPerNodeStruct& unbalancedLinkControl, // (possible) output
3133  const MLSamplingLevelOptions* currOptions, // input
3134  const SequenceOfVectors<P_V,P_M>& prevChain, // input
3135  double prevExponent, // input
3136  double currExponent, // input
3137  const ScalarSequence<double>& prevLogLikelihoodValues, // input
3138  const ScalarSequence<double>& prevLogTargetValues, // input
3139  std::vector<ExchangeInfoStruct>& exchangeStdVec, // (possible) input/output
3140  BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl) // (possible) output
3141 {
3142  int iRC = UQ_OK_RC;
3143  struct timeval timevalStep;
3144  iRC = gettimeofday(&timevalStep, NULL);
3145  if (iRC) {}; // just to remove compiler warning
3146 
3147  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3148  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3149  << ", level " << m_currLevel+LEVEL_REF_ID
3150  << ", step " << m_currStep
3151  << ": beginning step 7 of 11"
3152  << std::endl;
3153  }
3154 
3155  if (useBalancedChains) {
3156  prepareBalLinkedChains_inter0(currOptions, // input
3157  prevChain, // input
3158  prevExponent, // input
3159  currExponent, // input
3160  prevLogLikelihoodValues, // input
3161  prevLogTargetValues, // input
3162  exchangeStdVec, // input/output
3163  balancedLinkControl); // output
3164  }
3165  else {
3166  prepareUnbLinkedChains_inter0(indexOfFirstWeight, // input
3167  indexOfLastWeight, // input
3168  unifiedIndexCountersAtProc0Only, // input
3169  unbalancedLinkControl); // output
3170  }
3171 
3172  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3173  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3174  << ", level " << m_currLevel+LEVEL_REF_ID
3175  << ", step " << m_currStep
3176  << ": balancedLinkControl.balLinkedChains.size() = " << balancedLinkControl.balLinkedChains.size()
3177  << ", unbalancedLinkControl.unbLinkedChains.size() = " << unbalancedLinkControl.unbLinkedChains.size()
3178  << std::endl;
3179  }
3180 
3181  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3182  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3183  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3184  << ", level " << m_currLevel+LEVEL_REF_ID
3185  << ", step " << m_currStep
3186  << ", after " << stepRunTime << " seconds"
3187  << std::endl;
3188  }
3189 
3190  return;
3191 }
3192 //---------------------------------------------------
3193 template <class P_V,class P_M>
3194 void
3195 MLSampling<P_V,P_M>::generateSequence_Step08_all(
3196  BayesianJointPdf<P_V,P_M>& currPdf, // input/output
3197  GenericVectorRV<P_V,P_M>& currRv) // output
3198 {
3199  int iRC = UQ_OK_RC;
3200  struct timeval timevalStep;
3201  iRC = gettimeofday(&timevalStep, NULL);
3202  if (iRC) {}; // just to remove compiler warning
3203 
3204  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3205  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3206  << ", level " << m_currLevel+LEVEL_REF_ID
3207  << ", step " << m_currStep
3208  << ": beginning step 8 of 11"
3209  << std::endl;
3210  }
3211 
3212  currRv.setPdf(currPdf);
3213 
3214  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3215  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3216  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3217  << ", level " << m_currLevel+LEVEL_REF_ID
3218  << ", step " << m_currStep
3219  << ", after " << stepRunTime << " seconds"
3220  << std::endl;
3221  }
3222 
3223  return;
3224 }
3225 //---------------------------------------------------
3226 template <class P_V,class P_M>
3227 void
3228 MLSampling<P_V,P_M>::generateSequence_Step09_all(
3229  const SequenceOfVectors<P_V,P_M>& prevChain, // input
3230  double prevExponent, // input
3231  double currExponent, // input
3232  const ScalarSequence<double>& prevLogLikelihoodValues, // input
3233  const ScalarSequence<double>& prevLogTargetValues, // input
3234  unsigned int indexOfFirstWeight, // input
3235  unsigned int indexOfLastWeight, // input
3236  const std::vector<double>& unifiedWeightStdVectorAtProc0Only, // input
3237  const ScalarSequence<double>& weightSequence, // input
3238  double prevEta, // input
3239  const GenericVectorRV<P_V,P_M>& currRv, // input
3240  MLSamplingLevelOptions* currOptions, // input (changed temporarily internally)
3241  P_M& unifiedCovMatrix, // input/output
3242  double& currEta) // output
3243 {
3244  int iRC = UQ_OK_RC;
3245  struct timeval timevalStep;
3246  iRC = gettimeofday(&timevalStep, NULL);
3247  if (iRC) {}; // just to remove compiler warning
3248 
3249  if (currOptions->m_scaleCovMatrix == false) {
3250  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3251  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3252  << ", level " << m_currLevel+LEVEL_REF_ID
3253  << ", step " << m_currStep
3254  << ": skipping step 9 of 11"
3255  << std::endl;
3256  }
3257  }
3258  else {
3259  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3260  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3261  << ", level " << m_currLevel+LEVEL_REF_ID
3262  << ", step " << m_currStep
3263  << ": beginning step 9 of 11"
3264  << std::endl;
3265  }
3266 
3267  double beforeEta = prevEta;
3268  double beforeRejectionRate = 0.; // To be updated
3269  bool beforeRejectionRateIsBelowRange = true; // To be updated
3270 
3271  double nowEta = prevEta;
3272  double nowRejectionRate = 0.; // To be computed
3273  bool nowRejectionRateIsBelowRange = true; // To be computed
3274 
3275  std::vector<double> etas(2,0.);
3276  etas[0] = beforeEta;
3277  etas[1] = 1.;
3278 
3279  std::vector<double> rejs(2,0.);
3280  rejs[0] = 0.; // To be computed
3281  rejs[1] = 0.; // To be computed
3282 
3283  unsigned int nowAttempt = 0;
3284  bool testResult = false;
3285  double meanRejectionRate = .5*(currOptions->m_minRejectionRate + currOptions->m_maxRejectionRate);
3286  bool useMiddlePointLogicForEta = false;
3287  P_M nowCovMatrix(unifiedCovMatrix);
3288 #if 0 // KAUST, to check
3289  std::vector<double> unifiedWeightStdVectorAtProc0Only(0);
3290  weightSequence.getUnifiedContentsAtProc0Only(m_vectorSpace.numOfProcsForStorage() == 1,
3291  unifiedWeightStdVectorAtProc0Only);
3292 #endif
3293  do {
3294  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3295  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3296  << ", level " << m_currLevel+LEVEL_REF_ID
3297  << ", step " << m_currStep
3298  << ": entering loop for assessing rejection rate"
3299  << ", with nowAttempt = " << nowAttempt
3300  << ", nowRejectionRate = " << nowRejectionRate
3301  << std::endl;
3302  }
3303  nowCovMatrix = unifiedCovMatrix;
3304 
3305  if (nowRejectionRate < currOptions->m_minRejectionRate) {
3306  nowRejectionRateIsBelowRange = true;
3307  }
3308  else if (nowRejectionRate > currOptions->m_maxRejectionRate) {
3309  nowRejectionRateIsBelowRange = false;
3310  }
3311  else {
3312  queso_error_msg("nowRejectionRate should be out of the requested range at this point of the logic");
3313  }
3314 
3315  if (m_env.inter0Rank() >= 0) { // KAUST
3316  if (nowAttempt > 0) {
3317  if (useMiddlePointLogicForEta == false) {
3318  if (nowAttempt == 1) {
3319  // Ok, keep useMiddlePointLogicForEta = false
3320  }
3321  else if ((beforeRejectionRateIsBelowRange == true) &&
3322  (nowRejectionRateIsBelowRange == true)) {
3323  // Ok
3324  }
3325  else if ((beforeRejectionRateIsBelowRange == false) &&
3326  (nowRejectionRateIsBelowRange == false)) {
3327  // Ok
3328  }
3329  else if ((beforeRejectionRateIsBelowRange == true ) &&
3330  (nowRejectionRateIsBelowRange == false)) {
3331  useMiddlePointLogicForEta = true;
3332 
3333  // This is the first time the middle point logic will be used below
3334  etas[0] = std::min(beforeEta,nowEta);
3335  etas[1] = std::max(beforeEta,nowEta);
3336 
3337  if (etas[0] == beforeEta) {
3338  rejs[0] = beforeRejectionRate;
3339  rejs[1] = nowRejectionRate;
3340  }
3341  else {
3342  rejs[0] = nowRejectionRate;
3343  rejs[1] = beforeRejectionRate;
3344  }
3345  }
3346  else if ((beforeRejectionRateIsBelowRange == false) &&
3347  (nowRejectionRateIsBelowRange == true )) {
3348  useMiddlePointLogicForEta = true;
3349 
3350  // This is the first time the middle point logic will be used below
3351  etas[0] = std::min(beforeEta,nowEta);
3352  etas[1] = std::max(beforeEta,nowEta);
3353  }
3354  else {
3355  queso_error_msg("before and now range flags are inconsistent");
3356  }
3357  } // if (useMiddlePointLogicForEta == false)
3358 
3359  beforeEta = nowEta;
3360  beforeRejectionRate = nowRejectionRate;
3361  beforeRejectionRateIsBelowRange = nowRejectionRateIsBelowRange;
3362  if (useMiddlePointLogicForEta == false) {
3363  if (beforeRejectionRateIsBelowRange) nowEta *= 4.;
3364  else nowEta /= 4.;
3365  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3366  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3367  << ", level " << m_currLevel+LEVEL_REF_ID
3368  << ", step " << m_currStep
3369  << ": in loop for assessing rejection rate"
3370  << ", with nowAttempt = " << nowAttempt
3371  << ", useMiddlePointLogicForEta = false"
3372  << ", nowEta just updated to value (to be tested) " << nowEta
3373  << std::endl;
3374  }
3375  }
3376  else {
3377  if (nowRejectionRate > meanRejectionRate) {
3378  if (rejs[0] > meanRejectionRate) {
3379  etas[0] = nowEta;
3380  etas[1] = etas[1];
3381  }
3382  else {
3383  etas[0] = etas[0];
3384  etas[1] = nowEta;
3385  }
3386  }
3387  else {
3388  if (rejs[0] < meanRejectionRate) {
3389  etas[0] = nowEta;
3390  etas[1] = etas[1];
3391  }
3392  else {
3393  etas[0] = etas[0];
3394  etas[1] = nowEta;
3395  }
3396  }
3397  nowEta = .5*(etas[0] + etas[1]);
3398  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3399  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3400  << ", level " << m_currLevel+LEVEL_REF_ID
3401  << ", step " << m_currStep
3402  << ": in loop for assessing rejection rate"
3403  << ", with nowAttempt = " << nowAttempt
3404  << ", useMiddlePointLogicForEta = true"
3405  << ", nowEta just updated to value (to be tested) " << nowEta
3406  << ", etas[0] = " << etas[0]
3407  << ", etas[1] = " << etas[1]
3408  << std::endl;
3409  }
3410  }
3411  } // if (nowAttempt > 0)
3412  } // if (m_env.inter0Rank() >= 0) // KAUST
3413 
3414  nowCovMatrix *= nowEta;
3415 
3416  // prudencio 2010-12-09: logic 'originalSubNumSamples += 1' added because of the difference of results between GNU and INTEL compiled codes
3417  double doubSubNumSamples = (1.-meanRejectionRate)/meanRejectionRate/currOptions->m_covRejectionRate/currOptions->m_covRejectionRate; // e.g. 19.99...; or 20.0; or 20.1; or 20.9
3418  unsigned int originalSubNumSamples = 1 + (unsigned int) (doubSubNumSamples); // e.g. 20; or 21; or 21; or 21
3419  double auxDouble = (double) originalSubNumSamples; // e.g. 20.0; or 21.0; or 21.0; or 21.0
3420  if ((auxDouble - doubSubNumSamples) < 1.e-8) { // e.g. 0.00...01; or 1.0; or 0.9; or 0.1
3421  originalSubNumSamples += 1;
3422  }
3423 
3424  if (m_env.inter0Rank() >= 0) { // KAUST
3425  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3426  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3427  << ", level " << m_currLevel+LEVEL_REF_ID
3428  << ", step " << m_currStep
3429  << ": in loop for assessing rejection rate"
3430  << ", about to sample " << originalSubNumSamples << " indexes"
3431  << ", meanRejectionRate = " << meanRejectionRate
3432  << ", covRejectionRate = " << currOptions->m_covRejectionRate
3433  << std::endl;
3434  }
3435  } // KAUST
3436 
3437  std::vector<unsigned int> nowUnifiedIndexCountersAtProc0Only(0); // It will be resized by 'sampleIndexes_proc0()' below
3438  if (m_env.inter0Rank() >= 0) { // KAUST
3439  unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3440  sampleIndexes_proc0(tmpUnifiedNumSamples, // input
3441  unifiedWeightStdVectorAtProc0Only, // input
3442  nowUnifiedIndexCountersAtProc0Only); // output
3443 
3444  unsigned int auxUnifiedSize = weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3445  if (m_env.inter0Rank() == 0) {
3446  queso_require_equal_to_msg(nowUnifiedIndexCountersAtProc0Only.size(), auxUnifiedSize, "wrong output from sampleIndexesAtProc0() in step 9");
3447  }
3448 
3449  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3450  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3451  << ", level " << m_currLevel+LEVEL_REF_ID
3452  << ", step " << m_currStep
3453  << ": in loop for assessing rejection rate"
3454  << ", about to distribute sampled assessment indexes"
3455  << std::endl;
3456  }
3457  } // KAUST
3458 
3459  std::vector<ExchangeInfoStruct> exchangeStdVec(0);
3460  BalancedLinkedChainsPerNodeStruct<P_V> nowBalLinkControl;
3461  UnbalancedLinkedChainsPerNodeStruct nowUnbLinkControl; // KAUST
3462 
3463  // All processors should call this routine in order to have the same decision value
3464  bool useBalancedChains = decideOnBalancedChains_all(currOptions, // input
3465  indexOfFirstWeight, // input
3466  indexOfLastWeight, // input
3467  nowUnifiedIndexCountersAtProc0Only, // input
3468  exchangeStdVec); // output
3469 
3470  if (m_env.inter0Rank() >= 0) { // KAUST
3471  if (useBalancedChains) {
3472  prepareBalLinkedChains_inter0(currOptions, // input
3473  prevChain, // input
3474  prevExponent, // input
3475  currExponent, // input
3476  prevLogLikelihoodValues, // input
3477  prevLogTargetValues, // input
3478  exchangeStdVec, // input/output
3479  nowBalLinkControl); // output
3480  }
3481  else {
3482  prepareUnbLinkedChains_inter0(indexOfFirstWeight, // input
3483  indexOfLastWeight, // input
3484  nowUnifiedIndexCountersAtProc0Only, // input
3485  nowUnbLinkControl); // output
3486  }
3487  } // KAUST
3488 
3489  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3490  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3491  << ", level " << m_currLevel+LEVEL_REF_ID
3492  << ", step " << m_currStep
3493  << ": in loop for assessing rejection rate"
3494  << ", about to generate assessment chain"
3495  << std::endl;
3496  }
3497 
3498  SequenceOfVectors<P_V,P_M> nowChain(m_vectorSpace,
3499  0,
3500  m_options.m_prefix+"now_chain");
3501  double nowRunTime = 0.;
3502  unsigned int nowRejections = 0;
3503 
3504  // KAUST: all nodes should call here
3505  bool savedTotallyMute = currOptions->m_totallyMute; // HERE - ENHANCEMENT
3506  unsigned int savedRawChainSize = currOptions->m_rawChainSize; // Ok to use rawChainSize
3507 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3508  bool savedRawChainComputeStats = currOptions->m_rawChainComputeStats;
3509 #endif
3510  bool savedFilteredChainGenerate = currOptions->m_filteredChainGenerate;
3511  unsigned int savedDrMaxNumExtraStages = currOptions->m_drMaxNumExtraStages;
3512  unsigned int savedAmAdaptInterval = currOptions->m_amAdaptInterval;
3513 
3514  currOptions->m_totallyMute = true;
3515  if (m_env.displayVerbosity() >= 999999) {
3516  currOptions->m_totallyMute = false;
3517  }
3518  currOptions->m_rawChainSize = 0; // will be set inside generateXYZLinkedChains()
3519 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3520  currOptions->m_rawChainComputeStats = false;
3521 #endif
3522  currOptions->m_filteredChainGenerate = false;
3523  currOptions->m_drMaxNumExtraStages = 0;
3524  currOptions->m_amAdaptInterval = 0;
3525 
3526  // KAUST: all nodes in 'subComm' should call here, important
3527  if (useBalancedChains) {
3528  generateBalLinkedChains_all(*currOptions, // input, only m_rawChainSize changes
3529  nowCovMatrix, // input
3530  currRv, // input
3531  nowBalLinkControl, // input // Round Rock
3532  nowChain, // output
3533  nowRunTime, // output
3534  nowRejections, // output
3535  NULL, // output
3536  NULL); // output
3537  }
3538  else {
3539  generateUnbLinkedChains_all(*currOptions, // input, only m_rawChainSize changes
3540  nowCovMatrix, // input
3541  currRv, // input
3542  nowUnbLinkControl, // input // Round Rock
3543  indexOfFirstWeight, // input // Round Rock
3544  prevChain, // input // Round Rock
3545  prevExponent, // input
3546  currExponent, // input
3547  prevLogLikelihoodValues, // input
3548  prevLogTargetValues, // input
3549  nowChain, // output
3550  nowRunTime, // output
3551  nowRejections, // output
3552  NULL, // output
3553  NULL); // output
3554  }
3555 
3556  // KAUST: all nodes should call here
3557  currOptions->m_totallyMute = savedTotallyMute;
3558  currOptions->m_rawChainSize = savedRawChainSize;
3559 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3560  currOptions->m_rawChainComputeStats = savedRawChainComputeStats;
3561 #endif
3562  currOptions->m_filteredChainGenerate = savedFilteredChainGenerate; // FIX ME
3563  currOptions->m_drMaxNumExtraStages = savedDrMaxNumExtraStages;
3564  currOptions->m_amAdaptInterval = savedAmAdaptInterval;
3565 
3566  for (unsigned int i = 0; i < nowBalLinkControl.balLinkedChains.size(); ++i) {
3567  queso_require_msg(nowBalLinkControl.balLinkedChains[i].initialPosition, "Initial position pointer in step 9 should not be NULL");
3568  delete nowBalLinkControl.balLinkedChains[i].initialPosition;
3569  nowBalLinkControl.balLinkedChains[i].initialPosition = NULL;
3570  }
3571  nowBalLinkControl.balLinkedChains.clear();
3572 
3573  if (m_env.inter0Rank() >= 0) { // KAUST
3574  // If only one cov matrix is used, then the rejection should be assessed among all inter0Comm nodes // KAUST3
3575  unsigned int nowUnifiedRejections = 0;
3576  m_env.inter0Comm().Allreduce((void *) &nowRejections, (void *) &nowUnifiedRejections, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
3577  "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3578  "failed MPI.Allreduce() for now rejections");
3579 
3580 #if 0 // Round Rock 2009 12 29
3581  unsigned int tmpUnifiedNumSamples = 0;
3582  m_env.inter0Comm().Allreduce((void *) &tmpSubNumSamples, (void *) &tmpUnifiedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
3583  "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3584  "failed MPI.Allreduce() for num samples in step 9");
3585 #endif
3586 
3587  unsigned int tmpUnifiedNumSamples = originalSubNumSamples*m_env.inter0Comm().NumProc();
3588  nowRejectionRate = ((double) nowUnifiedRejections) / ((double) tmpUnifiedNumSamples);
3589 
3590  //bool aux1 = (nowRejectionRate == meanRejectionRate);
3591  bool aux2 = (nowRejectionRate >= currOptions->m_minRejectionRate)
3592  &&
3593  (nowRejectionRate <= currOptions->m_maxRejectionRate);
3594  testResult = aux2;
3595 
3596  // Make sure all nodes in 'inter0Comm' have the same value of 'testResult'
3597  if (MiscCheckForSameValueInAllNodes(testResult,
3598  0., // kept 'zero' on 2010/03/03
3599  m_env.inter0Comm(),
3600  "MLSampling<P_V,P_M>::generateSequence_Step09_all(), step 9, testResult") == false) {
3601  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3602  *m_env.subDisplayFile() << "WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3603  << ", level " << m_currLevel+LEVEL_REF_ID
3604  << ", step " << m_currStep
3605  << ": nowAttempt = " << nowAttempt
3606  << ", MiscCheck for 'testResult' detected a problem"
3607  << std::endl;
3608  }
3609  }
3610  } // if (m_env.inter0Rank() >= 0) { // KAUST
3611 
3612  // KAUST: all nodes in 'subComm' should have the same 'testResult'
3613  unsigned int tmpUint = (unsigned int) testResult;
3614  m_env.subComm().Bcast((void *) &tmpUint, (int) 1, RawValue_MPI_UNSIGNED, 0, // Yes, 'subComm', important
3615  "MLSampling<P_V,P_M>::generateSequence_Step09_all()",
3616  "failed MPI.Bcast() for testResult");
3617  testResult = (bool) tmpUint;
3618 
3619  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3620  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3621  << ", level " << m_currLevel+LEVEL_REF_ID
3622  << ", step " << m_currStep
3623  << ": in loop for assessing rejection rate"
3624  << ", nowAttempt = " << nowAttempt
3625  << ", beforeEta = " << beforeEta
3626  << ", etas[0] = " << etas[0]
3627  << ", nowEta = " << nowEta
3628  << ", etas[1] = " << etas[1]
3629  << ", minRejectionRate = " << currOptions->m_minRejectionRate
3630  << ", nowRejectionRate = " << nowRejectionRate
3631  << ", maxRejectionRate = " << currOptions->m_maxRejectionRate
3632  << std::endl;
3633  }
3634  nowAttempt++;
3635 
3636  if (m_env.inter0Rank() >= 0) { // KAUST
3637  // Make sure all nodes in 'inter0Comm' have the same value of 'nowEta'
3638  if (MiscCheckForSameValueInAllNodes(nowEta,
3639  1.0e-16, // changed from 'zero' on 2009/11/dd
3640  m_env.inter0Comm(),
3641  "MLSampling<P_V,P_M>::generateSequence_Step09_all(), step 9, nowEta") == false) {
3642  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3643  *m_env.subDisplayFile() << "WARNING, In MLSampling<P_V,P_M>::generateSequence()"
3644  << ", level " << m_currLevel+LEVEL_REF_ID
3645  << ", step " << m_currStep
3646  << ": nowAttempt = " << nowAttempt
3647  << ", MiscCheck for 'nowEta' detected a problem"
3648  << std::endl;
3649  }
3650  }
3651  }
3652  } while (testResult == false);
3653  currEta = nowEta;
3654  if (currEta != 1.) {
3655  unifiedCovMatrix *= currEta;
3656  if (m_numDisabledParameters > 0) { // gpmsa2
3657  for (unsigned int paramId = 0; paramId < m_vectorSpace.dimLocal(); ++paramId) {
3658  if (m_parameterEnabledStatus[paramId] == false) {
3659  for (unsigned int i = 0; i < m_vectorSpace.dimLocal(); ++i) {
3660  unifiedCovMatrix(i,paramId) = 0.;
3661  }
3662  for (unsigned int j = 0; j < m_vectorSpace.dimLocal(); ++j) {
3663  unifiedCovMatrix(paramId,j) = 0.;
3664  }
3665  unifiedCovMatrix(paramId,paramId) = 1.;
3666  }
3667  }
3668  }
3669  }
3670 
3671  unsigned int quantity1 = weightSequence.unifiedSequenceSize(m_vectorSpace.numOfProcsForStorage() == 1);
3672  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3673  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step09_all()"
3674  << ", level " << m_currLevel+LEVEL_REF_ID
3675  << ", step " << m_currStep
3676  << ": weightSequence.subSequenceSize() = " << weightSequence.subSequenceSize()
3677  << ", weightSequence.unifiedSequenceSize() = " << quantity1
3678  << ", currEta = " << currEta
3679  << ", assessed rejection rate = " << nowRejectionRate
3680  << std::endl;
3681  }
3682  }
3683 
3684  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3685  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3686  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3687  << ", level " << m_currLevel+LEVEL_REF_ID
3688  << ", step " << m_currStep
3689  << ", after " << stepRunTime << " seconds"
3690  << std::endl;
3691  }
3692 
3693  return;
3694 }
3695 //---------------------------------------------------
3696 template <class P_V,class P_M>
3697 void
3698 MLSampling<P_V,P_M>::generateSequence_Step10_all(
3699  MLSamplingLevelOptions& currOptions, // input (changed temporarily internally)
3700  const P_M& unifiedCovMatrix, // input
3701  const GenericVectorRV <P_V,P_M>& currRv, // input
3702  bool useBalancedChains, // input
3703  const UnbalancedLinkedChainsPerNodeStruct& unbalancedLinkControl, // input // Round Rock
3704  unsigned int indexOfFirstWeight, // input // Round Rock
3705  const SequenceOfVectors<P_V,P_M>& prevChain, // input // Round Rock
3706  double prevExponent, // input
3707  double currExponent, // input
3708  const ScalarSequence<double>& prevLogLikelihoodValues, // input
3709  const ScalarSequence<double>& prevLogTargetValues, // input
3710  const BalancedLinkedChainsPerNodeStruct<P_V>& balancedLinkControl, // input // Round Rock
3711  SequenceOfVectors <P_V,P_M>& currChain, // output
3712  double& cumulativeRawChainRunTime, // output
3713  unsigned int& cumulativeRawChainRejections, // output
3714  ScalarSequence <double>* currLogLikelihoodValues, // output
3715  ScalarSequence <double>* currLogTargetValues) // output
3716 {
3717  int iRC = UQ_OK_RC;
3718  struct timeval timevalStep;
3719  iRC = gettimeofday(&timevalStep, NULL);
3720  if (iRC) {}; // just to remove compiler warning
3721 
3722  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3723  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3724  << ", level " << m_currLevel+LEVEL_REF_ID
3725  << ", step " << m_currStep
3726  << ": beginning step 10 of 11"
3727  << ", currLogLikelihoodValues = " << currLogLikelihoodValues
3728  << std::endl;
3729  }
3730 
3731  // All nodes should call here
3732  bool savedTotallyMute = currOptions.m_totallyMute; // HERE - ENHANCEMENT
3733  unsigned int savedRawChainSize = currOptions.m_rawChainSize; // Ok to use rawChainSize
3734 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3735  bool savedRawChainComputeStats = currOptions.m_rawChainComputeStats;
3736 #endif
3737  bool savedFilteredChainGenerate = currOptions.m_filteredChainGenerate;
3738 
3739  currOptions.m_totallyMute = true;
3740  if (m_env.displayVerbosity() >= 999999) {
3741  currOptions.m_totallyMute = false;
3742  }
3743  currOptions.m_rawChainSize = 0; // will be set inside generateXYZLinkedChains()
3744 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3745  currOptions.m_rawChainComputeStats = false;
3746 #endif
3747  currOptions.m_filteredChainGenerate = false;
3748 
3749  // All nodes should call here
3750  if (useBalancedChains) {
3751  generateBalLinkedChains_all(currOptions, // input, only m_rawChainSize changes
3752  unifiedCovMatrix, // input
3753  currRv, // input
3754  balancedLinkControl, // input // Round Rock
3755  currChain, // output
3756  cumulativeRawChainRunTime, // output
3757  cumulativeRawChainRejections, // output
3758  currLogLikelihoodValues, // output // likelihood is important
3759  currLogTargetValues); // output
3760  }
3761  else {
3762  generateUnbLinkedChains_all(currOptions, // input, only m_rawChainSize changes
3763  unifiedCovMatrix, // input
3764  currRv, // input
3765  unbalancedLinkControl, // input // Round Rock
3766  indexOfFirstWeight, // input // Round Rock
3767  prevChain, // input // Round Rock
3768  prevExponent, // input
3769  currExponent, // input
3770  prevLogLikelihoodValues, // input
3771  prevLogTargetValues, // input
3772  currChain, // output
3773  cumulativeRawChainRunTime, // output
3774  cumulativeRawChainRejections, // output
3775  currLogLikelihoodValues, // output // likelihood is important
3776  currLogTargetValues); // output
3777  }
3778 
3779  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3780  double tmpValue = INFINITY;
3781  if (currLogLikelihoodValues) tmpValue = (*currLogLikelihoodValues)[0];
3782  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3783  << ", level " << m_currLevel+LEVEL_REF_ID
3784  << ", step " << m_currStep
3785  << ", after chain generatrion"
3786  << ", currLogLikelihoodValues[0] = " << tmpValue
3787  << std::endl;
3788  }
3789 
3790  // All nodes should call here
3791  currOptions.m_totallyMute = savedTotallyMute;
3792  currOptions.m_rawChainSize = savedRawChainSize;
3793 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3794  currOptions.m_rawChainComputeStats = savedRawChainComputeStats;
3795 #endif
3796  currOptions.m_filteredChainGenerate = savedFilteredChainGenerate; // FIX ME
3797 
3798  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3799  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3800  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3801  << ", level " << m_currLevel+LEVEL_REF_ID
3802  << ", step " << m_currStep
3803  << ", after " << stepRunTime << " seconds"
3804  << std::endl;
3805  }
3806 
3807  return;
3808 }
3809 //---------------------------------------------------
3810 template <class P_V,class P_M>
3811 void
3812 MLSampling<P_V,P_M>::generateSequence_Step11_inter0(
3813  const MLSamplingLevelOptions* currOptions, // input
3814  unsigned int unifiedRequestedNumSamples, // input
3815  unsigned int cumulativeRawChainRejections, // input
3816  SequenceOfVectors<P_V,P_M>& currChain, // input/output
3817  ScalarSequence<double>& currLogLikelihoodValues, // input/output
3818  ScalarSequence<double>& currLogTargetValues, // input/output
3819  unsigned int& unifiedNumberOfRejections) // output
3820 {
3821  int iRC = UQ_OK_RC;
3822  struct timeval timevalStep;
3823  iRC = gettimeofday(&timevalStep, NULL);
3824  if (iRC) {}; // just to remove compiler warning
3825 
3826  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3827  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
3828  << ", level " << m_currLevel+LEVEL_REF_ID
3829  << ", step " << m_currStep
3830  << ": beginning step 11 of 11"
3831  << std::endl;
3832  }
3833 
3834  //if (m_env.subComm().MyPID() == 0) std::cout << "Aqui 000" << std::endl;
3835 
3836 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3837  if (currOptions->m_rawChainComputeStats) {
3838  FilePtrSetStruct filePtrSet;
3839  m_env.openOutputFile(currOptions->m_dataOutputFileName,
3840  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
3841  currOptions->m_dataOutputAllowedSet,
3842  false,
3843  filePtrSet);
3844 
3845  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) { // output debug
3846  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3847  << ", level " << m_currLevel+LEVEL_REF_ID
3848  << ", step " << m_currStep
3849  << ", calling computeStatistics for raw chain"
3850  << ". Ofstream pointer value = " << filePtrSet.ofsVar
3851  << ", statistical options are"
3852  << "\n" << *currOptions->m_rawChainStatisticalOptionsObj
3853  << std::endl;
3854  }
3855  //m_env.syncPrintDebugMsg("At step 11, calling computeStatistics for raw chain",1,10,m_env.inter0Comm()); // output debug
3856  currChain.computeStatistics(*currOptions->m_rawChainStatisticalOptionsObj,
3857  filePtrSet.ofsVar);
3858 
3859  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
3860  }
3861  // Compute MLE and MAP
3862  // rr0
3863 #endif
3864  if (currOptions->m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
3865  currChain.unifiedWriteContents(currOptions->m_rawChainDataOutputFileName,
3866  currOptions->m_rawChainDataOutputFileType); // KAUST5
3867  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3868  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3869  << ", level " << m_currLevel+LEVEL_REF_ID
3870  << ", step " << m_currStep
3871  << ", before calling currLogLikelihoodValues.unifiedWriteContents()"
3872  << ", currLogLikelihoodValues[0] = " << currLogLikelihoodValues[0]
3873  << std::endl;
3874  }
3875  currLogLikelihoodValues.unifiedWriteContents(currOptions->m_rawChainDataOutputFileName,
3876  currOptions->m_rawChainDataOutputFileType);
3877  currLogTargetValues.unifiedWriteContents (currOptions->m_rawChainDataOutputFileName,
3878  currOptions->m_rawChainDataOutputFileType);
3879  }
3880 
3881  if (currOptions->m_filteredChainGenerate) {
3882  FilePtrSetStruct filePtrSet;
3883  m_env.openOutputFile(currOptions->m_dataOutputFileName,
3884  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
3885  currOptions->m_dataOutputAllowedSet,
3886  false,
3887  filePtrSet);
3888 
3889  unsigned int filterInitialPos = (unsigned int) (currOptions->m_filteredChainDiscardedPortion * (double) currChain.subSequenceSize());
3890  unsigned int filterSpacing = currOptions->m_filteredChainLag;
3891  if (filterSpacing == 0) {
3892  currChain.computeFilterParams(filePtrSet.ofsVar,
3893  filterInitialPos,
3894  filterSpacing);
3895  }
3896 
3897  // Filter positions from the converged portion of the chain
3898  currChain.filter(filterInitialPos,
3899  filterSpacing);
3900  currChain.setName(currOptions->m_prefix + "filtChain");
3901 
3902  currLogLikelihoodValues.filter(filterInitialPos,
3903  filterSpacing);
3904  currLogLikelihoodValues.setName(currOptions->m_prefix + "filtLogLikelihood");
3905 
3906  currLogTargetValues.filter(filterInitialPos,
3907  filterSpacing);
3908  currLogTargetValues.setName(currOptions->m_prefix + "filtLogTarget");
3909 
3910 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
3911  if (currOptions->m_filteredChainComputeStats) {
3912  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) { // output debug
3913  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence_Step()"
3914  << ", level " << m_currLevel+LEVEL_REF_ID
3915  << ", step " << m_currStep
3916  << ", calling computeStatistics for filtered chain"
3917  << ". Ofstream pointer value = " << filePtrSet.ofsVar
3918  << ", statistical options are"
3919  << "\n" << *currOptions->m_rawChainStatisticalOptionsObj
3920  << std::endl;
3921  }
3922 
3923  //m_env.syncPrintDebugMsg("At step 11, calling computeStatistics for filtered chain",1,10,m_env.inter0Comm()); // output debug
3924  currChain.computeStatistics(*currOptions->m_filteredChainStatisticalOptionsObj,
3925  filePtrSet.ofsVar);
3926 
3927  }
3928 #endif
3929  // Compute MLE and MAP
3930  // rr0
3931  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
3932 
3933  if (currOptions->m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
3934  currChain.unifiedWriteContents (currOptions->m_filteredChainDataOutputFileName,
3935  currOptions->m_filteredChainDataOutputFileType);
3936  currLogLikelihoodValues.unifiedWriteContents(currOptions->m_filteredChainDataOutputFileName,
3937  currOptions->m_filteredChainDataOutputFileType);
3938  currLogTargetValues.unifiedWriteContents (currOptions->m_filteredChainDataOutputFileName,
3939  currOptions->m_filteredChainDataOutputFileType);
3940  }
3941  } // if (currOptions->m_filteredChainGenerate)
3942 
3943  if (currOptions->m_filteredChainGenerate) {
3944  // Do not check
3945  }
3946  else {
3947  // Check if unified size of generated chain matches the unified requested size // KAUST
3948  unsigned int tmpSize = currChain.subSequenceSize();
3949  unsigned int unifiedGeneratedNumSamples = 0;
3950  m_env.inter0Comm().Allreduce((void *) &tmpSize, (void *) &unifiedGeneratedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
3951  "MLSampling<P_V,P_M>::generateSequence()",
3952  "failed MPI.Allreduce() for generated num samples in step 11");
3953  //std::cout << "unifiedGeneratedNumSamples = " << unifiedGeneratedNumSamples
3954  // << ", unifiedRequestedNumSamples = " << unifiedRequestedNumSamples
3955  // << std::endl;
3956  queso_require_equal_to_msg(unifiedGeneratedNumSamples, unifiedRequestedNumSamples, "currChain (linked one) has been generated with invalid size");
3957  }
3958 
3959  // Compute unified number of rejections
3960  m_env.inter0Comm().Allreduce((void *) &cumulativeRawChainRejections, (void *) &unifiedNumberOfRejections, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
3961  "MLSampling<P_V,P_M>::generateSequence()",
3962  "failed MPI.Allreduce() for number of rejections");
3963 
3964  double stepRunTime = MiscGetEllapsedSeconds(&timevalStep);
3965  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3966  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence_Step()"
3967  << ", level " << m_currLevel+LEVEL_REF_ID
3968  << ", step " << m_currStep
3969  << ", after " << stepRunTime << " seconds"
3970  << std::endl;
3971  }
3972 
3973  return;
3974 }
3975 
3976 // Default constructor -----------------------------
3977 template<class P_V,class P_M>
3978 MLSampling<P_V,P_M>::MLSampling(
3979  const char* prefix,
3980  const BaseVectorRV <P_V,P_M>& priorRv,
3981  const BaseScalarFunction<P_V,P_M>& likelihoodFunction)
3982  :
3983  m_env (priorRv.env()),
3984  m_priorRv (priorRv),
3985  m_likelihoodFunction(likelihoodFunction),
3986  m_vectorSpace (m_priorRv.imageSet().vectorSpace()),
3987  m_targetDomain (InstantiateIntersection(m_priorRv.pdf().domainSet(),m_likelihoodFunction.domainSet())),
3988  m_numDisabledParameters (0), // gpmsa2
3989  m_parameterEnabledStatus(m_vectorSpace.dimLocal(),true), // gpmsa2
3990  m_options (m_env,prefix),
3991  m_currLevel (0),
3992  m_currStep (0),
3993  m_debugExponent (0.),
3994  m_logEvidenceFactors(0),
3995  m_logEvidence (0.),
3996  m_meanLogLikelihood (0.),
3997  m_eig (0.)
3998 {
3999  if (m_env.subDisplayFile()) {
4000  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::constructor()"
4001  << std::endl;
4002  }
4003 
4004  if (m_env.subDisplayFile()) {
4005  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::constructor()"
4006  << std::endl;
4007  }
4008 }
4009 // Destructor ---------------------------------------
4010 template<class P_V,class P_M>
4011 MLSampling<P_V,P_M>::~MLSampling()
4012 {
4013  m_numDisabledParameters = 0; // gpmsa2
4014  m_parameterEnabledStatus.clear(); // gpmsa2
4015  if (m_targetDomain) delete m_targetDomain;
4016 }
4017 // Statistical methods-------------------------------
4018 /* This operation currently implements the PAMSSA algorithm (S. H. Cheung and E. E. Prudencio. Parallel adaptive multilevel
4019  * sampling algorithms for the Bayesian analysis of mathematical models. International Journal
4020  * for Uncertainty Quantification, 2(3):215237, 2012.)*/
4021 template <class P_V,class P_M>
4022 void
4023 MLSampling<P_V,P_M>::generateSequence(
4024  BaseVectorSequence<P_V,P_M>& workingChain,
4025  ScalarSequence<double>* workingLogLikelihoodValues,
4026  ScalarSequence<double>* workingLogTargetValues)
4027 {
4028  struct timeval timevalRoutineBegin;
4029  int iRC = 0;
4030  iRC = gettimeofday(&timevalRoutineBegin, NULL);
4031  if (iRC) {}; // just to remove compiler warning
4032 
4033  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4034  *m_env.subDisplayFile() << "Entering MLSampling<P_V,P_M>::generateSequence()"
4035  << ", at " << ctime(&timevalRoutineBegin.tv_sec)
4036  << ", after " << timevalRoutineBegin.tv_sec - m_env.timevalBegin().tv_sec
4037  << " seconds from queso environment instatiation..."
4038  << std::endl;
4039  }
4040 
4041  //***********************************************************
4042  // Declaration of Variables
4043  //***********************************************************
4044  double currExponent = 0.; // restate
4045  double currEta = 1.; // restate
4046  unsigned int currUnifiedRequestedNumSamples = 0;
4047  SequenceOfVectors<P_V,P_M> currChain (m_vectorSpace, // restate
4048  0,
4049  m_options.m_prefix+"curr_chain");
4050  ScalarSequence<double> currLogLikelihoodValues(m_env, // restate
4051  0,
4052  "");
4053  ScalarSequence<double> currLogTargetValues (m_env, // restate
4054  0,
4055  "");
4056 
4057  bool stopAtEndOfLevel = false;
4058  char levelPrefix[256];
4059 
4060  //***********************************************************
4061  // Take care of first level (level '0')
4062  //***********************************************************
4063  MLSamplingLevelOptions defaultLevelOptions(m_env,(m_options.m_prefix + "default_").c_str());
4064  defaultLevelOptions.scanOptionsValues(NULL);
4065 
4066  MLSamplingLevelOptions lastLevelOptions(m_env,(m_options.m_prefix + "last_").c_str());
4067  lastLevelOptions.scanOptionsValues(&defaultLevelOptions);
4068 
4069  if (m_options.m_restartInput_baseNameForFiles != ".") {
4070  restartML(currExponent, // output
4071  currEta, // output
4072  currChain, // output
4073  currLogLikelihoodValues, // output
4074  currLogTargetValues); // output
4075 
4076  if (currExponent == 1.) {
4077  if (lastLevelOptions.m_parameterDisabledSet.size() > 0) { // gpmsa2
4078  for (std::set<unsigned int>::iterator setIt = lastLevelOptions.m_parameterDisabledSet.begin(); setIt != lastLevelOptions.m_parameterDisabledSet.end(); ++setIt) {
4079  unsigned int paramId = *setIt;
4080  if (paramId < m_vectorSpace.dimLocal()) {
4081  m_numDisabledParameters++;
4082  m_parameterEnabledStatus[paramId] = false;
4083  }
4084  }
4085  }
4086 
4087 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4088  if (lastLevelOptions.m_rawChainComputeStats) {
4089  FilePtrSetStruct filePtrSet;
4090  m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4091  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
4092  lastLevelOptions.m_dataOutputAllowedSet,
4093  false,
4094  filePtrSet);
4095 
4096  currChain.computeStatistics(*lastLevelOptions.m_rawChainStatisticalOptionsObj,
4097  filePtrSet.ofsVar);
4098 
4099  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
4100  }
4101 #endif
4102  // Compute MLE and MAP
4103  // rr0
4104 
4105  if (lastLevelOptions.m_rawChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
4106  currChain.unifiedWriteContents (lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4107  currLogLikelihoodValues.unifiedWriteContents(lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4108  currLogTargetValues.unifiedWriteContents (lastLevelOptions.m_rawChainDataOutputFileName,lastLevelOptions.m_rawChainDataOutputFileType);
4109  }
4110 
4111  if (lastLevelOptions.m_filteredChainGenerate) {
4112  FilePtrSetStruct filePtrSet;
4113  m_env.openOutputFile(lastLevelOptions.m_dataOutputFileName,
4114  UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT, // Yes, always ".m"
4115  lastLevelOptions.m_dataOutputAllowedSet,
4116  false,
4117  filePtrSet);
4118 
4119  unsigned int filterInitialPos = (unsigned int) (lastLevelOptions.m_filteredChainDiscardedPortion * (double) currChain.subSequenceSize());
4120  unsigned int filterSpacing = lastLevelOptions.m_filteredChainLag;
4121  if (filterSpacing == 0) {
4122  currChain.computeFilterParams(filePtrSet.ofsVar,
4123  filterInitialPos,
4124  filterSpacing);
4125  }
4126 
4127  // Filter positions from the converged portion of the chain
4128  currChain.filter(filterInitialPos,
4129  filterSpacing);
4130  currChain.setName(lastLevelOptions.m_prefix + "filtChain");
4131 
4132  currLogLikelihoodValues.filter(filterInitialPos,
4133  filterSpacing);
4134  currLogLikelihoodValues.setName(lastLevelOptions.m_prefix + "filtLogLikelihood");
4135 
4136  currLogTargetValues.filter(filterInitialPos,
4137  filterSpacing);
4138  currLogTargetValues.setName(lastLevelOptions.m_prefix + "filtLogTarget");
4139 
4140 #ifdef QUESO_USES_SEQUENCE_STATISTICAL_OPTIONS
4141  if (lastLevelOptions.m_filteredChainComputeStats) {
4142  currChain.computeStatistics(*lastLevelOptions.m_filteredChainStatisticalOptionsObj,
4143  filePtrSet.ofsVar);
4144  }
4145 #endif
4146  // Compute MLE and MAP
4147  // rr0
4148  m_env.closeFile(filePtrSet,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT);
4149 
4150  if (lastLevelOptions.m_filteredChainDataOutputFileName != UQ_MH_SG_FILENAME_FOR_NO_FILE) {
4151  currChain.unifiedWriteContents (lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4152  currLogLikelihoodValues.unifiedWriteContents(lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4153  currLogTargetValues.unifiedWriteContents (lastLevelOptions.m_filteredChainDataOutputFileName,lastLevelOptions.m_filteredChainDataOutputFileType);
4154  }
4155  } // if (lastLevelOptions.m_filteredChainGenerate)
4156  }
4157  }
4158  else {
4159  sprintf(levelPrefix,"%d_",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
4160  MLSamplingLevelOptions currOptions(m_env,(m_options.m_prefix + levelPrefix).c_str());
4161  currOptions.scanOptionsValues(&defaultLevelOptions);
4162 
4163  if (currOptions.m_parameterDisabledSet.size() > 0) { // gpmsa2
4164  for (std::set<unsigned int>::iterator setIt = currOptions.m_parameterDisabledSet.begin(); setIt != currOptions.m_parameterDisabledSet.end(); ++setIt) {
4165  unsigned int paramId = *setIt;
4166  if (paramId < m_vectorSpace.dimLocal()) {
4167  m_numDisabledParameters++;
4168  m_parameterEnabledStatus[paramId] = false;
4169  }
4170  }
4171  }
4172 
4173  generateSequence_Level0_all(currOptions, // input
4174  currUnifiedRequestedNumSamples, // output
4175  currChain, // output
4176  currLogLikelihoodValues, // output
4177  currLogTargetValues); // output
4178 
4179  stopAtEndOfLevel = currOptions.m_stopAtEnd;
4180  bool performCheckpoint = stopAtEndOfLevel;
4181  if (m_options.m_restartOutput_levelPeriod > 0) {
4182  performCheckpoint = performCheckpoint || ( ((m_currLevel + 1) % m_options.m_restartOutput_levelPeriod) == 0 );
4183  }
4184  if (performCheckpoint) {
4185  checkpointML(currExponent, // input
4186  currEta, // input
4187  currChain, // input
4188  currLogLikelihoodValues, // input
4189  currLogTargetValues); // input
4190  }
4191  }
4192  //std::cout << "In QUESO: end of level 0. Exiting on purpose" << std::endl;
4193  //exit(1);
4194 
4195  double minLogLike = -INFINITY;
4196  double maxLogLike = INFINITY;
4197  currLogLikelihoodValues.subMinMaxExtra(0,
4198  currLogLikelihoodValues.subSequenceSize(),
4199  minLogLike,
4200  maxLogLike);
4201  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4202  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4203  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4204  << ", sub minLogLike = " << minLogLike
4205  << ", sub maxLogLike = " << maxLogLike
4206  << std::endl;
4207  }
4208 
4209  m_env.fullComm().Barrier();
4210 
4211  minLogLike = -INFINITY;
4212  maxLogLike = INFINITY;
4213  currLogLikelihoodValues.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4214  0,
4215  currLogLikelihoodValues.subSequenceSize(),
4216  minLogLike,
4217  maxLogLike);
4218  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4219  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4220  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4221  << ", unified minLogLike = " << minLogLike
4222  << ", unified maxLogLike = " << maxLogLike
4223  << std::endl;
4224  }
4225 
4226  //***********************************************************
4227  // Take care of next levels
4228  //***********************************************************
4229  while ((currExponent < 1. ) && // begin level while
4230  (stopAtEndOfLevel == false)) {
4231  m_currLevel++; // restate
4232 
4233  struct timeval timevalLevelBegin;
4234  iRC = 0;
4235  iRC = gettimeofday(&timevalLevelBegin, NULL);
4236 
4237  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4238  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4239  << ": beginning level " << m_currLevel+LEVEL_REF_ID
4240  << ", at " << ctime(&timevalLevelBegin.tv_sec)
4241  << ", after " << timevalLevelBegin.tv_sec - timevalRoutineBegin.tv_sec
4242  << " seconds from entering the routine"
4243  << ", after " << timevalLevelBegin.tv_sec - m_env.timevalBegin().tv_sec
4244  << " seconds from queso environment instatiation"
4245  << std::endl;
4246  }
4247 
4248  iRC = UQ_OK_RC;
4249  struct timeval timevalLevel;
4250  iRC = gettimeofday(&timevalLevel, NULL);
4251  double cumulativeRawChainRunTime = 0.;
4252  unsigned int cumulativeRawChainRejections = 0;
4253 
4254  bool tryExponentEta = true; // gpmsa1
4255  double failedExponent = 0.; // gpmsa1
4256  double failedEta = 0.; // gpmsa1
4257 
4258  MLSamplingLevelOptions* currOptions = NULL; // step 1
4259  SequenceOfVectors<P_V,P_M>* prevChain = NULL; // step 2
4260  ScalarSequence<double> prevLogLikelihoodValues(m_env, 0, "");
4261  ScalarSequence<double> prevLogTargetValues (m_env, 0, "");
4262  double prevExponent = 0.; // step 3
4263  unsigned int indexOfFirstWeight = 0; // step 2
4264  unsigned int indexOfLastWeight = 0; // step 2
4265  P_M* unifiedCovMatrix = NULL; // step 4
4266  bool useBalancedChains = false; // step 6
4267  BalancedLinkedChainsPerNodeStruct<P_V>* balancedLinkControl = NULL; // step 7
4268  UnbalancedLinkedChainsPerNodeStruct* unbalancedLinkControl = NULL; // step 7
4269  BayesianJointPdf<P_V,P_M>* currPdf = NULL; // step 8
4270  GenericVectorRV<P_V,P_M>* currRv = NULL; // step 8
4271 
4272  unsigned int exponentEtaTriedAmount = 0;
4273  while (tryExponentEta) { // gpmsa1
4274  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4275  *m_env.subDisplayFile() << "In IMLSampling<P_V,P_M>::generateSequence()"
4276  << ", level " << m_currLevel+LEVEL_REF_ID
4277  << ", beginning 'do-while(tryExponentEta)"
4278  << ": failedExponent = " << failedExponent
4279  << ", failedEta = " << failedEta
4280  << std::endl;
4281  }
4282 
4283  //***********************************************************
4284  // Step 1 of 11: read options
4285  //***********************************************************
4286  m_currStep = 1;
4287  sprintf(levelPrefix,"%d_",m_currLevel+LEVEL_REF_ID); // Yes, '+0'
4288  currOptions = new MLSamplingLevelOptions(m_env,(m_options.m_prefix + levelPrefix).c_str());
4289  currOptions->scanOptionsValues(&defaultLevelOptions);
4290 
4291  if (currOptions->m_parameterDisabledSet.size() > 0) { // gpmsa2
4292  for (std::set<unsigned int>::iterator setIt = currOptions->m_parameterDisabledSet.begin(); setIt != currOptions->m_parameterDisabledSet.end(); ++setIt) {
4293  unsigned int paramId = *setIt;
4294  if (paramId < m_vectorSpace.dimLocal()) {
4295  m_numDisabledParameters++;
4296  m_parameterEnabledStatus[paramId] = false;
4297  }
4298  }
4299  }
4300 
4301  if (m_env.inter0Rank() >= 0) {
4302  generateSequence_Step01_inter0(currOptions, // input
4303  currUnifiedRequestedNumSamples); // output
4304  }
4305 
4306  //***********************************************************
4307  // Step 2 of 11: save [chain and corresponding target pdf values] from previous level
4308  //***********************************************************
4309  m_currStep = 2;
4310  prevExponent = currExponent;
4311  double prevEta = currEta;
4312  unsigned int prevUnifiedRequestedNumSamples = currUnifiedRequestedNumSamples;
4313  prevChain = new SequenceOfVectors<P_V,P_M>(m_vectorSpace,
4314  0,
4315  m_options.m_prefix+"prev_chain");
4316 
4317  indexOfFirstWeight = 0;
4318  indexOfLastWeight = 0;
4319 
4320  //std::cout << "m_env.inter0Rank() = " << m_env.inter0Rank() << std::endl;
4321  if (m_env.inter0Rank() >= 0) {
4322  generateSequence_Step02_inter0(currOptions, // input
4323  currChain, // input/output // restate
4324  currLogLikelihoodValues, // input/output // restate
4325 
4326  currLogTargetValues, // input/output // restate
4327  *prevChain, // output
4328  prevLogLikelihoodValues, // output
4329  prevLogTargetValues, // output
4330  indexOfFirstWeight, // output
4331  indexOfLastWeight); // output
4332  }
4333 
4334  //***********************************************************
4335  // Step 3 of 11: compute [currExponent and sequence of weights] for current level
4336  // update 'm_logEvidenceFactors'
4337  //***********************************************************
4338  m_currStep = 3;
4339  ScalarSequence<double> weightSequence(m_env,prevLogLikelihoodValues.subSequenceSize(),"");
4340  if (m_env.inter0Rank() >= 0) {
4341  generateSequence_Step03_inter0(currOptions, // input
4342  prevLogLikelihoodValues, // input
4343  prevExponent, // input
4344  failedExponent, // input // gpmsa1
4345  currExponent, // output
4346  weightSequence); // output
4347  }
4348 
4349  // All nodes in 'subComm' should have the same 'currExponent'
4350  m_env.subComm().Bcast((void *) &currExponent, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm', important
4351  "MLSampling<P_V,P_M>::generateSequence()",
4352  "failed MPI.Bcast() for currExponent");
4353  m_debugExponent = currExponent;
4354 
4355  if (currExponent == 1.) {
4356  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4357  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4358  << ", level " << m_currLevel+LEVEL_REF_ID
4359  << ", step " << m_currStep
4360  << ": copying 'last' level options to current options" // In all nodes of 'subComm', important
4361  << std::endl;
4362  }
4363  delete currOptions;
4364  currOptions = &lastLevelOptions;
4365 
4366  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4367  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4368  << ", level " << m_currLevel+LEVEL_REF_ID
4369  << ", step " << m_currStep
4370  << ": after copying 'last' level options to current options, the current options are"
4371  << "\n" << *currOptions
4372  << std::endl;
4373  }
4374 
4375  if (m_env.inter0Rank() >= 0) {
4376  // It is necessary to recompute 'currUnifiedRequestedNumSamples' because
4377  // 'currOptions' has just been replaced by 'lastLevelOptions'
4378  unsigned int tmpSize = currOptions->m_rawChainSize;
4379  m_env.inter0Comm().Allreduce((void *) &tmpSize, (void *) &currUnifiedRequestedNumSamples, (int) 1, RawValue_MPI_UNSIGNED, RawValue_MPI_SUM,
4380  "MLSampling<P_V,P_M>::generateSequence()",
4381  "failed MPI.Allreduce() for requested num samples in step 3");
4382  }
4383  }
4384 
4385  //***********************************************************
4386  // Step 4 of 11: create covariance matrix for current level
4387  //***********************************************************
4388  m_currStep = 4;
4389  P_V oneVec(m_vectorSpace.zeroVector());
4390  oneVec.cwSet(1.);
4391 
4392  unifiedCovMatrix = NULL;
4393  if (m_env.inter0Rank() >= 0) {
4394  unifiedCovMatrix = m_vectorSpace.newMatrix();
4395  }
4396  else {
4397  unifiedCovMatrix = new P_M(oneVec);
4398  }
4399 
4400  if (m_env.inter0Rank() >= 0) {
4401  generateSequence_Step04_inter0(*prevChain, // input
4402  weightSequence, // input
4403  *unifiedCovMatrix); // output
4404  }
4405 
4406  //***********************************************************
4407  // Step 5 of 11: create *unified* finite distribution for current level
4408  //***********************************************************
4409  m_currStep = 5;
4410  std::vector<unsigned int> unifiedIndexCountersAtProc0Only(0);
4411  std::vector<double> unifiedWeightStdVectorAtProc0Only(0); // KAUST, to check
4412  if (m_env.inter0Rank() >= 0) {
4413  generateSequence_Step05_inter0(currUnifiedRequestedNumSamples, // input
4414  weightSequence, // input
4415  unifiedIndexCountersAtProc0Only, // output
4416  unifiedWeightStdVectorAtProc0Only); // output
4417  }
4418 
4419  //***********************************************************
4420  // Step 6 of 11: decide on using balanced chains or not
4421  //***********************************************************
4422  m_currStep = 6;
4423  useBalancedChains = false;
4424  std::vector<ExchangeInfoStruct> exchangeStdVec(0);
4425  // All processors should call this routine in order to have the same decision value
4426  generateSequence_Step06_all(currOptions, // input
4427  indexOfFirstWeight, // input
4428  indexOfLastWeight, // input
4429  unifiedIndexCountersAtProc0Only, // input
4430  useBalancedChains, // output
4431  exchangeStdVec); // output
4432 
4433  //***********************************************************
4434  // Step 7 of 11: plan for number of linked chains for each node so that all
4435  // nodes generate the closest possible to the same number of positions
4436  //***********************************************************
4437  m_currStep = 7;
4438  balancedLinkControl = new BalancedLinkedChainsPerNodeStruct<P_V>();
4439  unbalancedLinkControl = new UnbalancedLinkedChainsPerNodeStruct ();
4440  if (m_env.inter0Rank() >= 0) {
4441  generateSequence_Step07_inter0(useBalancedChains, // input
4442  indexOfFirstWeight, // input
4443  indexOfLastWeight, // input
4444  unifiedIndexCountersAtProc0Only, // input
4445  *unbalancedLinkControl, // (possible) output
4446  currOptions, // input
4447  *prevChain, // input
4448  prevExponent, // input
4449  currExponent, // input
4450  prevLogLikelihoodValues, // input
4451  prevLogTargetValues, // input
4452  exchangeStdVec, // (possible) input/output
4453  *balancedLinkControl); // (possible) output
4454  }
4455 
4456  // aqui: prevChain might not be needed anymore. Delete it to save memory
4457 
4458  //***********************************************************
4459  // Step 8 of 11: create vector RV for current level
4460  //***********************************************************
4461  m_currStep = 8;
4462  currPdf = new BayesianJointPdf<P_V,P_M> (m_options.m_prefix.c_str(),
4463  m_priorRv.pdf(),
4464  m_likelihoodFunction,
4465  currExponent,
4466  *m_targetDomain);
4467 
4468  currRv = new GenericVectorRV<P_V,P_M> (m_options.m_prefix.c_str(),
4469  *m_targetDomain);
4470 
4471  // All nodes should set 'currRv'
4472  generateSequence_Step08_all(*currPdf,
4473  *currRv);
4474 
4475  //***********************************************************
4476  // Step 9 of 11: scale unified covariance matrix until min <= rejection rate <= max
4477  //***********************************************************
4478  m_currStep = 9;
4479  generateSequence_Step09_all(*prevChain, // input
4480  prevExponent, // input
4481  currExponent, // input
4482  prevLogLikelihoodValues, // input
4483  prevLogTargetValues, // input
4484  indexOfFirstWeight, // input
4485  indexOfLastWeight, // input
4486  unifiedWeightStdVectorAtProc0Only, // input
4487  weightSequence, // input
4488  prevEta, // input
4489  *currRv, // input
4490  currOptions, // input (changed temporarily internally)
4491  *unifiedCovMatrix, // input/output
4492  currEta); // output
4493 
4494  tryExponentEta = false; // gpmsa1
4495  if ((currOptions->m_minAcceptableEta > 0. ) && // gpmsa1
4496  (currEta < currOptions->m_minAcceptableEta)) {
4497  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4498  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4499  << ", level " << m_currLevel+LEVEL_REF_ID
4500  << ", preparing to retry ExponentEta"
4501  << ": currExponent = " << currExponent
4502  << ", currEta = " << currEta
4503  << std::endl;
4504  }
4505  exponentEtaTriedAmount++;
4506  tryExponentEta = true;
4507  failedExponent = currExponent;
4508  failedEta = currEta;
4509 
4510  // "Return" to previous level
4511  delete currRv; // Step 8
4512  currRv = NULL; // Step 8
4513  delete currPdf; // Step 8
4514  currPdf = NULL; // Step 8
4515 
4516  delete balancedLinkControl; // Step 7
4517  balancedLinkControl = NULL; // Step 7
4518  delete unbalancedLinkControl; // Step 7
4519  unbalancedLinkControl = NULL; // Step 7
4520 
4521  delete unifiedCovMatrix; // Step 4
4522  unifiedCovMatrix = NULL; // Step 4
4523 
4524  currExponent = prevExponent;
4525  currEta = 1.; // prevEta;
4526  currUnifiedRequestedNumSamples = prevUnifiedRequestedNumSamples;
4527 
4528  currChain.clear(); // Step 2
4529  currChain = (*prevChain); // Step 2
4530  delete prevChain; // Step 2
4531  prevChain = NULL; // Step 2
4532 
4533  currLogLikelihoodValues = prevLogLikelihoodValues;
4534  currLogTargetValues = prevLogTargetValues;
4535  }
4536  } // while (tryExponentEta) // gpmsa1
4537 
4538  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) { // gpmsa1
4539  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4540  << ", level " << m_currLevel+LEVEL_REF_ID
4541  << ", exited 'do-while(tryExponentEta)"
4542  << ", failedExponent = " << failedExponent
4543  << ", failedEta = " << failedEta
4544  << std::endl;
4545  }
4546 
4547  //***********************************************************
4548  // Step 10 of 11: sample vector RV of current level
4549  //***********************************************************
4550  m_currStep = 10;
4551  // All nodes should call here
4552  generateSequence_Step10_all(*currOptions, // input (changed temporarily internally)
4553  *unifiedCovMatrix, // input
4554  *currRv, // input
4555  useBalancedChains, // input
4556  *unbalancedLinkControl, // input // Round Rock
4557  indexOfFirstWeight, // input // Round Rock
4558  *prevChain, // input // Round Rock
4559  prevExponent, // input
4560  currExponent, // input
4561  prevLogLikelihoodValues, // input
4562  prevLogTargetValues, // input
4563  *balancedLinkControl, // input // Round Rock
4564  currChain, // output
4565  cumulativeRawChainRunTime, // output
4566  cumulativeRawChainRejections, // output
4567  &currLogLikelihoodValues, // output // likelihood is important
4568  &currLogTargetValues); // output);
4569 
4570  //***********************************************************
4571  // Perform checkpoint if necessary
4572  //***********************************************************
4573  stopAtEndOfLevel = currOptions->m_stopAtEnd;
4574  bool performCheckpoint = stopAtEndOfLevel;
4575  if (m_options.m_restartOutput_levelPeriod > 0) {
4576  performCheckpoint = performCheckpoint || ( ((m_currLevel + 1) % m_options.m_restartOutput_levelPeriod) == 0 );
4577  if (currExponent == 1.) {
4578  performCheckpoint = true;
4579  }
4580  }
4581  if (performCheckpoint) {
4582  checkpointML(currExponent, // input
4583  currEta, // input
4584  currChain, // input
4585  currLogLikelihoodValues, // input
4586  currLogTargetValues); // input
4587  }
4588 
4589  //***********************************************************
4590  // Just free some memory
4591  //***********************************************************
4592  {
4593  delete unifiedCovMatrix;
4594 
4595  for (unsigned int i = 0; i < balancedLinkControl->balLinkedChains.size(); ++i) {
4596  queso_require_msg(balancedLinkControl->balLinkedChains[i].initialPosition, "Initial position pointer in step 9 should not be NULL");
4597  delete balancedLinkControl->balLinkedChains[i].initialPosition;
4598  balancedLinkControl->balLinkedChains[i].initialPosition = NULL;
4599  }
4600  balancedLinkControl->balLinkedChains.clear();
4601  }
4602 
4603  //***********************************************************
4604  // Step 11 of 11: filter chain if requested
4605  //***********************************************************
4606  m_currStep = 11;
4607  unsigned int unifiedNumberOfRejections = 0;
4608  if (m_env.inter0Rank() >= 0) {
4609  generateSequence_Step11_inter0(currOptions, // input
4610  currUnifiedRequestedNumSamples, // input
4611  cumulativeRawChainRejections, // input
4612  currChain, // input/output
4613  currLogLikelihoodValues, // input/output
4614  currLogTargetValues, // input/output
4615  unifiedNumberOfRejections); // output
4616  }
4617 
4618  minLogLike = -INFINITY;
4619  maxLogLike = INFINITY;
4620  currLogLikelihoodValues.subMinMaxExtra(0,
4621  currLogLikelihoodValues.subSequenceSize(),
4622  minLogLike,
4623  maxLogLike);
4624  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4625  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4626  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4627  << ", sub minLogLike = " << minLogLike
4628  << ", sub maxLogLike = " << maxLogLike
4629  << std::endl;
4630  }
4631 
4632  m_env.fullComm().Barrier();
4633 
4634  minLogLike = -INFINITY;
4635  maxLogLike = INFINITY;
4636  currLogLikelihoodValues.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4637  0,
4638  currLogLikelihoodValues.subSequenceSize(),
4639  minLogLike,
4640  maxLogLike);
4641  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4642  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4643  << ": at end of level " << m_currLevel+LEVEL_REF_ID
4644  << ", unified minLogLike = " << minLogLike
4645  << ", unified maxLogLike = " << maxLogLike
4646  << std::endl;
4647  }
4648 
4649  //***********************************************************
4650  // Prepare to end current level
4651  //***********************************************************
4652  double levelRunTime = MiscGetEllapsedSeconds(&timevalLevel);
4653  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4654  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4655  << ": ending level " << m_currLevel+LEVEL_REF_ID
4656  << ", having generated " << currChain.subSequenceSize()
4657  << " chain positions"
4658  << ", cumulativeRawChainRunTime = " << cumulativeRawChainRunTime << " seconds"
4659  << ", total level time = " << levelRunTime << " seconds"
4660  << ", cumulativeRawChainRejections = " << cumulativeRawChainRejections
4661  << " (" << 100.*((double) cumulativeRawChainRejections)/((double) currOptions->m_rawChainSize)
4662  << "% at this processor)"
4663  << " (" << 100.*((double) unifiedNumberOfRejections)/((double) currUnifiedRequestedNumSamples)
4664  << "% over all processors)"
4665  << ", stopAtEndOfLevel = " << stopAtEndOfLevel
4666  << std::endl;
4667  }
4668 
4669  if (m_env.inter0Rank() >= 0) {
4670  double minCumulativeRawChainRunTime = 0.;
4671  m_env.inter0Comm().Allreduce((void *) &cumulativeRawChainRunTime, (void *) &minCumulativeRawChainRunTime, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_MIN,
4672  "MLSampling<P_V,P_M>::generateSequence()",
4673  "failed MPI.Allreduce() for min cumulative raw chain run time");
4674 
4675  double maxCumulativeRawChainRunTime = 0.;
4676  m_env.inter0Comm().Allreduce((void *) &cumulativeRawChainRunTime, (void *) &maxCumulativeRawChainRunTime, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_MAX,
4677  "MLSampling<P_V,P_M>::generateSequence()",
4678  "failed MPI.Allreduce() for max cumulative raw chain run time");
4679 
4680  double avgCumulativeRawChainRunTime = 0.;
4681  m_env.inter0Comm().Allreduce((void *) &cumulativeRawChainRunTime, (void *) &avgCumulativeRawChainRunTime, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
4682  "MLSampling<P_V,P_M>::generateSequence()",
4683  "failed MPI.Allreduce() for sum cumulative raw chain run time");
4684  avgCumulativeRawChainRunTime /= ((double) m_env.inter0Comm().NumProc());
4685 
4686  double minLevelRunTime = 0.;
4687  m_env.inter0Comm().Allreduce((void *) &levelRunTime, (void *) &minLevelRunTime, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_MIN,
4688  "MLSampling<P_V,P_M>::generateSequence()",
4689  "failed MPI.Allreduce() for min level run time");
4690 
4691  double maxLevelRunTime = 0.;
4692  m_env.inter0Comm().Allreduce((void *) &levelRunTime, (void *) &maxLevelRunTime, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_MAX,
4693  "MLSampling<P_V,P_M>::generateSequence()",
4694  "failed MPI.Allreduce() for max level run time");
4695 
4696  double avgLevelRunTime = 0.;
4697  m_env.inter0Comm().Allreduce((void *) &levelRunTime, (void *) &avgLevelRunTime, (int) 1, RawValue_MPI_DOUBLE, RawValue_MPI_SUM,
4698  "MLSampling<P_V,P_M>::generateSequence()",
4699  "failed MPI.Allreduce() for sum level run time");
4700  avgLevelRunTime /= ((double) m_env.inter0Comm().NumProc());
4701 
4702  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4703  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4704  << ", level " << m_currLevel+LEVEL_REF_ID
4705  << ": min cumul seconds = " << minCumulativeRawChainRunTime
4706  << ", avg cumul seconds = " << avgCumulativeRawChainRunTime
4707  << ", max cumul seconds = " << maxCumulativeRawChainRunTime
4708  << ", min level seconds = " << minLevelRunTime
4709  << ", avg level seconds = " << avgLevelRunTime
4710  << ", max level seconds = " << maxLevelRunTime
4711  << std::endl;
4712  }
4713  }
4714 
4715  if (currExponent != 1.) delete currOptions;
4716 
4717  struct timeval timevalLevelEnd;
4718  iRC = 0;
4719  iRC = gettimeofday(&timevalLevelEnd, NULL);
4720 
4721  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4722  *m_env.subDisplayFile() << "Getting at the end of level " << m_currLevel+LEVEL_REF_ID
4723  << ", as part of a 'while' on levels"
4724  << ", at " << ctime(&timevalLevelEnd.tv_sec)
4725  << ", after " << timevalLevelEnd.tv_sec - timevalRoutineBegin.tv_sec
4726  << " seconds from entering the routine"
4727  << ", after " << timevalLevelEnd.tv_sec - m_env.timevalBegin().tv_sec
4728  << " seconds from queso environment instatiation"
4729  << std::endl;
4730  }
4731  } // end of level while
4732 
4733 
4734  // m_env.worldRank(),
4735  // "MLSampling<P_V,P_M>::generateSequence()",
4736  // "exponent has not achieved value '1' even after exiting level loop");
4737 
4738  //***********************************************************
4739  // Compute information gain
4740  // ln( \pi(D|M) ) = E[ln( \pi(D|\theta,M) )] - E[ln( \pi(\theta|D,M) / \pi(\theta|M) )]
4741  //***********************************************************
4742  if (m_env.inter0Rank() >= 0) { // KAUST
4743  queso_require_equal_to_msg(m_currLevel, m_logEvidenceFactors.size(), "invalid m_currLevel at the exit of the level loop");
4744  m_logEvidence = 0.;
4745  for (unsigned int i = 0; i < m_logEvidenceFactors.size(); ++i) {
4746  m_logEvidence += m_logEvidenceFactors[i];
4747  }
4748 
4749 #if 1 // prudenci-2012-07-06
4750  m_meanLogLikelihood = currLogLikelihoodValues.unifiedMeanPlain(m_vectorSpace.numOfProcsForStorage() == 1);
4751 #else
4752  m_meanLogLikelihood = currLogLikelihoodValues.unifiedMeanExtra(m_vectorSpace.numOfProcsForStorage() == 1,
4753  0,
4754  currLogLikelihoodValues.subSequenceSize());
4755 #endif
4756 
4757  m_eig = m_meanLogLikelihood - m_logEvidence;
4758 
4759  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4760  *m_env.subDisplayFile() << "In MLSampling<P_V,P_M>::generateSequence()"
4761  << ", log(evidence) = " << m_logEvidence
4762  << ", evidence = " << exp(m_logEvidence)
4763  << ", meanLogLikelihood = " << m_meanLogLikelihood
4764  << ", eig = " << m_eig
4765  << std::endl;
4766  }
4767  }
4768 
4769  m_env.subComm().Bcast((void *) &m_logEvidence, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm'
4770  "MLSampling<P_V,P_M>::generateSequence()",
4771  "failed MPI.Bcast() for m_logEvidence");
4772 
4773  m_env.subComm().Bcast((void *) &m_meanLogLikelihood, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm'
4774  "MLSampling<P_V,P_M>::generateSequence()",
4775  "failed MPI.Bcast() for m_meanLogLikelihood");
4776 
4777  m_env.subComm().Bcast((void *) &m_eig, (int) 1, RawValue_MPI_DOUBLE, 0, // Yes, 'subComm'
4778  "MLSampling<P_V,P_M>::generateSequence()",
4779  "failed MPI.Bcast() for m_eig");
4780 
4781  //***********************************************************
4782  // Prepare to return
4783  //***********************************************************
4784  workingChain.clear();
4785  workingChain.resizeSequence(currChain.subSequenceSize());
4786  P_V auxVec(m_vectorSpace.zeroVector());
4787  for (unsigned int i = 0; i < workingChain.subSequenceSize(); ++i) {
4788  if (m_env.inter0Rank() >= 0) {
4789  currChain.getPositionValues(i,auxVec);
4790  }
4791  workingChain.setPositionValues(i,auxVec);
4792  }
4793 
4794  if (workingLogLikelihoodValues) *workingLogLikelihoodValues = currLogLikelihoodValues;
4795  if (workingLogTargetValues ) *workingLogTargetValues = currLogTargetValues;
4796 
4797  struct timeval timevalRoutineEnd;
4798  iRC = 0;
4799  iRC = gettimeofday(&timevalRoutineEnd, NULL);
4800 
4801  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
4802  *m_env.subDisplayFile() << "Leaving MLSampling<P_V,P_M>::generateSequence()"
4803  << ", at " << ctime(&timevalRoutineEnd.tv_sec)
4804  << ", after " << timevalRoutineEnd.tv_sec - timevalRoutineBegin.tv_sec
4805  << " seconds from entering the routine"
4806  << ", after " << timevalRoutineEnd.tv_sec - m_env.timevalBegin().tv_sec
4807  << " seconds from queso environment instatiation"
4808  << std::endl;
4809  }
4810 
4811  return;
4812 }
4813 
4814 template<class P_V,class P_M>
4815 std::ostream& operator<<(std::ostream& os, const MLSampling<P_V,P_M>& obj)
4816 {
4817  obj.print(os);
4818 
4819  return os;
4820 }
4821 
4822 template <class P_V,class P_M>
4823 double MLSampling<P_V,P_M>::logEvidence() const
4824 {
4825  return m_logEvidence;
4826 }
4827 
4828 template <class P_V,class P_M>
4829 double MLSampling<P_V,P_M>::meanLogLikelihood() const
4830 {
4831  return m_meanLogLikelihood;
4832 }
4833 
4834 template <class P_V,class P_M>
4835 double MLSampling<P_V,P_M>::eig() const
4836 {
4837  return m_eig;
4838 }
4839 
4840 } // End namespace QUESO
4841 
4842 template class QUESO::MLSampling<QUESO::GslVector, QUESO::GslMatrix>;
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
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:52
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::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:3094
RawValue_MPI_IN_PLACE
#define RawValue_MPI_IN_PLACE
Definition: MpiComm.h:44
QUESO::ScalarSequence::setName
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
Definition: ScalarSequence.C:115
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:3127
QUESO::MLSamplingLevelOptions::m_filteredChainDiscardedPortion
double m_filteredChainDiscardedPortion
Initial discarded portion for chain filtering.
Definition: MLSamplingLevelOptions.h:269
QUESO::MLSamplingLevelOptions::m_loadBalanceAlgorithmId
unsigned int m_loadBalanceAlgorithmId
Perform load balancing with chosen algorithm (0 = no balancing).
Definition: MLSamplingLevelOptions.h:173
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::MLSampling::meanLogLikelihood
double meanLogLikelihood() const
Method to calculate the mean of the logarithm of the likelihood.
Definition: MLSampling.C:4829
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:3978
QUESO::ScalarSequence::clear
void clear()
Clears the sequence of scalars.
Definition: ScalarSequence.C:123
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:3228
QUESO::ScalarSequence::unifiedMaxPlain
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
Definition: ScalarSequence.C:380
QUESO::MLSamplingLevelOptions::m_dataOutputFileName
std::string m_dataOutputFileName
Name of generic output file.
Definition: MLSamplingLevelOptions.h:161
QUESO::FiniteDistribution::sample
unsigned int sample() const
Samples.
Definition: FiniteDistribution.C:140
QUESO::ExchangeInfoStruct::numberOfPositions
unsigned int numberOfPositions
Definition: MLSampling.h:71
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:2438
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:614
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:2920
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:874
QUESO::RawType_MPI_Status
MPI_Status RawType_MPI_Status
Definition: MpiComm.h:42
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:3812
QUESO::ScalarSequence::unifiedMeanPlain
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Definition: ScalarSequence.C:405
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
QUESO::MetropolisHastingsSG
A templated class that represents a Metropolis-Hastings generator of samples.
Definition: MetropolisHastingsSG.h:123
QUESO::SequenceOfVectors::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: SequenceOfVectors.C:1923
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:1643
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::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::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:81
RawValue_MPI_MIN
#define RawValue_MPI_MIN
Definition: MpiComm.h:50
ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
#define ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
Definition: MLSampling.h:46
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:72
QUESO::MLSamplingLevelOptions::m_minRejectionRate
double m_minRejectionRate
minimum allowed attempted rejection rate at current level
Definition: MLSamplingLevelOptions.h:188
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_totallyMute
bool m_totallyMute
Whether or not to be totally mute (no printout message).
Definition: MLSamplingLevelOptions.h:201
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:1510
queso_require_less_equal_msg
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:89
QUESO::MHRawChainInfoStruct
A struct that represents a Metropolis-Hastings sample.
Definition: MetropolisHastingsSG.h:58
QUESO::MLSamplingLevelOptions::m_filteredChainGenerate
bool m_filteredChainGenerate
Whether or not to generate filtered chain.
Definition: MLSamplingLevelOptions.h:266
QUESO::MLSamplingLevelOptions::scanOptionsValues
void scanOptionsValues(const MLSamplingLevelOptions *defaultOptions)
Definition: MLSamplingLevelOptions.C:438
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_error_msg
#define queso_error_msg(msg)
Definition: asserts.h:47
QUESO::ScalarSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
Definition: ScalarSequence.C:143
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:3195
QUESO::MLSampling::logEvidence
double logEvidence() const
Method to calculate the logarithm of the evidence.
Definition: MLSampling.C:4823
QUESO::MLSamplingLevelOptions::m_minAcceptableEta
double m_minAcceptableEta
minimum acceptable eta,
Definition: MLSamplingLevelOptions.h:198
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::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:2111
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:2497
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::MLSampling::eig
double eig() const
Calculates the expected information gain value, EIG.
Definition: MLSampling.C:4835
QUESO::MLSamplingLevelOptions::m_rawChainDataOutputFileName
std::string m_rawChainDataOutputFileName
Name of output file for raw chain.
Definition: MLSamplingLevelOptions.h:247
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:2004
QUESO::MLSamplingLevelOptions::m_rawChainSize
unsigned int m_rawChainSize
Size of raw chain.
Definition: MLSamplingLevelOptions.h:232
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:2484
QUESO::MLSampling::~MLSampling
~MLSampling()
Destructor.
Definition: MLSampling.C:4011
QUESO::MLSamplingLevelOptions::m_maxRejectionRate
double m_maxRejectionRate
maximum allowed attempted rejection rate at current level.
Definition: MLSamplingLevelOptions.h:191
QUESO::MLSamplingLevelOptions::m_dataOutputAllowedSet
std::set< unsigned int > m_dataOutputAllowedSet
subEnvs that will write to generic output file.
Definition: MLSamplingLevelOptions.h:167
QUESO::UnbalancedLinkedChainsPerNodeStruct::unbLinkedChains
std::vector< UnbalancedLinkedChainControlStruct > unbLinkedChains
Definition: MLSampling.h:101
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::MLSamplingLevelOptions::m_parameterDisabledSet
std::set< unsigned int > m_parameterDisabledSet
Definition: MLSamplingLevelOptions.h:219
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
QUESO::BaseVectorSequence::clear
void clear()
Reset the values and the size of the sequence of vectors.
Definition: VectorSequence.C:133
QUESO::MLSamplingLevelOptions::m_covRejectionRate
double m_covRejectionRate
c.o.v. for judging attempted rejection rate at current level.
Definition: MLSamplingLevelOptions.h:194
QUESO::MLSamplingLevelOptions::m_rawChainDataOutputFileType
std::string m_rawChainDataOutputFileType
Type of output file for raw chain.
Definition: MLSamplingLevelOptions.h:250
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::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:282
queso_require_msg
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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:1539
UQ_MH_SG_FILENAME_FOR_NO_FILE
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
Definition: MetropolisHastingsSGOptions.h:36
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:274
QUESO::MLSamplingLevelOptions::m_amAdaptInterval
unsigned int m_amAdaptInterval
&#39;am&#39; adaptation interval.
Definition: MLSamplingLevelOptions.h:322
queso_require_less_msg
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:87
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2026
QUESO::ExchangeInfoStruct::originalIndexOfInitialPosition
unsigned int originalIndexOfInitialPosition
Definition: MLSampling.h:69
QUESO::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2632
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:2291
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:2616
QUESO::GenericVectorRV::setPdf
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
Definition: GenericVectorRV.C:90
UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:89
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:3698
QUESO::MLSamplingLevelOptions::m_filteredChainDataOutputFileType
std::string m_filteredChainDataOutputFileType
Type of output file for filtered chain.
Definition: MLSamplingLevelOptions.h:278
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:76
RawValue_MPI_MAX
#define RawValue_MPI_MAX
Definition: MpiComm.h:51
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:4023
QUESO::MLSamplingLevelOptions::m_scaleCovMatrix
bool m_scaleCovMatrix
Whether or not scale proposal covariance matrix.
Definition: MLSamplingLevelOptions.h:185
QUESO::MLSamplingLevelOptions::m_filteredChainLag
unsigned int m_filteredChainLag
Spacing for chain filtering.
Definition: MLSamplingLevelOptions.h:272
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::MLSamplingLevelOptions::m_filteredChainDataOutputFileName
std::string m_filteredChainDataOutputFileName
Name of output file for filtered chain.
Definition: MLSamplingLevelOptions.h:275
QUESO::SequenceOfVectors::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
Definition: SequenceOfVectors.C:103
QUESO::MLSampling::justBalance_proc0
void justBalance_proc0(const MLSamplingLevelOptions *currOptions, std::vector< ExchangeInfoStruct > &exchangeStdVec)
Definition: MLSampling.C:1495
QUESO::MLSampling
A templated class that represents a Multilevel generator of samples.
Definition: MLSampling.h:122
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::MetropolisHastingsSG::getRawChainInfo
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
Definition: MetropolisHastingsSG.C:1430
QUESO::BaseVectorSequence::resizeSequence
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
QUESO::MLSamplingLevelOptions::m_initialValuesOfDisabledParameters
std::vector< double > m_initialValuesOfDisabledParameters
Definition: MLSamplingLevelOptions.h:221
QUESO::MLSamplingLevelOptions::m_loadBalanceTreshold
double m_loadBalanceTreshold
Perform load balancing if load unbalancing ratio &gt; threshold.
Definition: MLSamplingLevelOptions.h:176
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::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:3014
queso_require_greater_equal_msg
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
RawValue_MPI_CHAR
#define RawValue_MPI_CHAR
Definition: MpiComm.h:46
QUESO::SequenceOfVectors::resizeSequence
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
Definition: SequenceOfVectors.C:110
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::MLSampling::m_env
const BaseEnvironment & m_env
Queso enviroment.
Definition: MLSampling.h:470
QUESO::MLSamplingLevelOptions::m_stopAtEnd
bool m_stopAtEnd
Stop at end of such level.
Definition: MLSamplingLevelOptions.h:158
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::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
QUESO::FiniteDistribution
A templated class for a finite distribution.
Definition: FiniteDistribution.h:45
QUESO::ScalarSequence::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:136
QUESO::MiscCheckForSameValueInAllNodes
bool MiscCheckForSameValueInAllNodes(T &inputValue, double acceptableTreshold, const MpiComm &comm, const char *whereString)
Definition: Miscellaneous.C:479
RawValue_MPI_DOUBLE
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:48
LEVEL_REF_ID
#define LEVEL_REF_ID
Definition: MLSamplingLevelOptions.h:28
QUESO::MLSamplingLevelOptions::m_drMaxNumExtraStages
unsigned int m_drMaxNumExtraStages
&#39;dr&#39; maximum number of extra stages.
Definition: MLSamplingLevelOptions.h:306
QUESO::MLSamplingLevelOptions::m_minEffectiveSizeRatio
double m_minEffectiveSizeRatio
Minimum allowed effective size ratio wrt previous level.
Definition: MLSamplingLevelOptions.h:179
RawValue_MPI_UNSIGNED
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:49
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2429
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
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:2833
QUESO::SequenceOfVectors::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Definition: SequenceOfVectors.C:1377
QUESO::MLSamplingLevelOptions::m_maxEffectiveSizeRatio
double m_maxEffectiveSizeRatio
Maximum allowed effective size ratio wrt previous level.
Definition: MLSamplingLevelOptions.h:182
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.
Generated on Thu Jun 11 2015 13:52:32 for queso-0.53.0 by  doxygen 1.8.5