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