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