queso-0.56.0
MLSampling.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <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::MLSamplingLevelOptions::m_filteredChainDiscardedPortion
double m_filteredChainDiscardedPortion
Initial discarded portion for chain filtering.
Definition: MLSamplingLevelOptions.h:274
QUESO::MLSampling::logEvidence
double logEvidence() const
Method to calculate the logarithm of the evidence.
Definition: MLSampling.C:4830
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::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::MLSamplingLevelOptions::m_maxEffectiveSizeRatio
double m_maxEffectiveSizeRatio
Maximum allowed effective size ratio wrt previous level.
Definition: MLSamplingLevelOptions.h:184
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:89
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_filteredChainGenerate
bool m_filteredChainGenerate
Whether or not to generate filtered chain.
Definition: MLSamplingLevelOptions.h:271
QUESO::BaseVectorSequence::setName
void setName(const std::string &newName)
Changes the name of the sequence of vectors.
Definition: VectorSequence.C:125
QUESO::BaseVectorSequence::subSequenceSize
virtual unsigned int subSequenceSize() const =0
Size of the sub-sequence of vectors. See template specialization.
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:282
QUESO::RawType_MPI_Status
int RawType_MPI_Status
Definition: MpiComm.h:62
QUESO::MiscCheckForSameValueInAllNodes
bool MiscCheckForSameValueInAllNodes(T &inputValue, double acceptableTreshold, const MpiComm &comm, const char *whereString)
Definition: Miscellaneous.C:479
RawValue_MPI_MIN
#define RawValue_MPI_MIN
Definition: MpiComm.h:69
QUESO::ExchangeInfoStruct::originalIndexOfInitialPosition
unsigned int originalIndexOfInitialPosition
Definition: MLSampling.h:69
QUESO::BaseEnvironment::subDisplayFile
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
queso_require_greater_equal_msg
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
QUESO::MLSamplingLevelOptions::m_covRejectionRate
double m_covRejectionRate
c.o.v. for judging attempted rejection rate at current level.
Definition: MLSamplingLevelOptions.h:196
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:1530
QUESO::ScalarSequence::unifiedMeanPlain
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Definition: ScalarSequence.C:407
QUESO::MLSamplingLevelOptions::m_minEffectiveSizeRatio
double m_minEffectiveSizeRatio
Minimum allowed effective size ratio wrt previous level.
Definition: MLSamplingLevelOptions.h:181
QUESO::FiniteDistribution
A templated class for a finite distribution.
Definition: FiniteDistribution.h:45
QUESO::MLSampling::meanLogLikelihood
double meanLogLikelihood() const
Method to calculate the mean of the logarithm of the likelihood.
Definition: MLSampling.C:4836
QUESO::MLSamplingLevelOptions::m_dataOutputFileName
std::string m_dataOutputFileName
Name of generic output file.
Definition: MLSamplingLevelOptions.h:163
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:2897
QUESO::MLSamplingLevelOptions::m_minRejectionRate
double m_minRejectionRate
minimum allowed attempted rejection rate at current level
Definition: MLSamplingLevelOptions.h:190
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::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_initialPositionUsePreviousLevelLikelihood
bool m_initialPositionUsePreviousLevelLikelihood
Use previous level likelihood for initial chain position instead of re-computing it from target pdf...
Definition: MLSamplingLevelOptions.h:218
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::ScalarSequence::clear
void clear()
Clears the sequence of scalars.
Definition: ScalarSequence.C:123
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:2488
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::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::SequenceOfVectors::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: SequenceOfVectors.C:2042
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::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::BaseVectorSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize() const
Calculates the size of the unified sequence of vectors.
Definition: VectorSequence.C:70
QUESO::MetropolisHastingsSG::getRawChainInfo
void getRawChainInfo(MHRawChainInfoStruct &info) const
Gets information from the raw chain.
Definition: MetropolisHastingsSG.C:1313
QUESO::MLSamplingLevelOptions::m_rawChainSize
unsigned int m_rawChainSize
Size of raw chain.
Definition: MLSamplingLevelOptions.h:234
QUESO::MLSamplingLevelOptions::m_scaleCovMatrix
bool m_scaleCovMatrix
Whether or not scale proposal covariance matrix.
Definition: MLSamplingLevelOptions.h:187
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::sampleIndexes_proc0
void sampleIndexes_proc0(unsigned int unifiedRequestedNumSamples, const std::vector< double > &unifiedWeightStdVectorAtProc0Only, std::vector< unsigned int > &unifiedIndexCountersAtProc0Only)
Definition: MLSampling.C:95
QUESO::MLSamplingLevelOptions::m_filteredChainLag
unsigned int m_filteredChainLag
Spacing for chain filtering.
Definition: MLSamplingLevelOptions.h:277
QUESO::SequenceOfVectors::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
Definition: SequenceOfVectors.C:103
QUESO::BaseVectorSequence::resizeSequence
virtual void resizeSequence(unsigned int newSubSequenceSize)=0
Resize the sequence. See template specialization.
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
queso_require_less_msg
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
QUESO::SequenceOfVectors::resizeSequence
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
Definition: SequenceOfVectors.C:110
QUESO::GenericVectorRV::setPdf
void setPdf(BaseJointPdf< V, M > &pdf)
Sets the PDF of this vector RV to pdf.
Definition: GenericVectorRV.C:90
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::ScalarSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
Definition: ScalarSequence.C:143
QUESO::ScalarSequence::setName
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
Definition: ScalarSequence.C:115
RawValue_MPI_MAX
#define RawValue_MPI_MAX
Definition: MpiComm.h:70
QUESO::MLSamplingLevelOptions::m_loadBalanceTreshold
double m_loadBalanceTreshold
Perform load balancing if load unbalancing ratio &gt; threshold.
Definition: MLSamplingLevelOptions.h:178
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:76
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2420
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::MetropolisHastingsSG
A templated class that represents a Metropolis-Hastings generator of samples.
Definition: MetropolisHastingsSG.h:126
QUESO::FiniteDistribution::sample
unsigned int sample() const
Samples.
Definition: FiniteDistribution.C:140
QUESO::MLSampling::m_env
const BaseEnvironment & m_env
Queso enviroment.
Definition: MLSampling.h:470
RawValue_MPI_DOUBLE
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
RawValue_MPI_UNSIGNED
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:68
QUESO::MLSamplingLevelOptions::m_drMaxNumExtraStages
unsigned int m_drMaxNumExtraStages
&#39;dr&#39; maximum number of extra stages.
Definition: MLSamplingLevelOptions.h:311
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::MLSamplingLevelOptions::m_minAcceptableEta
double m_minAcceptableEta
minimum acceptable eta,
Definition: MLSamplingLevelOptions.h:200
QUESO::MLSamplingLevelOptions::m_parameterDisabledSet
std::set< unsigned int > m_parameterDisabledSet
Definition: MLSamplingLevelOptions.h:221
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::MLSamplingLevelOptions::m_totallyMute
bool m_totallyMute
Whether or not to be totally mute (no printout message).
Definition: MLSamplingLevelOptions.h:203
RawValue_MPI_CHAR
#define RawValue_MPI_CHAR
Definition: MpiComm.h:65
QUESO::MLSamplingLevelOptions::m_rawChainDataOutputFileName
std::string m_rawChainDataOutputFileName
Name of output file for raw chain.
Definition: MLSamplingLevelOptions.h:249
QUESO::UnbalancedLinkedChainsPerNodeStruct::unbLinkedChains
std::vector< UnbalancedLinkedChainControlStruct > unbLinkedChains
Definition: MLSampling.h:101
QUESO::MLSamplingLevelOptions::m_maxRejectionRate
double m_maxRejectionRate
maximum allowed attempted rejection rate at current level.
Definition: MLSamplingLevelOptions.h:193
UQ_MH_SG_FILENAME_FOR_NO_FILE
#define UQ_MH_SG_FILENAME_FOR_NO_FILE
Definition: MetropolisHastingsSGOptions.h:39
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::MLSamplingLevelOptions
This class provides options for each level of the Multilevel sequence generator if no input file is a...
Definition: MLSamplingLevelOptions.h:119
QUESO::MLSamplingLevelOptions::m_stopAtEnd
bool m_stopAtEnd
Stop at end of such level.
Definition: MLSamplingLevelOptions.h:160
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::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
reason
if the work is an executable linked with the with the complete machine readable work that uses the as object code and or source so that the user can modify the Library and then relink to produce a modified executable containing the modified rather than copying library functions into the if the user installs as long as the modified version is interface compatible with the version that the work was made with c Accompany the work with a written valid for at least three to give the same user the materials specified in for a charge no more than the cost of performing this distribution d If distribution of the work is made by offering access to copy from a designated offer equivalent access to copy the above specified materials from the same place e Verify that the user has already received a copy of these materials or that you have already sent this user a copy For an the required form of the work that uses the Library must include any data and utility programs needed for reproducing the executable from it as a special the materials to be distributed need not include anything that is normally and so on of the operating system on which the executable unless that component itself accompanies the executable It may happen that this requirement contradicts the license restrictions of other proprietary libraries that do not normally accompany the operating system Such a contradiction means you cannot use both them and the Library together in an executable that you distribute You may place library facilities that are a work based on the Library side by side in a single library together with other library facilities not covered by this and distribute such a combined provided that the separate distribution of the work based on the Library and of the other library facilities is otherwise and provided that you do these two uncombined with any other library facilities This must be distributed under the terms of the Sections above b Give prominent notice with the combined library of the fact that part of it is a work based on the and explaining where to find the accompanying uncombined form of the same work You may not link or distribute the Library except as expressly provided under this License Any attempt otherwise to link or distribute the Library is and will automatically terminate your rights under this License parties who have received or from you under this License will not have their licenses terminated so long as such parties remain in full compliance You are not required to accept this since you have not signed it nothing else grants you permission to modify or distribute the Library or its derivative works These actions are prohibited by law if you do not accept this License by modifying or distributing the you indicate your acceptance of this License to do and all its terms and conditions for distributing or modifying the Library or works based on it Each time you redistribute the the recipient automatically receives a license from the original licensor to link with or modify the Library subject to these terms and conditions You may not impose any further restrictions on the recipients exercise of the rights granted herein You are not responsible for enforcing compliance by third parties with this License as a consequence of a court judgment or allegation of patent infringement or for any other reason(not limited to patent issues)
QUESO::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::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2672
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:85
QUESO::MHRawChainInfoStruct
A struct that represents a Metropolis-Hastings sample.
Definition: MetropolisHastingsSG.h:61
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_dataOutputAllowedSet
std::set< unsigned int > m_dataOutputAllowedSet
subEnvs that will write to generic output file.
Definition: MLSamplingLevelOptions.h:169
QUESO::SequenceOfVectors::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Definition: SequenceOfVectors.C:1493
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
LEVEL_REF_ID
#define LEVEL_REF_ID
Definition: MLSamplingLevelOptions.h:28
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_require_less_equal_msg
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:77
QUESO::ExchangeInfoStruct::numberOfPositions
unsigned int numberOfPositions
Definition: MLSampling.h:71
QUESO::MLSamplingLevelOptions::m_filteredChainDataOutputFileType
std::string m_filteredChainDataOutputFileType
Type of output file for filtered chain.
Definition: MLSamplingLevelOptions.h:283
QUESO::BaseVectorSequence::clear
void clear()
Reset the values and the size of the sequence of vectors.
Definition: VectorSequence.C:133
QUESO::MLSamplingLevelOptions::m_filteredChainDataOutputFileName
std::string m_filteredChainDataOutputFileName
Name of output file for filtered chain.
Definition: MLSamplingLevelOptions.h:280
UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
queso_require_msg
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
QUESO::MLSamplingLevelOptions::m_rawChainDataOutputFileType
std::string m_rawChainDataOutputFileType
Type of output file for raw chain.
Definition: MLSamplingLevelOptions.h:255
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
QUESO::BalancedLinkedChainsPerNodeStruct::balLinkedChains
std::vector< BalancedLinkedChainControlStruct< P_V > > balLinkedChains
Definition: MLSampling.h:88
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::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::MLSamplingLevelOptions::scanOptionsValues
void scanOptionsValues(const MLSamplingLevelOptions *defaultOptions)
Definition: MLSamplingLevelOptions.C:585
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::MLSampling::~MLSampling
~MLSampling()
Destructor.
Definition: MLSampling.C:4018
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:1501
QUESO::ScalarSequence::unifiedMeanExtra
T unifiedMeanExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos) const
Finds the mean value of the unified sequence of numPos positions starting at position initialPos...
Definition: ScalarSequence.C:876
queso_error_msg
#define queso_error_msg(msg)
Definition: asserts.h:47
QUESO::MLSamplingLevelOptions::m_amAdaptInterval
unsigned int m_amAdaptInterval
&#39;am&#39; adaptation interval.
Definition: MLSamplingLevelOptions.h:327
QUESO::MLSampling::justBalance_proc0
void justBalance_proc0(const MLSamplingLevelOptions *currOptions, std::vector< ExchangeInfoStruct > &exchangeStdVec)
Definition: MLSampling.C:1500
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::MLSamplingLevelOptions::m_loadBalanceAlgorithmId
unsigned int m_loadBalanceAlgorithmId
Perform load balancing with chosen algorithm (0 = no balancing).
Definition: MLSamplingLevelOptions.h:175
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:860
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::ScalarSequence::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:136
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::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::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::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::MLSamplingLevelOptions::m_initialValuesOfDisabledParameters
std::vector< double > m_initialValuesOfDisabledParameters
Definition: MLSamplingLevelOptions.h:223
QUESO::matrixProduct
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2036
ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
#define ML_CHECKPOINT_FIXED_AMOUNT_OF_DATA
Definition: MLSampling.h:46
Generated on Tue Nov 29 2016 10:53:11 for queso-0.56.0 by  doxygen 1.8.5