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