queso-0.57.0
ScalarFunctionSynchronizer.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/Environment.h>
26 #include <queso/ScalarFunction.h>
27 #include <queso/ScalarFunctionSynchronizer.h>
28 #include <queso/BayesianJointPdf.h>
29 #include <queso/JointPdf.h>
30 #include <queso/GslVector.h>
31 #include <queso/GslMatrix.h>
32 
33 namespace QUESO {
34 
35 // Default constructor
36 template <class V, class M>
37 ScalarFunctionSynchronizer<V,M>::ScalarFunctionSynchronizer(
38  const BaseScalarFunction<V,M>& inputFunction,
39  const V& auxVec)
40  : m_env(inputFunction.domainSet().env()),
41  m_scalarFunction(inputFunction),
42  m_bayesianJointPdfPtr(dynamic_cast<const BayesianJointPdf<V,M>* >(&m_scalarFunction)),
43  m_auxVec(auxVec)
44 {
45 }
46 
47 // Destructor
48 template <class V, class M>
49 ScalarFunctionSynchronizer<V,M>::~ScalarFunctionSynchronizer()
50 {
51 }
52 
53 // Math methods
54 template<class V,class M>
55 const VectorSet<V,M>& ScalarFunctionSynchronizer<V,M>::domainSet() const
56 {
57  return m_scalarFunction.domainSet();
58 }
59 
60 // Sync methods
61 template <class V,class M>
62 double ScalarFunctionSynchronizer<V,M>::callFunction(const V* vecValues,
63  const V* vecDirection,
64  V* gradVector,
65  M* hessianMatrix,
66  V* hessianEffect,
67  double* extraOutput1,
68  double* extraOutput2) const
69 {
70  double result = 0.;
71 
72  if ((m_env.numSubEnvironments() < (unsigned int) m_env.fullComm().NumProc()) &&
73  (m_auxVec.numOfProcsForStorage() == 1 )) {
74  bool stayInRoutine = true;
75  do {
76  const V* internalValues = NULL;
77  const V* internalDirection = NULL;
78  V* internalGrad = NULL;
79  M* internalHessian = NULL;
80  V* internalEffect = NULL;
81 
83  // Broadcast 1 of 3
85  // bufferChar[0] = '0' or '1' (vecValues is NULL or not)
86  // bufferChar[1] = '0' or '1' (vecDirection is NULL or not)
87  // bufferChar[2] = '0' or '1' (gradVector is NULL or not)
88  // bufferChar[3] = '0' or '1' (hessianMatrix is NULL or not)
89  // bufferChar[4] = '0' or '1' (hessianEffect is NULL or not)
90  std::vector<char> bufferChar(5,'0');
91 
92  if (m_env.subRank() == 0) {
93  internalValues = vecValues;
94  internalDirection = vecDirection;
95  internalGrad = gradVector;
96  internalHessian = hessianMatrix;
97  internalEffect = hessianEffect;
98 
99  if (internalValues != NULL) bufferChar[0] = '1';
100  if (internalDirection != NULL) bufferChar[1] = '1';
101  if (internalGrad != NULL) bufferChar[2] = '1';
102  if (internalHessian != NULL) bufferChar[3] = '1';
103  if (internalEffect != NULL) bufferChar[4] = '1';
104  }
105 
106  m_env.subComm().syncPrintDebugMsg("In ScalarFunctionSynchronizer<V,M>::callFunction(), just before char Bcast()",3,3000000);
107  //if (m_env.subId() != 0) while (true) sleep(1);
108 
109  int count = (int) bufferChar.size();
110  m_env.subComm().Bcast((void *) &bufferChar[0], count, RawValue_MPI_CHAR, 0,
111  "ScalarFunctionSynchronizer<V,M>::callFunction()",
112  "failed broadcast 1 of 3");
113 
114  m_env.subComm().syncPrintDebugMsg("In ScalarFunctionSynchronizer<V,M>::callFunction(), just after char Bcast()",3,3000000);
115  //std::cout << "char contents = " << bufferChar[0] << " " << bufferChar[1] << " " << bufferChar[2] << " " << bufferChar[3] << " " << bufferChar[4]
116  // << std::endl;
117 
118  if (bufferChar[0] == '1') {
120  // Broadcast 2 of 3
122 
123  // bufferDouble[0...] = contents for (eventual) vecValues
124  std::vector<double> bufferDouble(m_auxVec.sizeLocal(),0.);
125 
126  if (m_env.subRank() == 0) {
127  for (unsigned int i = 0; i < internalValues->sizeLocal(); ++i) {
128  bufferDouble[i] = (*internalValues)[i];
129  }
130  }
131 
132  //m_env.fullComm().Barrier();
133  //for (int i = 0; i < m_env.fullComm().NumProc(); ++i) {
134  // if (i == m_env.fullRank()) {
135  // std::cout << " In ScalarFunctionSynchronizer<V,M>::callFunction(), just before double Bcast()"
136  // << ": fullRank " << m_env.fullRank()
137  // << ", subEnvironment " << m_env.subId()
138  // << ", subRank " << m_env.subRank()
139  // << ": buffer related to first double Bcast() is ready to be broadcasted"
140  // << " and has size " << bufferDouble.size()
141  // << std::endl;
142  // if (m_env.subRank() == 0) {
143  // std::cout << "Buffer contents are";
144  // for (unsigned int i = 0; i < bufferDouble.size(); ++i) {
145  // std::cout << " " << bufferDouble[i];
146  // }
147  // std::cout << std::endl;
148  // }
149  // }
150  // m_env.fullComm().Barrier();
151  //}
152  //if (m_env.fullRank() == 0) std::cout << "Sleeping 3 seconds..."
153  // << std::endl;
154  //sleep(3);
155 
156  count = (int) bufferDouble.size();
157  m_env.subComm().Bcast((void *) &bufferDouble[0], count, RawValue_MPI_DOUBLE, 0,
158  "ScalarFunctionSynchronizer<V,M>::callFunction()",
159  "failed broadcast 2 of 3");
160 
161  if (m_env.subRank() != 0) {
162  V tmpVec(m_auxVec);
163  for (unsigned int i = 0; i < tmpVec.sizeLocal(); ++i) {
164  tmpVec[i] = bufferDouble[i];
165  }
166  internalValues = new V(tmpVec);
167  //if (vecValues) *vecValues = tmpVec; // prudencio 2010-08-01
168  }
169 
170  if (bufferChar[1] == '1') {
172  // Broadcast 3 of 3
174  // bufferDouble[0...] = contents for (eventual) vecDirection
175 
176  if (m_env.subRank() == 0) {
177  for (unsigned int i = 0; i < internalDirection->sizeLocal(); ++i) {
178  bufferDouble[i] = (*internalDirection)[i];
179  }
180  }
181 
182  count = (int) bufferDouble.size();
183  m_env.subComm().Bcast((void *) &bufferDouble[0], count, RawValue_MPI_DOUBLE, 0,
184  "ScalarFunctionSynchronizer<V,M>::callFunction()",
185  "failed broadcast 3 of 3");
186 
187  if (m_env.subRank() != 0) {
188  V tmpVec(m_auxVec);
189  for (unsigned int i = 0; i < tmpVec.sizeLocal(); ++i) {
190  tmpVec[i] = bufferDouble[i];
191  }
192  internalDirection = new V(tmpVec);
193  }
194  }
195 
197  // All processors now call 'scalarFunction()'
199  if (m_env.subRank() != 0) {
200  if (bufferChar[2] == '1') internalGrad = new V(m_auxVec);
201  if (bufferChar[3] == '1') internalHessian = new M(m_auxVec);
202  if (bufferChar[4] == '1') internalEffect = new V(m_auxVec);
203  }
204 
205  m_env.subComm().syncPrintDebugMsg("In ScalarFunctionSynchronizer<V,M>::callFunction(), just before actual lnValue()",3,3000000);
206  m_env.subComm().Barrier();
207  result = m_scalarFunction.lnValue(*internalValues, // input
208  internalDirection, // input
209  internalGrad, // output
210  internalHessian, // output
211  internalEffect); // output
212  if (extraOutput1) {
213  if (m_bayesianJointPdfPtr) {
214  *extraOutput1 = m_bayesianJointPdfPtr->lastComputedLogPrior();
215  }
216  }
217  if (extraOutput2) {
218  if (m_bayesianJointPdfPtr) {
219  *extraOutput2 = m_bayesianJointPdfPtr->lastComputedLogLikelihood();
220  }
221  }
222  } // if (bufferChar[0] == '1')
223 
225  // Prepare to exit routine or to stay in it
227  if (m_env.subRank() == 0) {
228  stayInRoutine = false; // Always for processor 0
229  }
230  else {
231  if (internalValues != NULL) delete internalValues;
232  if (internalDirection != NULL) delete internalDirection;
233  if (internalGrad != NULL) delete internalGrad;
234  if (internalHessian != NULL) delete internalHessian;
235  if (internalEffect != NULL) delete internalEffect;
236 
237  stayInRoutine = (vecValues == NULL) && (bufferChar[0] == '1');
238  //if (!stayInRoutine) std::cout << "Fullrank " << m_env.fullRank() << " is leaving scalarFunctionSync()" << std::endl;
239  }
240  } while (stayInRoutine);
241  }
242  else {
243  queso_require_msg(vecValues, "vecValues should not be NULL");
244 
245  m_env.subComm().Barrier();
246  result = m_scalarFunction.lnValue(*vecValues,
247  vecDirection,
248  gradVector,
249  hessianMatrix,
250  hessianEffect);
251  if (extraOutput1) {
252  if (m_bayesianJointPdfPtr) {
253  *extraOutput1 = m_bayesianJointPdfPtr->lastComputedLogPrior();
254  }
255  }
256  if (extraOutput2) {
257  if (m_bayesianJointPdfPtr) {
258  *extraOutput2 = m_bayesianJointPdfPtr->lastComputedLogLikelihood();
259  }
260  }
261  }
262 
263  return result;
264 }
265 
266 // Sync methods
267 template <class V,class M>
268 double ScalarFunctionSynchronizer<V,M>::callFunction(const V* vecValues,
269  double* extraOutput1,
270  double* extraOutput2) const
271 {
272  double result = 0.;
273 
274  if ((m_env.numSubEnvironments() < (unsigned int) m_env.fullComm().NumProc()) &&
275  (m_auxVec.numOfProcsForStorage() == 1 )) {
276  bool stayInRoutine = true;
277  do {
278  const V* internalValues = NULL;
279  const V* internalDirection = NULL;
280  V* internalGrad = NULL;
281  M* internalHessian = NULL;
282  V* internalEffect = NULL;
283 
285  // Broadcast 1 of 3
287  // bufferChar[0] = '0' or '1' (vecValues is NULL or not)
288  // bufferChar[1] = '0' or '1' (vecDirection is NULL or not)
289  // bufferChar[2] = '0' or '1' (gradVector is NULL or not)
290  // bufferChar[3] = '0' or '1' (hessianMatrix is NULL or not)
291  // bufferChar[4] = '0' or '1' (hessianEffect is NULL or not)
292  std::vector<char> bufferChar(5,'0');
293 
294  if (m_env.subRank() == 0) {
295  internalValues = vecValues;
296 
297  if (internalValues != NULL) bufferChar[0] = '1';
298  if (internalDirection != NULL) bufferChar[1] = '1';
299  if (internalGrad != NULL) bufferChar[2] = '1';
300  if (internalHessian != NULL) bufferChar[3] = '1';
301  if (internalEffect != NULL) bufferChar[4] = '1';
302  }
303 
304  m_env.subComm().syncPrintDebugMsg("In ScalarFunctionSynchronizer<V,M>::callFunction(), just before char Bcast()",3,3000000);
305  //if (m_env.subId() != 0) while (true) sleep(1);
306 
307  int count = (int) bufferChar.size();
308  m_env.subComm().Bcast((void *) &bufferChar[0], count, RawValue_MPI_CHAR, 0,
309  "ScalarFunctionSynchronizer<V,M>::callFunction()",
310  "failed broadcast 1 of 3");
311 
312  m_env.subComm().syncPrintDebugMsg("In ScalarFunctionSynchronizer<V,M>::callFunction(), just after char Bcast()",3,3000000);
313  //std::cout << "char contents = " << bufferChar[0] << " " << bufferChar[1] << " " << bufferChar[2] << " " << bufferChar[3] << " " << bufferChar[4]
314  // << std::endl;
315 
316  if (bufferChar[0] == '1') {
318  // Broadcast 2 of 3
320 
321  // bufferDouble[0...] = contents for (eventual) vecValues
322  std::vector<double> bufferDouble(m_auxVec.sizeLocal(),0.);
323 
324  if (m_env.subRank() == 0) {
325  for (unsigned int i = 0; i < internalValues->sizeLocal(); ++i) {
326  bufferDouble[i] = (*internalValues)[i];
327  }
328  }
329 
330  //m_env.fullComm().Barrier();
331  //for (int i = 0; i < m_env.fullComm().NumProc(); ++i) {
332  // if (i == m_env.fullRank()) {
333  // std::cout << " In ScalarFunctionSynchronizer<V,M>::callFunction(), just before double Bcast()"
334  // << ": fullRank " << m_env.fullRank()
335  // << ", subEnvironment " << m_env.subId()
336  // << ", subRank " << m_env.subRank()
337  // << ": buffer related to first double Bcast() is ready to be broadcasted"
338  // << " and has size " << bufferDouble.size()
339  // << std::endl;
340  // if (m_env.subRank() == 0) {
341  // std::cout << "Buffer contents are";
342  // for (unsigned int i = 0; i < bufferDouble.size(); ++i) {
343  // std::cout << " " << bufferDouble[i];
344  // }
345  // std::cout << std::endl;
346  // }
347  // }
348  // m_env.fullComm().Barrier();
349  //}
350  //if (m_env.fullRank() == 0) std::cout << "Sleeping 3 seconds..."
351  // << std::endl;
352  //sleep(3);
353 
354  count = (int) bufferDouble.size();
355  m_env.subComm().Bcast((void *) &bufferDouble[0], count, RawValue_MPI_DOUBLE, 0,
356  "ScalarFunctionSynchronizer<V,M>::callFunction()",
357  "failed broadcast 2 of 3");
358 
359  if (m_env.subRank() != 0) {
360  V tmpVec(m_auxVec);
361  for (unsigned int i = 0; i < tmpVec.sizeLocal(); ++i) {
362  tmpVec[i] = bufferDouble[i];
363  }
364  internalValues = new V(tmpVec);
365  //if (vecValues) *vecValues = tmpVec; // prudencio 2010-08-01
366  }
367 
368  if (bufferChar[1] == '1') {
370  // Broadcast 3 of 3
372  // bufferDouble[0...] = contents for (eventual) vecDirection
373 
374  if (m_env.subRank() == 0) {
375  for (unsigned int i = 0; i < internalDirection->sizeLocal(); ++i) {
376  bufferDouble[i] = (*internalDirection)[i];
377  }
378  }
379 
380  count = (int) bufferDouble.size();
381  m_env.subComm().Bcast((void *) &bufferDouble[0], count, RawValue_MPI_DOUBLE, 0,
382  "ScalarFunctionSynchronizer<V,M>::callFunction()",
383  "failed broadcast 3 of 3");
384 
385  if (m_env.subRank() != 0) {
386  V tmpVec(m_auxVec);
387  for (unsigned int i = 0; i < tmpVec.sizeLocal(); ++i) {
388  tmpVec[i] = bufferDouble[i];
389  }
390  internalDirection = new V(tmpVec);
391  }
392  }
393 
395  // All processors now call 'scalarFunction()'
397  if (m_env.subRank() != 0) {
398  if (bufferChar[2] == '1') internalGrad = new V(m_auxVec);
399  if (bufferChar[3] == '1') internalHessian = new M(m_auxVec);
400  if (bufferChar[4] == '1') internalEffect = new V(m_auxVec);
401  }
402 
403  m_env.subComm().syncPrintDebugMsg("In ScalarFunctionSynchronizer<V,M>::callFunction(), just before actual lnValue()",3,3000000);
404  m_env.subComm().Barrier();
405  result = m_scalarFunction.lnValue(*internalValues);
406  if (extraOutput1) {
407  if (m_bayesianJointPdfPtr) {
408  *extraOutput1 = m_bayesianJointPdfPtr->lastComputedLogPrior();
409  }
410  }
411  if (extraOutput2) {
412  if (m_bayesianJointPdfPtr) {
413  *extraOutput2 = m_bayesianJointPdfPtr->lastComputedLogLikelihood();
414  }
415  }
416  } // if (bufferChar[0] == '1')
417 
419  // Prepare to exit routine or to stay in it
421  if (m_env.subRank() == 0) {
422  stayInRoutine = false; // Always for processor 0
423  }
424  else {
425  if (internalValues != NULL) delete internalValues;
426  if (internalDirection != NULL) delete internalDirection;
427  if (internalGrad != NULL) delete internalGrad;
428  if (internalHessian != NULL) delete internalHessian;
429  if (internalEffect != NULL) delete internalEffect;
430 
431  stayInRoutine = (vecValues == NULL) && (bufferChar[0] == '1');
432  //if (!stayInRoutine) std::cout << "Fullrank " << m_env.fullRank() << " is leaving scalarFunctionSync()" << std::endl;
433  }
434  } while (stayInRoutine);
435  }
436  else {
437  queso_require_msg(vecValues, "vecValues should not be NULL");
438 
439  m_env.subComm().Barrier();
440  result = m_scalarFunction.lnValue(*vecValues);
441  if (extraOutput1) {
442  if (m_bayesianJointPdfPtr) {
443  *extraOutput1 = m_bayesianJointPdfPtr->lastComputedLogPrior();
444  }
445  }
446  if (extraOutput2) {
447  if (m_bayesianJointPdfPtr) {
448  *extraOutput2 = m_bayesianJointPdfPtr->lastComputedLogLikelihood();
449  }
450  }
451  }
452 
453  return result;
454 }
455 
456 } // End namespace QUESO
457 
458 template class QUESO::ScalarFunctionSynchronizer<QUESO::GslVector, QUESO::GslMatrix>;
QUESO::ScalarFunctionSynchronizer::ScalarFunctionSynchronizer
ScalarFunctionSynchronizer(const BaseScalarFunction< V, M > &inputFunction, const V &auxVec)
Default constructor.
Definition: ScalarFunctionSynchronizer.C:37
QUESO::ScalarFunctionSynchronizer::domainSet
const VectorSet< V, M > & domainSet() const
Access to the domain set of the scalar function which will be synchronized.
Definition: ScalarFunctionSynchronizer.C:55
QUESO::VectorSet
A templated class for handling sets.
Definition: VectorSet.h:52
QUESO::BayesianJointPdf
A class for handling Bayesian joint PDFs.
Definition: ScalarFunctionSynchronizer.h:36
QUESO::BaseScalarFunction
A templated (base) class for handling scalar functions.
Definition: ScalarFunction.h:52
QUESO::ScalarFunctionSynchronizer::callFunction
double callFunction(const V *vecValues, const V *vecDirection, V *gradVector, M *hessianMatrix, V *hessianEffect, double *extraOutput1, double *extraOutput2) const
Calls the scalar function which will be synchronized.
Definition: ScalarFunctionSynchronizer.C:62
Generated on Sat Apr 22 2017 14:04:35 for queso-0.57.0 by  doxygen 1.8.5