queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ScalarSequence.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <algorithm>
26 #include <queso/ScalarSequence.h>
27 
28 namespace QUESO {
29 
30 // Default constructor -----------------------------
31 template <class T>
32 ScalarSequence<T>::ScalarSequence(
33  const BaseEnvironment& env,
34  unsigned int subSequenceSize,
35  const std::string& name)
36  :
37  m_env (env),
38  m_name (name),
39  m_seq (subSequenceSize,0.),
40  m_subMinPlain (NULL),
41  m_unifiedMinPlain (NULL),
42  m_subMaxPlain (NULL),
43  m_unifiedMaxPlain (NULL),
44  m_subMeanPlain (NULL),
45  m_unifiedMeanPlain (NULL),
46  m_subMedianPlain (NULL),
47  m_unifiedMedianPlain (NULL),
48  m_subSampleVariancePlain (NULL),
49  m_unifiedSampleVariancePlain(NULL)
50 {
51 }
52 // Destructor ---------------------------------------
53 template <class T>
54 ScalarSequence<T>::~ScalarSequence()
55 {
56  deleteStoredScalars();
57 }
58 // Set methods --------------------------------------
59 template <class T>
60 ScalarSequence<T>&
61 ScalarSequence<T>::operator= (const ScalarSequence<T>& rhs)
62 {
63  this->copy(rhs);
64  return *this;
65 }
66 // Access methods -----------------------------------
67 template <class T>
68 const T&
69 ScalarSequence<T>::operator[](unsigned int posId) const
70 {
71  if (posId >= this->subSequenceSize()) {
72  std::cerr << "In ScalarSequence<T>::operator[]() const"
73  << ": posId = " << posId
74  << ", this->subSequenceSize() = " << this->subSequenceSize()
75  << std::endl;
76  }
77  queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
78 
79  return m_seq[posId];
80 }
81 //---------------------------------------------------
82 template <class T>
83 T&
84 ScalarSequence<T>::operator[](unsigned int posId)
85 {
86  if (posId >= this->subSequenceSize()) {
87  std::cerr << "In ScalarSequence<T>::operator[]()"
88  << ": posId = " << posId
89  << ", this->subSequenceSize() = " << this->subSequenceSize()
90  << std::endl;
91  }
92  queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
93 
94  deleteStoredScalars();
95 
96  return m_seq[posId];
97 }
98 // Sequence methods ---------------------------------
99 template <class T>
100 const BaseEnvironment&
101 ScalarSequence<T>::env() const
102 {
103  return m_env;
104 }
105 // --------------------------------------------------
106 template <class T>
107 const std::string&
108 ScalarSequence<T>::name() const
109 {
110  return m_name;
111 }
112 // --------------------------------------------------
113 template <class T>
114 void
115 ScalarSequence<T>::setName(const std::string& newName)
116 {
117  m_name = newName;
118  return;
119 }
120 // --------------------------------------------------
121 template <class T>
122 void
123 ScalarSequence<T>::clear()
124 {
125  unsigned int numPos = this->subSequenceSize();
126  if (numPos) {
127  this->resetValues(0,numPos);
128  this->resizeSequence(0);
129  }
130 
131  return;
132 }
133 // --------------------------------------------------
134 template <class T>
135 unsigned int
136 ScalarSequence<T>::subSequenceSize() const
137 {
138  return m_seq.size();
139 }
140 // --------------------------------------------------
141 template <class T>
142 unsigned int
143 ScalarSequence<T>::unifiedSequenceSize(bool useOnlyInter0Comm) const
144 {
145  if (m_env.numSubEnvironments() == 1) {
146  return this->subSequenceSize();
147  }
148 
149  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
150 
151  unsigned int unifiedNumSamples = 0;
152  if (useOnlyInter0Comm) {
153  if (m_env.inter0Rank() >= 0) {
154  unsigned int subNumSamples = this->subSequenceSize();
155  m_env.inter0Comm().template Allreduce<unsigned int>(&subNumSamples, &unifiedNumSamples, (int) 1, RawValue_MPI_SUM,
156  "ScalarSequence<T>::unifiedSequenceSize()",
157  "failed MPI.Allreduce() for unifiedSequenceSize()");
158  }
159  else {
160  // Node not in the 'inter0' communicator
161  unifiedNumSamples = this->subSequenceSize();
162  }
163  }
164  else {
165  queso_error_msg("parallel vectors not supported yet");
166  }
167 
168  return unifiedNumSamples;
169 }
170 // --------------------------------------------------
171 template <class T>
172 void
173 ScalarSequence<T>::resizeSequence(unsigned int newSequenceSize)
174 {
175  if (newSequenceSize != this->subSequenceSize()) {
176  m_seq.resize(newSequenceSize,0.);
177  std::vector<T>(m_seq).swap(m_seq);
178  deleteStoredScalars();
179  }
180 
181  return;
182 }
183 // --------------------------------------------------
184 template <class T>
185 void
186 ScalarSequence<T>::resetValues(unsigned int initialPos, unsigned int numPos)
187 {
188  if (this->subSequenceSize() == 0) return;
189 
190  bool bRC = ((initialPos < this->subSequenceSize()) &&
191  (0 < numPos ) &&
192  ((initialPos+numPos) <= this->subSequenceSize()));
193  queso_require_msg(bRC, "invalid input data");
194 
195  for (unsigned int j = 0; j < numPos; ++j) {
196  m_seq[initialPos+j] = 0.;
197  }
198 
199  deleteStoredScalars();
200 
201  return;
202 }
203 // --------------------------------------------------
204 template <class T>
205 void
206 ScalarSequence<T>::erasePositions(unsigned int initialPos, unsigned int numPos)
207 {
208  if (this->subSequenceSize() == 0) return;
209 
210  bool bRC = ((initialPos < this->subSequenceSize()) &&
211  (0 < numPos ) &&
212  ((initialPos+numPos) <= this->subSequenceSize()));
213  queso_require_msg(bRC, "invalid input data");
214 
215  seqScalarPositionIteratorTypedef posIteratorBegin = m_seq.begin();
216  if (initialPos < this->subSequenceSize()) std::advance(posIteratorBegin,initialPos);
217  else posIteratorBegin = m_seq.end();
218 
219  unsigned int posEnd = initialPos + numPos;
220  seqScalarPositionIteratorTypedef posIteratorEnd = m_seq.begin();
221  if (posEnd < this->subSequenceSize()) std::advance(posIteratorEnd,posEnd);
222  else posIteratorEnd = m_seq.end();
223 
224  unsigned int oldSequenceSize = this->subSequenceSize();
225  m_seq.erase(posIteratorBegin,posIteratorEnd);
226  queso_require_equal_to_msg((oldSequenceSize - numPos), this->subSequenceSize(), "(oldSequenceSize - numPos) != this->subSequenceSize()");
227 
228  deleteStoredScalars();
229 
230  return;
231 }
232 // --------------------------------------------------
233 template <class T>
234 void
235 ScalarSequence<T>::getUnifiedContentsAtProc0Only(
236  bool useOnlyInter0Comm,
237  std::vector<T>& outputVec) const
238 {
239  // The logic (numSubEnvs == 1) does *not* apply here because 'outputVec' needs to be filled
240  //if (m_env.numSubEnvironments() == 1) {
241  // // No need to do anything
242  // return;
243  //}
244 
245  if (useOnlyInter0Comm) {
246  if (m_env.inter0Rank() >= 0) {
247  int auxSubSize = (int) this->subSequenceSize();
248  unsigned int auxUnifiedSize = this->unifiedSequenceSize(useOnlyInter0Comm);
249  outputVec.resize(auxUnifiedSize,0.);
250 
251  //******************************************************************
252  // Use MPI_Gatherv for the case different nodes have different amount of data // KAUST4
253  //******************************************************************
254  std::vector<int> recvcnts(m_env.inter0Comm().NumProc(),0); // '0' is NOT the correct value for recvcnts[0]
255  m_env.inter0Comm().template Gather<int>(&auxSubSize, 1, &recvcnts[0], (int) 1, 0,
256  "ScalarSequence<T>::getUnifiedContentsAtProc0Only()",
257  "failed MPI.Gather()");
258  if (m_env.inter0Rank() == 0) {
259  //recvcnts[0] = (int) this->subSequenceSize(); // FIX ME: really necessary????
260  queso_require_equal_to_msg(recvcnts[0], (int) this->subSequenceSize(), "failed MPI.Gather() result at proc 0");
261  }
262 
263  std::vector<int> displs(m_env.inter0Comm().NumProc(),0);
264  for (unsigned int r = 1; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // Yes, from '1' on
265  displs[r] = displs[r-1] + recvcnts[r-1];
266  }
267 
268 #if 0 // for debug only
269  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
270  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
271  *m_env.subDisplayFile() << " auxSubSize = " << auxSubSize
272  << ", recvcnts[" << r << "] = " << recvcnts[r]
273  << ", displs[" << r << "] = " << displs[r]
274  << ", m_seq.size() = " << m_seq.size()
275  << ", outputVec.size() = " << outputVec.size()
276  << std::endl;
277  }
278  for (unsigned int i = 0; i < m_seq.size(); ++i) {
279  *m_env.subDisplayFile() << " (before gatherv) m_seq[" << i << "]= " << m_seq[i]
280  << std::endl;
281  }
282  }
283 #endif
284 
285  m_env.inter0Comm().template Gatherv<double>(&m_seq[0], auxSubSize,
286  &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], 0,
287  "ScalarSequence<T>::getUnifiedContentsAtProc0Only()",
288  "failed MPI.Gatherv()");
289 
290 #if 0 // for debug only
291  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
292  for (unsigned int i = 0; i < m_seq.size(); ++i) {
293  *m_env.subDisplayFile() << " (after gatherv) m_seq[" << i << "]= " << m_seq[i]
294  << std::endl;
295  }
296  for (unsigned int i = 0; i < outputVec.size(); ++i) {
297  *m_env.subDisplayFile() << " (after gatherv) outputVec[" << i << "]= " << outputVec[i]
298  << std::endl;
299  }
300  }
301 #endif
302 
303 #if 0 // for debug only
304  if (m_env.inter0Rank() == 0) {
305  for (unsigned int i = 0; i < auxSubSize; ++i) {
306  outputVec[i] = m_seq[i];
307  }
308  m_env.inter0Comm().Gatherv(RawValue_MPI_IN_PLACE, auxSubSize, RawValue_MPI_DOUBLE, (void *) &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, 0,
309  "ScalarSequence<T>::getUnifiedContentsAtProc0Only(1)",
310  "failed MPI.Gatherv()");
311  }
312  else {
313  m_env.inter0Comm().Gatherv((void *) &m_seq[0], auxSubSize, RawValue_MPI_DOUBLE, (void *) &outputVec[0], (int *) &recvcnts[0], (int *) &displs[0], RawValue_MPI_DOUBLE, 0,
314  "ScalarSequence<T>::getUnifiedContentsAtProc0Only(2)",
315  "failed MPI.Gatherv()");
316  }
317  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
318  for (unsigned int i = 0; i < m_seq.size(); ++i) {
319  *m_env.subDisplayFile() << " (after 2nd gatherv) m_seq[" << i << "]= " << m_seq[i]
320  << std::endl;
321  }
322  for (unsigned int i = 0; i < outputVec.size(); ++i) {
323  *m_env.subDisplayFile() << " (after 2nd gatherv) outputVec[" << i << "]= " << outputVec[i]
324  << std::endl;
325  }
326  }
327 #endif
328  }
329  else {
330  // Node not in the 'inter0' communicator
331  // No need to do anything
332  }
333  }
334  else {
335  queso_error_msg("parallel vectors not supported yet");
336  }
337 
338  return;
339 }
340 // --------------------------------------------------
341 template <class T>
342 const T&
343 ScalarSequence<T>::subMinPlain() const
344 {
345  if (m_subMinPlain == NULL) {
346  m_subMinPlain = new T(0.);
347  if (m_subMaxPlain == NULL) m_subMaxPlain = new T(0.);
348  subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
349  }
350 
351  return *m_subMinPlain;
352 }
353 // --------------------------------------------------
354 template <class T>
355 const T&
356 ScalarSequence<T>::unifiedMinPlain(bool useOnlyInter0Comm) const
357 {
358  if (m_unifiedMinPlain == NULL) {
359  m_unifiedMinPlain = new T(0.);
360  if (m_unifiedMaxPlain == NULL) m_unifiedMaxPlain = new T(0.);
361  unifiedMinMaxExtra(useOnlyInter0Comm,0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
362  }
363 
364  return *m_unifiedMinPlain;
365 }
366 // --------------------------------------------------
367 template <class T>
368 const T&
369 ScalarSequence<T>::subMaxPlain() const
370 {
371  if (m_subMaxPlain == NULL) {
372  if (m_subMinPlain == NULL) m_subMinPlain = new T(0.);
373  m_subMaxPlain = new T(0.);
374  subMinMaxExtra(0,this->subSequenceSize(),*m_subMinPlain,*m_subMaxPlain);
375  }
376 
377  return *m_subMaxPlain;
378 }
379 // --------------------------------------------------
380 template <class T>
381 const T&
382 ScalarSequence<T>::unifiedMaxPlain(bool useOnlyInter0Comm) const
383 {
384  if (m_unifiedMaxPlain == NULL) {
385  if (m_unifiedMinPlain == NULL) m_unifiedMinPlain = new T(0.);
386  m_unifiedMaxPlain = new T(0.);
387  unifiedMinMaxExtra(useOnlyInter0Comm,0,this->subSequenceSize(),*m_unifiedMinPlain,*m_unifiedMaxPlain);
388  }
389 
390  return *m_unifiedMaxPlain;
391 }
392 // --------------------------------------------------
393 template <class T>
394 const T&
395 ScalarSequence<T>::subMeanPlain() const
396 {
397  if (m_subMeanPlain == NULL) {
398  m_subMeanPlain = new T(0.);
399  *m_subMeanPlain = subMeanExtra(0,subSequenceSize());
400  }
401 
402  return *m_subMeanPlain;
403 }
404 // --------------------------------------------------
405 template <class T>
406 const T&
407 ScalarSequence<T>::unifiedMeanPlain(bool useOnlyInter0Comm) const
408 {
409  if (m_unifiedMeanPlain == NULL) {
410  m_unifiedMeanPlain = new T(0.);
411  *m_unifiedMeanPlain = unifiedMeanExtra(useOnlyInter0Comm,0,subSequenceSize());
412  }
413 
414  return *m_unifiedMeanPlain;
415 }
416 // --------------------------------------------------
417 template <class T>
418 const T&
419 ScalarSequence<T>::subMedianPlain() const
420 {
421  if (m_subMedianPlain == NULL) {
422  m_subMedianPlain = new T(0.);
423  *m_subMedianPlain = subMedianExtra(0,subSequenceSize());
424  }
425 
426  return *m_subMedianPlain;
427 }
428 // --------------------------------------------------
429 template <class T>
430 const T&
431 ScalarSequence<T>::unifiedMedianPlain(bool useOnlyInter0Comm) const
432 {
433  if (m_unifiedMedianPlain == NULL) {
434  m_unifiedMedianPlain = new T(0.);
435  *m_unifiedMedianPlain = unifiedMedianExtra(useOnlyInter0Comm,0,subSequenceSize());
436  }
437 
438  return *m_unifiedMedianPlain;
439 }
440 // --------------------------------------------------
441 template <class T>
442 const T&
443 ScalarSequence<T>::subSampleVariancePlain() const
444 {
445  if (m_subSampleVariancePlain == NULL) {
446  m_subSampleVariancePlain = new T(0.);
447  *m_subSampleVariancePlain = subSampleVarianceExtra(0,subSequenceSize(),subMeanPlain());
448  }
449 
450  return *m_subSampleVariancePlain;
451 }
452 // --------------------------------------------------
453 template <class T>
454 const T&
455 ScalarSequence<T>::unifiedSampleVariancePlain(bool useOnlyInter0Comm) const
456 {
457  if (m_unifiedSampleVariancePlain == NULL) {
458  m_unifiedSampleVariancePlain = new T(0.);
459  *m_unifiedSampleVariancePlain = unifiedSampleVarianceExtra(useOnlyInter0Comm,0,subSequenceSize(),unifiedMeanPlain(useOnlyInter0Comm));
460  }
461 
462  return *m_unifiedSampleVariancePlain;
463 }
464 // --------------------------------------------------
465 template <class T>
466 void
467 ScalarSequence<T>::deleteStoredScalars()
468 {
469  if (m_subMinPlain) {
470  delete m_subMinPlain;
471  m_subMinPlain = NULL;
472  }
473  if (m_unifiedMinPlain) {
474  delete m_unifiedMinPlain;
475  m_unifiedMinPlain = NULL;
476  }
477  if (m_subMaxPlain) {
478  delete m_subMaxPlain;
479  m_subMaxPlain = NULL;
480  }
481  if (m_unifiedMaxPlain) {
482  delete m_unifiedMaxPlain;
483  m_unifiedMaxPlain = NULL;
484  }
485  if (m_subMeanPlain) {
486  delete m_subMeanPlain;
487  m_subMeanPlain = NULL;
488  }
489  if (m_unifiedMeanPlain) {
490  delete m_unifiedMeanPlain;
491  m_unifiedMeanPlain = NULL;
492  }
493  if (m_subMedianPlain) {
494  delete m_subMedianPlain;
495  m_subMedianPlain = NULL;
496  }
497  if (m_unifiedMedianPlain) {
498  delete m_unifiedMedianPlain;
499  m_unifiedMedianPlain = NULL;
500  }
501  if (m_subSampleVariancePlain) {
502  delete m_subSampleVariancePlain;
503  m_subSampleVariancePlain = NULL;
504  }
505  if (m_unifiedSampleVariancePlain) {
506  delete m_unifiedSampleVariancePlain;
507  m_unifiedSampleVariancePlain = NULL;
508  }
509 
510  return;
511 }
512 // --------------------------------------------------
513 template <class T>
514 void
515 ScalarSequence<T>::setGaussian(const T& meanValue, const T& stdDev)
516 {
517  unsigned int maxJ = this->subSequenceSize();
518  for (unsigned int j = 0; j < maxJ; ++j) {
519  m_seq[j] = meanValue + m_env.rngObject()->gaussianSample(stdDev);
520  }
521 
522  deleteStoredScalars();
523 }
524 // --------------------------------------------------
525 template <class T>
526 void
527 ScalarSequence<T>::setUniform(const T& a, const T& b)
528 {
529  unsigned int maxJ = this->subSequenceSize();
530  for (unsigned int j = 0; j < maxJ; ++j) {
531  m_seq[j] = a + (b-a)*m_env.rngObject()->uniformSample();
532  }
533 
534  deleteStoredScalars();
535 }
536 // --------------------------------------------------
537 template <class T>
538 void
539 ScalarSequence<T>::subUniformlySampledCdf(
540  unsigned int numEvaluationPoints,
541  T& minDomainValue,
542  T& maxDomainValue,
543  std::vector<T>& cdfValues) const
544 {
545  T tmpMinValue;
546  T tmpMaxValue;
547  std::vector<T> centers(numEvaluationPoints,0.);
548  std::vector<unsigned int> bins (numEvaluationPoints,0);
549 
550  subMinMaxExtra(0, // initialPos
551  this->subSequenceSize(),
552  tmpMinValue,
553  tmpMaxValue);
554  subHistogram(0, // initialPos,
555  tmpMinValue,
556  tmpMaxValue,
557  centers,
558  bins);
559 
560  minDomainValue = centers[0];
561  maxDomainValue = centers[centers.size()-1];
562 
563  unsigned int sumOfBins = 0;
564  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
565  sumOfBins += bins[i];
566  }
567 
568  cdfValues.clear();
569  cdfValues.resize(numEvaluationPoints);
570  unsigned int partialSum = 0;
571  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
572  partialSum += bins[i];
573  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
574  }
575 
576  return;
577 }
578 // --------------------------------------------------
579 template <class T>
580 void
581 ScalarSequence<T>::unifiedUniformlySampledCdf(
582  bool useOnlyInter0Comm,
583  unsigned int numEvaluationPoints,
584  T& unifiedMinDomainValue,
585  T& unifiedMaxDomainValue,
586  std::vector<T>& unifiedCdfValues) const
587 {
588  if (m_env.numSubEnvironments() == 1) {
589  return this->subUniformlySampledCdf(numEvaluationPoints,
590  unifiedMinDomainValue,
591  unifiedMaxDomainValue,
592  unifiedCdfValues);
593  }
594 
595  // KAUST2
596  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
597 
598  if (useOnlyInter0Comm) {
599  if (m_env.inter0Rank() >= 0) {
600  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
601  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedUniformlySampledCdf()"
602  << std::endl;
603  }
604 
605  T unifiedTmpMinValue;
606  T unifiedTmpMaxValue;
607  std::vector<T> unifiedCenters(numEvaluationPoints,0.);
608  std::vector<unsigned int> unifiedBins (numEvaluationPoints,0);
609 
610  this->unifiedMinMaxExtra(useOnlyInter0Comm,
611  0, // initialPos
612  this->subSequenceSize(),
613  unifiedTmpMinValue,
614  unifiedTmpMaxValue);
615  this->unifiedHistogram(useOnlyInter0Comm,
616  0, // initialPos
617  unifiedTmpMinValue,
618  unifiedTmpMaxValue,
619  unifiedCenters,
620  unifiedBins);
621 
622  unifiedMinDomainValue = unifiedCenters[0];
623  unifiedMaxDomainValue = unifiedCenters[unifiedCenters.size()-1];
624 
625  unsigned int unifiedTotalSumOfBins = 0;
626  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
627  unifiedTotalSumOfBins += unifiedBins[i];
628  }
629 
630  std::vector<unsigned int> unifiedPartialSumsOfBins(numEvaluationPoints,0);
631  unifiedPartialSumsOfBins[0] = unifiedBins[0];
632  for (unsigned int i = 1; i < numEvaluationPoints; ++i) { // Yes, from '1'
633  unifiedPartialSumsOfBins[i] = unifiedPartialSumsOfBins[i-1] + unifiedBins[i];
634  }
635 
636  unifiedCdfValues.clear();
637  unifiedCdfValues.resize(numEvaluationPoints);
638  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
639  unifiedCdfValues[i] = ((T) unifiedPartialSumsOfBins[i])/((T) unifiedTotalSumOfBins);
640  }
641 
642  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
643  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
644  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedUniformlySampledCdf()"
645  << ": i = " << i
646  << ", unifiedTmpMinValue = " << unifiedTmpMinValue
647  << ", unifiedTmpMaxValue = " << unifiedTmpMaxValue
648  << ", unifiedBins = " << unifiedBins[i]
649  << ", unifiedCdfValue = " << unifiedCdfValues[i]
650  << ", unifiedPartialSumsOfBins = " << unifiedPartialSumsOfBins[i]
651  << ", unifiedTotalSumOfBins = " << unifiedTotalSumOfBins
652  << std::endl;
653  }
654  }
655 
656  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
657  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedUniformlySampledCdf()"
658  << std::endl;
659  }
660  }
661  else {
662  // Node not in the 'inter0' communicator
663  this->subUniformlySampledCdf(numEvaluationPoints,
664  unifiedMinDomainValue,
665  unifiedMaxDomainValue,
666  unifiedCdfValues);
667  }
668  }
669  else {
670  queso_error_msg("parallel vectors not supported yet");
671  }
672 
673  //m_env.fullComm().Barrier();
674 
675  return;
676 }
677 // --------------------------------------------------
678 template <class T>
679 void
680 ScalarSequence<T>::subBasicCdf(
681  unsigned int numEvaluationPoints,
682  UniformOneDGrid<T>*& gridValues,
683  std::vector<T>& cdfValues) const
684 {
685  T tmpMinValue;
686  T tmpMaxValue;
687  std::vector<unsigned int> bins(numEvaluationPoints,0);
688 
689  subMinMaxExtra(0, // initialPos
690  this->subSequenceSize(),
691  tmpMinValue,
692  tmpMaxValue);
693  subBasicHistogram(0, // initialPos,
694  tmpMinValue,
695  tmpMaxValue,
696  gridValues,
697  bins);
698 
699  unsigned int sumOfBins = 0;
700  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
701  sumOfBins += bins[i];
702  }
703 
704  cdfValues.clear();
705  cdfValues.resize(numEvaluationPoints);
706  unsigned int partialSum = 0;
707  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
708  partialSum += bins[i];
709  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
710  }
711 
712  return;
713 }
714 // --------------------------------------------------
715 template <class T>
716 void
717 ScalarSequence<T>::subWeightCdf(
718  unsigned int numEvaluationPoints,
719  std::vector<T>& gridValues,
720  std::vector<T>& cdfValues) const
721 {
722  T tmpMinValue;
723  T tmpMaxValue;
724  std::vector<unsigned int> bins(numEvaluationPoints,0);
725  gridValues.resize (numEvaluationPoints,0.);
726  cdfValues.resize (numEvaluationPoints,0.);
727 
728  subMinMaxExtra(0, // initialPos
729  this->subSequenceSize(),
730  tmpMinValue,
731  tmpMaxValue);
732 
733  if (tmpMinValue == tmpMaxValue) {
734  if (tmpMinValue < -1.e-12) {
735  tmpMinValue += tmpMinValue*(1.e-8);
736  }
737  else if (tmpMinValue > 1.e-12) {
738  tmpMinValue -= tmpMinValue*(1.e-8);
739  }
740  else {
741  tmpMinValue = 1.e-12;
742  }
743  }
744 
745  subWeightHistogram(0, // initialPos,
746  tmpMinValue,
747  tmpMaxValue,
748  gridValues,
749  bins);
750 
751  unsigned int sumOfBins = 0;
752  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
753  sumOfBins += bins[i];
754  }
755 
756  cdfValues.clear();
757  cdfValues.resize(numEvaluationPoints);
758  unsigned int partialSum = 0;
759  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
760  partialSum += bins[i];
761  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
762  }
763 
764  return;
765 }
766 // --------------------------------------------------
767 template <class T>
768 void
769 ScalarSequence<T>::subWeightCdf(
770  unsigned int numEvaluationPoints,
771  UniformOneDGrid<T>*& gridValues,
772  std::vector<T>& cdfValues) const
773 {
774  T tmpMinValue;
775  T tmpMaxValue;
776  std::vector<unsigned int> bins(numEvaluationPoints,0);
777 
778  subMinMaxExtra(0, // initialPos
779  this->subSequenceSize(),
780  tmpMinValue,
781  tmpMaxValue);
782  subWeightHistogram(0, // initialPos,
783  tmpMinValue,
784  tmpMaxValue,
785  gridValues,
786  bins);
787 
788  unsigned int sumOfBins = 0;
789  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
790  sumOfBins += bins[i];
791  }
792 
793  cdfValues.clear();
794  cdfValues.resize(numEvaluationPoints);
795  unsigned int partialSum = 0;
796  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
797  partialSum += bins[i];
798  cdfValues[i] = ((T) partialSum)/((T) sumOfBins);
799  }
800 
801  return;
802 }
803 // --------------------------------------------------
804 template <class T>
805 T
806 ScalarSequence<T>::subMeanExtra(
807  unsigned int initialPos,
808  unsigned int numPos) const
809 {
810  if (this->subSequenceSize() == 0) return 0.;
811 
812  bool bRC = ((initialPos < this->subSequenceSize()) &&
813  (0 < numPos ) &&
814  ((initialPos+numPos) <= this->subSequenceSize()));
815  if (bRC == false) {
816  std::cerr << "In ScalarSequence<T>::subMeanExtra()"
817  << ": ERROR at fullRank " << m_env.fullRank()
818  << ", initialPos = " << initialPos
819  << ", numPos = " << numPos
820  << ", this->subSequenceSize() = " << this->subSequenceSize()
821  << std::endl;
822  if (m_env.subDisplayFile()) {
823  *m_env.subDisplayFile() << "In ScalarSequence<T>::subMeanExtra()"
824  << ": ERROR at fullRank " << m_env.fullRank()
825  << ", initialPos = " << initialPos
826  << ", numPos = " << numPos
827  << ", this->subSequenceSize() = " << this->subSequenceSize()
828  << std::endl;
829  }
830  }
831  queso_require_msg(bRC, "invalid input data");
832 
833  unsigned int finalPosPlus1 = initialPos + numPos;
834  T tmpSum = 0.;
835  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
836  tmpSum += m_seq[j];
837  }
838 
839  return tmpSum/(T) numPos;
840 }
841 // --------------------------------------------------
842 template <class T>
843 T
844 ScalarSequence<T>::unifiedMeanExtra(
845  bool useOnlyInter0Comm,
846  unsigned int initialPos,
847  unsigned int numPos) const
848 {
849  if (m_env.numSubEnvironments() == 1) {
850  return this->subMeanExtra(initialPos,
851  numPos);
852  }
853 
854  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
855 
856  T unifiedMeanValue = 0.;
857  if (useOnlyInter0Comm) {
858  if (m_env.inter0Rank() >= 0) {
859  bool bRC = ((initialPos < this->subSequenceSize()) &&
860  (0 < numPos ) &&
861  ((initialPos+numPos) <= this->subSequenceSize()));
862  queso_require_msg(bRC, "invalid input data");
863 
864  unsigned int finalPosPlus1 = initialPos + numPos;
865  T localSum = 0.;
866  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
867  localSum += m_seq[j];
868  }
869 
870  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
871  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
872  << ": initialPos = " << initialPos
873  << ", numPos = " << numPos
874  << ", before MPI.Allreduce"
875  << std::endl;
876  }
877  //std::cout << m_env.inter0Comm().MyPID()
878  // << std::endl;
879  //sleep(1);
880  unsigned int unifiedNumPos = 0;
881  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
882  "ScalarSequence<T>::unifiedMeanExtra()",
883  "failed MPI.Allreduce() for numPos");
884 
885  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
886  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
887  << ": numPos = " << numPos
888  << ", unifiedNumPos = " << unifiedNumPos
889  << std::endl;
890  }
891 
892  m_env.inter0Comm().template Allreduce<double>(&localSum, &unifiedMeanValue, (int) 1, RawValue_MPI_SUM,
893  "ScalarSequence<T>::unifiedMeanExtra()",
894  "failed MPI.Allreduce() for sum");
895 
896  unifiedMeanValue /= ((T) unifiedNumPos);
897 
898  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
899  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMeanExtra()"
900  << ": localSum = " << localSum
901  << ", unifiedMeanValue = " << unifiedMeanValue
902  << std::endl;
903  }
904  }
905  else {
906  // Node not in the 'inter0' communicator
907  this->subMeanExtra(initialPos,
908  numPos);
909  }
910  }
911  else {
912  queso_error_msg("parallel vectors not supported yet");
913  }
914 
915  //m_env.fullComm().Barrier();
916 
917  return unifiedMeanValue;
918 }
919 // --------------------------------------------------
920 template <class T>
921 T
922 ScalarSequence<T>::subMedianExtra(
923  unsigned int initialPos,
924  unsigned int numPos) const
925 {
926  if (this->subSequenceSize() == 0) return 0.;
927 
928  bool bRC = ((initialPos < this->subSequenceSize()) &&
929  (0 < numPos ) &&
930  ((initialPos+numPos) <= this->subSequenceSize()));
931  if (bRC == false) {
932  std::cerr << "In ScalarSequence<T>::subMedianExtra()"
933  << ": ERROR at fullRank " << m_env.fullRank()
934  << ", initialPos = " << initialPos
935  << ", numPos = " << numPos
936  << ", this->subSequenceSize() = " << this->subSequenceSize()
937  << std::endl;
938  if (m_env.subDisplayFile()) {
939  *m_env.subDisplayFile() << "In ScalarSequence<T>::subMedianExtra()"
940  << ": ERROR at fullRank " << m_env.fullRank()
941  << ", initialPos = " << initialPos
942  << ", numPos = " << numPos
943  << ", this->subSequenceSize() = " << this->subSequenceSize()
944  << std::endl;
945  }
946  }
947  queso_require_msg(bRC, "invalid input data");
948 
949  ScalarSequence sortedSequence(m_env,0,"");
950  sortedSequence.resizeSequence(numPos);
951  this->extractScalarSeq(initialPos,
952  1,
953  numPos,
954  sortedSequence);
955  sortedSequence.subSort();
956 
957  unsigned int tmpPos = (unsigned int) (0.5 * (double) numPos);
958  T resultValue = sortedSequence[tmpPos];
959 
960  return resultValue;
961 }
962 // --------------------------------------------------
963 template <class T>
964 T
965 ScalarSequence<T>::unifiedMedianExtra( // rr0
966  bool useOnlyInter0Comm,
967  unsigned int initialPos,
968  unsigned int numPos) const
969 {
970  if (m_env.numSubEnvironments() == 1) {
971  return this->subMedianExtra(initialPos,
972  numPos);
973  }
974 
975  // As of 07/Jul/2012, this routine does *not* require sub sequences to have equal size. Good.
976 
977  T unifiedMedianValue = 0.;
978  if (useOnlyInter0Comm) {
979  if (m_env.inter0Rank() >= 0) {
980  bool bRC = ((initialPos < this->subSequenceSize()) &&
981  (0 < numPos ) &&
982  ((initialPos+numPos) <= this->subSequenceSize()));
983  queso_require_msg(bRC, "invalid input data");
984 
985  ScalarSequence unifiedSortedSequence(m_env,0,"");
986  this->unifiedSort(useOnlyInter0Comm,
987  initialPos,
988  unifiedSortedSequence);
989  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
990  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMedianExtra()"
991  << ", unifiedMedianValue = " << unifiedMedianValue
992  << std::endl;
993  }
994  }
995  else {
996  // Node not in the 'inter0' communicator
997  this->subMedianExtra(initialPos,
998  numPos);
999  }
1000  }
1001  else {
1002  queso_error_msg("parallel vectors not supported yet");
1003  }
1004 
1005  //m_env.fullComm().Barrier();
1006 
1007  return unifiedMedianValue;
1008 }
1009 // --------------------------------------------------
1010 template <class T>
1011 T
1012 ScalarSequence<T>::subSampleVarianceExtra(
1013  unsigned int initialPos,
1014  unsigned int numPos,
1015  const T& meanValue) const
1016 {
1017  if (this->subSequenceSize() == 0) return 0.;
1018 
1019  bool bRC = ((initialPos < this->subSequenceSize()) &&
1020  (0 < numPos ) &&
1021  ((initialPos+numPos) <= this->subSequenceSize()));
1022  queso_require_msg(bRC, "invalid input data");
1023 
1024  unsigned int finalPosPlus1 = initialPos + numPos;
1025  T diff;
1026  T samValue = 0.;
1027  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1028  diff = m_seq[j] - meanValue;
1029  samValue += diff*diff;
1030  }
1031 
1032  samValue /= (((T) numPos) - 1.);
1033 
1034  return samValue;
1035 }
1036 // --------------------------------------------------
1037 template <class T>
1038 T
1039 ScalarSequence<T>::unifiedSampleVarianceExtra(
1040  bool useOnlyInter0Comm,
1041  unsigned int initialPos,
1042  unsigned int numPos,
1043  const T& unifiedMeanValue) const
1044 {
1045  if (m_env.numSubEnvironments() == 1) {
1046  return this->subSampleVarianceExtra(initialPos,
1047  numPos,
1048  unifiedMeanValue);
1049  }
1050 
1051  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1052 
1053  T unifiedSamValue = 0.;
1054  if (useOnlyInter0Comm) {
1055  if (m_env.inter0Rank() >= 0) {
1056  bool bRC = ((initialPos < this->subSequenceSize()) &&
1057  (0 < numPos ) &&
1058  ((initialPos+numPos) <= this->subSequenceSize()));
1059  queso_require_msg(bRC, "invalid input data");
1060 
1061  unsigned int finalPosPlus1 = initialPos + numPos;
1062  T diff;
1063  T localSamValue = 0.;
1064  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1065  diff = m_seq[j] - unifiedMeanValue;
1066  localSamValue += diff*diff;
1067  }
1068 
1069  unsigned int unifiedNumPos = 0;
1070  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
1071  "ScalarSequence<T>::unifiedSampleVarianceExtra()",
1072  "failed MPI.Allreduce() for numPos");
1073 
1074  m_env.inter0Comm().template Allreduce<double>(&localSamValue, &unifiedSamValue, (int) 1, RawValue_MPI_SUM,
1075  "ScalarSequence<T>::unifiedSampleVarianceExtra()",
1076  "failed MPI.Allreduce() for samValue");
1077 
1078  unifiedSamValue /= (((T) unifiedNumPos) - 1.);
1079  }
1080  else {
1081  // Node not in the 'inter0' communicator
1082  this->subSampleVarianceExtra(initialPos,
1083  numPos,
1084  unifiedMeanValue);
1085  }
1086  }
1087  else {
1088  queso_error_msg("parallel vectors not supported yet");
1089  }
1090 
1091  //m_env.fullComm().Barrier();
1092 
1093  return unifiedSamValue;
1094 }
1095 // --------------------------------------------------
1096 template <class T>
1097 T
1098 ScalarSequence<T>::subSampleStd(
1099  unsigned int initialPos,
1100  unsigned int numPos,
1101  const T& meanValue) const
1102 {
1103  if (this->subSequenceSize() == 0) return 0.;
1104 
1105  bool bRC = ((initialPos < this->subSequenceSize()) &&
1106  (0 < numPos ) &&
1107  ((initialPos+numPos) <= this->subSequenceSize()));
1108  queso_require_msg(bRC, "invalid input data");
1109 
1110  unsigned int finalPosPlus1 = initialPos + numPos;
1111  T diff;
1112  T stdValue = 0.;
1113  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1114  diff = m_seq[j] - meanValue;
1115  stdValue += diff*diff;
1116  }
1117 
1118  stdValue /= (((T) numPos) - 1.);
1119  stdValue = sqrt(stdValue);
1120 
1121  return stdValue;
1122 }
1123 // --------------------------------------------------
1124 template <class T>
1125 T
1126 ScalarSequence<T>::unifiedSampleStd(
1127  bool useOnlyInter0Comm,
1128  unsigned int initialPos,
1129  unsigned int numPos,
1130  const T& unifiedMeanValue) const
1131 {
1132  if (m_env.numSubEnvironments() == 1) {
1133  return this->subSampleStd(initialPos,
1134  numPos,
1135  unifiedMeanValue);
1136  }
1137 
1138  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1139 
1140  T unifiedStdValue = 0.;
1141  if (useOnlyInter0Comm) {
1142  if (m_env.inter0Rank() >= 0) {
1143  bool bRC = ((initialPos < this->subSequenceSize()) &&
1144  (0 < numPos ) &&
1145  ((initialPos+numPos) <= this->subSequenceSize()));
1146  queso_require_msg(bRC, "invalid input data");
1147 
1148  unsigned int finalPosPlus1 = initialPos + numPos;
1149  T diff;
1150  T localStdValue = 0.;
1151  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1152  diff = m_seq[j] - unifiedMeanValue;
1153  localStdValue += diff*diff;
1154  }
1155 
1156  unsigned int unifiedNumPos = 0;
1157  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
1158  "ScalarSequence<T>::unifiedSampleStd()",
1159  "failed MPI.Allreduce() for numPos");
1160 
1161  m_env.inter0Comm().template Allreduce<double>(&localStdValue, &unifiedStdValue, (int) 1, RawValue_MPI_SUM,
1162  "ScalarSequence<T>::unifiedSampleStd()",
1163  "failed MPI.Allreduce() for stdValue");
1164 
1165  unifiedStdValue /= (((T) unifiedNumPos) - 1.);
1166  unifiedStdValue = sqrt(unifiedStdValue);
1167  }
1168  else {
1169  // Node not in the 'inter0' communicator
1170  this->subSampleStd(initialPos,
1171  numPos,
1172  unifiedMeanValue);
1173  }
1174  }
1175  else {
1176  queso_error_msg("parallel vectors not supported yet");
1177  }
1178 
1179  //m_env.fullComm().Barrier();
1180 
1181  return unifiedStdValue;
1182 }
1183 // --------------------------------------------------
1184 template <class T>
1185 T
1186 ScalarSequence<T>::subPopulationVariance(
1187  unsigned int initialPos,
1188  unsigned int numPos,
1189  const T& meanValue) const
1190 {
1191  if (this->subSequenceSize() == 0) return 0.;
1192 
1193  bool bRC = ((initialPos < this->subSequenceSize()) &&
1194  (0 < numPos ) &&
1195  ((initialPos+numPos) <= this->subSequenceSize()));
1196  queso_require_msg(bRC, "invalid input data");
1197 
1198  unsigned int finalPosPlus1 = initialPos + numPos;
1199  T diff;
1200  T popValue = 0.;
1201  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1202  diff = m_seq[j] - meanValue;
1203  popValue += diff*diff;
1204  }
1205 
1206  popValue /= (T) numPos;
1207 
1208  return popValue;
1209 }
1210 // --------------------------------------------------
1211 template <class T>
1212 T
1213 ScalarSequence<T>::unifiedPopulationVariance(
1214  bool useOnlyInter0Comm,
1215  unsigned int initialPos,
1216  unsigned int numPos,
1217  const T& unifiedMeanValue) const
1218 {
1219  if (m_env.numSubEnvironments() == 1) {
1220  return this->subPopulationVariance(initialPos,
1221  numPos,
1222  unifiedMeanValue);
1223  }
1224 
1225  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1226 
1227  T unifiedPopValue = 0.;
1228  if (useOnlyInter0Comm) {
1229  if (m_env.inter0Rank() >= 0) {
1230  bool bRC = ((initialPos < this->subSequenceSize()) &&
1231  (0 < numPos ) &&
1232  ((initialPos+numPos) <= this->subSequenceSize()));
1233  queso_require_msg(bRC, "invalid input data");
1234 
1235  unsigned int finalPosPlus1 = initialPos + numPos;
1236  T diff;
1237  T localPopValue = 0.;
1238  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1239  diff = m_seq[j] - unifiedMeanValue;
1240  localPopValue += diff*diff;
1241  }
1242 
1243  unsigned int unifiedNumPos = 0;
1244  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
1245  "ScalarSequence<T>::unifiedPopulationVariance()",
1246  "failed MPI.Allreduce() for numPos");
1247 
1248  m_env.inter0Comm().template Allreduce<double>(&localPopValue, &unifiedPopValue, (int) 1, RawValue_MPI_SUM,
1249  "ScalarSequence<T>::unifiedPopulationVariance()",
1250  "failed MPI.Allreduce() for popValue");
1251 
1252  unifiedPopValue /= ((T) unifiedNumPos);
1253  }
1254  else {
1255  // Node not in the 'inter0' communicator
1256  this->subPopulationVariance(initialPos,
1257  numPos,
1258  unifiedMeanValue);
1259  }
1260  }
1261  else {
1262  queso_error_msg("parallel vectors not supported yet");
1263  }
1264 
1265  //m_env.fullComm().Barrier();
1266 
1267  return unifiedPopValue;
1268 }
1269 // --------------------------------------------------
1270 template <class T>
1271 T
1272 ScalarSequence<T>::autoCovariance(
1273  unsigned int initialPos,
1274  unsigned int numPos,
1275  const T& meanValue,
1276  unsigned int lag) const
1277 {
1278  bool bRC = ((initialPos < this->subSequenceSize()) &&
1279  (0 < numPos ) &&
1280  ((initialPos+numPos) <= this->subSequenceSize()) &&
1281  (lag < numPos )); // lag should not be too large
1282  queso_require_msg(bRC, "invalid input data");
1283 
1284  unsigned int loopSize = numPos - lag;
1285  unsigned int finalPosPlus1 = initialPos + loopSize;
1286  T diff1;
1287  T diff2;
1288  T covValue = 0.;
1289  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
1290  diff1 = m_seq[j ] - meanValue;
1291  diff2 = m_seq[j+lag] - meanValue;
1292  covValue += diff1*diff2;
1293  }
1294 
1295  covValue /= (T) loopSize;
1296 
1297  return covValue;
1298 }
1299 // --------------------------------------------------
1300 template <class T>
1301 T
1302 ScalarSequence<T>::autoCorrViaDef(
1303  unsigned int initialPos,
1304  unsigned int numPos,
1305  unsigned int lag) const
1306 {
1307  bool bRC = ((initialPos < this->subSequenceSize()) &&
1308  (0 < numPos ) &&
1309  ((initialPos+numPos) <= this->subSequenceSize()) &&
1310  (lag < numPos )); // lag should not be too large
1311  queso_require_msg(bRC, "invalid input data");
1312 
1313  T meanValue = this->subMeanExtra(initialPos,
1314  numPos);
1315 
1316  T covValueZero = this->autoCovariance(initialPos,
1317  numPos,
1318  meanValue,
1319  0); // lag
1320 
1321  T corrValue = this->autoCovariance(initialPos,
1322  numPos,
1323  meanValue,
1324  lag);
1325 
1326  return corrValue/covValueZero;
1327 }
1328 // --------------------------------------------------
1329 template <class T>
1330 void
1331 ScalarSequence<T>::autoCorrViaFft(
1332  unsigned int initialPos,
1333  unsigned int numPos,
1334  unsigned int maxLag,
1335  std::vector<T>& autoCorrs) const
1336 {
1337  double tmp = log((double) numPos)/log(2.);
1338  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1339  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1340  unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
1341 
1342  std::vector<double> rawDataVec(numPos,0.);
1343  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1344  Fft<T> fftObj(m_env);
1345 
1346  // Forward FFT
1347  this->extractRawData(initialPos,
1348  1, // spacing
1349  numPos,
1350  rawDataVec);
1351  T meanValue = this->subMeanExtra(initialPos,
1352  numPos);
1353  for (unsigned int j = 0; j < numPos; ++j) {
1354  rawDataVec[j] -= meanValue; // IMPORTANT
1355  }
1356 
1357  rawDataVec.resize(fftSize,0.);
1358 
1359  //if (m_env.subDisplayFile()) {
1360  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1361  // << ": about to call fftObj.forward()"
1362  // << " with rawDataVec.size() = " << rawDataVec.size()
1363  // << ", fftSize = " << fftSize
1364  // << ", resultData.size() = " << resultData.size()
1365  // << std::endl;
1366  //}
1367  fftObj.forward(rawDataVec,fftSize,resultData);
1368 
1369  // Inverse FFT
1370  for (unsigned int j = 0; j < fftSize; ++j) {
1371  rawDataVec[j] = std::norm(resultData[j]);
1372  }
1373  //if (m_env.subDisplayFile()) {
1374  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1375  // << ": about to call fftObj.inverse()"
1376  // << " with rawDataVec.size() = " << rawDataVec.size()
1377  // << ", fftSize = " << fftSize
1378  // << ", resultData.size() = " << resultData.size()
1379  // << std::endl;
1380  //}
1381  fftObj.inverse(rawDataVec,fftSize,resultData);
1382  //if (m_env.subDisplayFile()) {
1383  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1384  // << ": returned succesfully from fftObj.inverse()"
1385  // << std::endl;
1386  //}
1387 
1388  // Prepare return data
1389  autoCorrs.resize(maxLag+1,0.); // Yes, +1
1390  for (unsigned int j = 0; j < autoCorrs.size(); ++j) {
1391  double ratio = ((double) j)/((double) (numPos-1));
1392  autoCorrs[j] = ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1393  }
1394 
1395  return;
1396 }
1397 // --------------------------------------------------
1398 template <class T>
1399 void
1400 ScalarSequence<T>::autoCorrViaFft(
1401  unsigned int initialPos,
1402  unsigned int numPos,
1403  unsigned int numSum,
1404  T& autoCorrsSum) const
1405 {
1406  //if (m_env.subDisplayFile()) {
1407  // *m_env.subDisplayFile() << "Entering ScalarSequence<T>::autoCorrViaFft(), for sum"
1408  // << ": initialPos = " << initialPos
1409  // << ", numPos = " << numPos
1410  // << std::endl;
1411  //}
1412 
1413  double tmp = log((double) numPos)/log(2.);
1414  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1415  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1416  unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
1417 
1418  std::vector<double> rawDataVec(numPos,0.);
1419  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1420  Fft<T> fftObj(m_env);
1421 
1422  // Forward FFT
1423  this->extractRawData(initialPos,
1424  1, // spacing
1425  numPos,
1426  rawDataVec);
1427  T meanValue = this->subMeanExtra(initialPos,
1428  numPos);
1429  for (unsigned int j = 0; j < numPos; ++j) {
1430  rawDataVec[j] -= meanValue; // IMPORTANT
1431  }
1432  rawDataVec.resize(fftSize,0.);
1433 
1434  //if (m_env.subDisplayFile()) {
1435  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
1436  // << ": about to call fftObj.forward()"
1437  // << " with rawDataVec.size() = " << rawDataVec.size()
1438  // << ", fftSize = " << fftSize
1439  // << ", resultData.size() = " << resultData.size()
1440  // << std::endl;
1441  //}
1442  fftObj.forward(rawDataVec,fftSize,resultData);
1443 
1444  // Inverse FFT
1445  for (unsigned int j = 0; j < fftSize; ++j) {
1446  rawDataVec[j] = std::norm(resultData[j]);
1447  }
1448  fftObj.inverse(rawDataVec,fftSize,resultData);
1449 
1450  //if (m_env.subDisplayFile()) {
1451  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
1452  // << ": computed auto covariance for lag 0 = " << resultData[0].real()/((double) (numPos))
1453  // << ", computed resultData[0].imag() = " << resultData[0].imag()
1454  // << std::endl;
1455  //}
1456 
1457  // Prepare return data
1458  autoCorrsSum = 0.;
1459  for (unsigned int j = 0; j < numSum; ++j) { // Yes, begin at lag '0'
1460  double ratio = ((double) j)/((double) (numPos-1));
1461  autoCorrsSum += ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1462  }
1463 
1464  return;
1465 }
1466 // --------------------------------------------------
1467 template <class T>
1468 void
1469 ScalarSequence<T>::subMinMaxExtra(
1470  unsigned int initialPos,
1471  unsigned int numPos,
1472  T& minValue,
1473  T& maxValue) const
1474 {
1475  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
1476 
1477  seqScalarPositionConstIteratorTypedef pos1 = m_seq.begin();
1478  std::advance(pos1,initialPos);
1479 
1480  seqScalarPositionConstIteratorTypedef pos2 = m_seq.begin();
1481  std::advance(pos2,initialPos+numPos);
1482 
1483  if ((initialPos+numPos) == this->subSequenceSize()) {
1484  queso_require_msg(!(pos2 != m_seq.end()), "invalid state");
1485  }
1486 
1487  seqScalarPositionConstIteratorTypedef pos;
1488  pos = std::min_element(pos1, pos2);
1489  minValue = *pos;
1490  pos = std::max_element(pos1, pos2);
1491  maxValue = *pos;
1492 
1493  return;
1494 }
1495 // --------------------------------------------------
1496 template <class T>
1497 void
1498 ScalarSequence<T>::unifiedMinMaxExtra(
1499  bool useOnlyInter0Comm,
1500  unsigned int initialPos,
1501  unsigned int numPos,
1502  T& unifiedMinValue,
1503  T& unifiedMaxValue) const
1504 {
1505  if (m_env.numSubEnvironments() == 1) {
1506  return this->subMinMaxExtra(initialPos,
1507  numPos,
1508  unifiedMinValue,
1509  unifiedMaxValue);
1510  }
1511 
1512  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1513 
1514  if (useOnlyInter0Comm) {
1515  if (m_env.inter0Rank() >= 0) {
1516  // Find local min and max
1517  T minValue;
1518  T maxValue;
1519  this->subMinMaxExtra(initialPos,
1520  numPos,
1521  minValue,
1522  maxValue);
1523 
1524  // Get overall min
1525  std::vector<double> sendBuf(1,0.);
1526  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1527  sendBuf[i] = minValue;
1528  }
1529  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMinValue, (int) sendBuf.size(), RawValue_MPI_MIN,
1530  "ScalarSequence<T>::unifiedMinMaxExtra()",
1531  "failed MPI.Allreduce() for min");
1532 
1533  // Get overall max
1534  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1535  sendBuf[i] = maxValue;
1536  }
1537  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMaxValue, (int) sendBuf.size(), RawValue_MPI_MAX,
1538  "ScalarSequence<T>::unifiedMinMaxExtra()",
1539  "failed MPI.Allreduce() for max");
1540 
1541  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1542  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMinMaxExtra()"
1543  << ": localMinValue = " << minValue
1544  << ", localMaxValue = " << maxValue
1545  << ", unifiedMinValue = " << unifiedMinValue
1546  << ", unifiedMaxValue = " << unifiedMaxValue
1547  << std::endl;
1548  }
1549  }
1550  else {
1551  // Node not in the 'inter0' communicator
1552  this->subMinMaxExtra(initialPos,
1553  numPos,
1554  unifiedMinValue,
1555  unifiedMaxValue);
1556  }
1557  }
1558  else {
1559  queso_error_msg("parallel vectors not supported yet");
1560  }
1561 
1562  //m_env.fullComm().Barrier();
1563 
1564  return;
1565 }
1566 // --------------------------------------------------
1567 template <class T>
1568 void
1569 ScalarSequence<T>::subHistogram(
1570  unsigned int initialPos,
1571  const T& minHorizontalValue,
1572  const T& maxHorizontalValue,
1573  std::vector<T>& centers,
1574  std::vector<unsigned int>& bins) const
1575 {
1576  queso_require_equal_to_msg(centers.size(), bins.size(), "vectors 'centers' and 'bins' have different sizes");
1577 
1578  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1579 
1580  if (initialPos) {}; // just to remove compiler warning
1581 
1582  for (unsigned int j = 0; j < bins.size(); ++j) {
1583  centers[j] = 0.;
1584  bins[j] = 0;
1585  }
1586 
1587  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1588 
1589  double minCenter = minHorizontalValue - horizontalDelta/2.;
1590  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1591  for (unsigned int j = 0; j < centers.size(); ++j) {
1592  double factor = ((double) j)/(((double) centers.size()) - 1.);
1593  centers[j] = (1. - factor) * minCenter + factor * maxCenter;
1594  }
1595 
1596  unsigned int dataSize = this->subSequenceSize();
1597  for (unsigned int j = 0; j < dataSize; ++j) {
1598  double value = m_seq[j];
1599  if (value < minHorizontalValue) {
1600  bins[0]++;
1601  }
1602  else if (value >= maxHorizontalValue) {
1603  bins[bins.size()-1]++;
1604  }
1605  else {
1606  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1607  bins[index]++;
1608  }
1609  }
1610 
1611  return;
1612 }
1613 // --------------------------------------------------
1614 template <class T>
1615 void
1616 ScalarSequence<T>::unifiedHistogram(
1617  bool useOnlyInter0Comm,
1618  unsigned int initialPos,
1619  const T& unifiedMinHorizontalValue,
1620  const T& unifiedMaxHorizontalValue,
1621  std::vector<T>& unifiedCenters,
1622  std::vector<unsigned int>& unifiedBins) const
1623 {
1624  if (m_env.numSubEnvironments() == 1) {
1625  return this->subHistogram(initialPos,
1626  unifiedMinHorizontalValue,
1627  unifiedMaxHorizontalValue,
1628  unifiedCenters,
1629  unifiedBins);
1630  }
1631 
1632  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
1633 
1634  if (useOnlyInter0Comm) {
1635  if (m_env.inter0Rank() >= 0) {
1636  queso_require_equal_to_msg(unifiedCenters.size(), unifiedBins.size(), "vectors 'unifiedCenters' and 'unifiedBins' have different sizes");
1637 
1638  queso_require_greater_equal_msg(unifiedBins.size(), 3, "number of 'unifiedBins' is too small: should be at least 3");
1639 
1640  for (unsigned int j = 0; j < unifiedBins.size(); ++j) {
1641  unifiedCenters[j] = 0.;
1642  unifiedBins[j] = 0;
1643  }
1644 
1645  double unifiedHorizontalDelta = (unifiedMaxHorizontalValue - unifiedMinHorizontalValue)/(((double) unifiedBins.size()) - 2.); // IMPORTANT: -2
1646 
1647  double unifiedMinCenter = unifiedMinHorizontalValue - unifiedHorizontalDelta/2.;
1648  double unifiedMaxCenter = unifiedMaxHorizontalValue + unifiedHorizontalDelta/2.;
1649  for (unsigned int j = 0; j < unifiedCenters.size(); ++j) {
1650  double factor = ((double) j)/(((double) unifiedCenters.size()) - 1.);
1651  unifiedCenters[j] = (1. - factor) * unifiedMinCenter + factor * unifiedMaxCenter;
1652  }
1653 
1654  std::vector<unsigned int> localBins(unifiedBins.size(),0);
1655  unsigned int dataSize = this->subSequenceSize();
1656  for (unsigned int j = 0; j < dataSize; ++j) {
1657  double value = m_seq[j];
1658  if (value < unifiedMinHorizontalValue) {
1659  localBins[0]++;
1660  }
1661  else if (value >= unifiedMaxHorizontalValue) {
1662  localBins[localBins.size()-1]++;
1663  }
1664  else {
1665  unsigned int index = 1 + (unsigned int) ((value - unifiedMinHorizontalValue)/unifiedHorizontalDelta);
1666  localBins[index]++;
1667  }
1668  }
1669 
1670  m_env.inter0Comm().template Allreduce<unsigned int>(&localBins[0], &unifiedBins[0], (int) localBins.size(), RawValue_MPI_SUM,
1671  "ScalarSequence<T>::unifiedHistogram()",
1672  "failed MPI.Allreduce() for bins");
1673 
1674  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1675  for (unsigned int i = 0; i < unifiedCenters.size(); ++i) {
1676  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedHistogram()"
1677  << ": i = " << i
1678  << ", unifiedMinHorizontalValue = " << unifiedMinHorizontalValue
1679  << ", unifiedMaxHorizontalValue = " << unifiedMaxHorizontalValue
1680  << ", unifiedCenters = " << unifiedCenters[i]
1681  << ", unifiedBins = " << unifiedBins[i]
1682  << std::endl;
1683  }
1684  }
1685  }
1686  else {
1687  // Node not in the 'inter0' communicator
1688  this->subHistogram(initialPos,
1689  unifiedMinHorizontalValue,
1690  unifiedMaxHorizontalValue,
1691  unifiedCenters,
1692  unifiedBins);
1693  }
1694  }
1695  else {
1696  queso_error_msg("parallel vectors not supported yet");
1697  }
1698 
1699  //m_env.fullComm().Barrier();
1700 
1701  return;
1702 }
1703 // --------------------------------------------------
1704 template <class T>
1705 void
1706 ScalarSequence<T>::subBasicHistogram(
1707  unsigned int initialPos,
1708  const T& minHorizontalValue,
1709  const T& maxHorizontalValue,
1710  UniformOneDGrid<T>*& gridValues,
1711  std::vector<unsigned int>& bins) const
1712 {
1713  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1714 
1715  for (unsigned int j = 0; j < bins.size(); ++j) {
1716  bins[j] = 0;
1717  }
1718 
1719  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1720  double minCenter = minHorizontalValue - horizontalDelta/2.;
1721  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1722  gridValues = new UniformOneDGrid<T>(m_env,
1723  "",
1724  bins.size(),
1725  minCenter,
1726  maxCenter);
1727 
1728  unsigned int dataSize = this->subSequenceSize();
1729  for (unsigned int j = 0; j < dataSize; ++j) {
1730  double value = m_seq[j];
1731  if (value < minHorizontalValue) {
1732  bins[0]++;
1733  }
1734  else if (value >= maxHorizontalValue) {
1735  bins[bins.size()-1]++;
1736  }
1737  else {
1738  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1739  bins[index]++;
1740  }
1741  }
1742 
1743  return;
1744 }
1745 // --------------------------------------------------
1746 template <class T>
1747 void
1748 ScalarSequence<T>::subWeightHistogram(
1749  unsigned int initialPos,
1750  const T& minHorizontalValue,
1751  const T& maxHorizontalValue,
1752  UniformOneDGrid<T>*& gridValues,
1753  std::vector<unsigned int>& bins) const
1754 {
1755  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1756 
1757  for (unsigned int j = 0; j < bins.size(); ++j) {
1758  bins[j] = 0;
1759  }
1760 
1761  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1762  double minCenter = minHorizontalValue - horizontalDelta/2.;
1763  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1764  gridValues = new UniformOneDGrid<T>(m_env,
1765  "",
1766  bins.size(),
1767  minCenter,
1768  maxCenter);
1769 
1770  unsigned int dataSize = this->subSequenceSize();
1771  for (unsigned int j = 0; j < dataSize; ++j) {
1772  double value = m_seq[j];
1773  if (value < minHorizontalValue) {
1774  bins[0]++;
1775  }
1776  else if (value >= maxHorizontalValue) {
1777  bins[bins.size()-1]++;
1778  }
1779  else {
1780  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1781  bins[index]++;
1782  }
1783  }
1784 
1785  return;
1786 }
1787 // --------------------------------------------------
1788 template <class T>
1789 void
1790 ScalarSequence<T>::subWeightHistogram(
1791  unsigned int initialPos,
1792  const T& minHorizontalValue,
1793  const T& maxHorizontalValue,
1794  std::vector<T>& gridValues,
1795  std::vector<unsigned int>& bins) const
1796 {
1797  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1798 
1799  for (unsigned int j = 0; j < bins.size(); ++j) {
1800  bins[j] = 0;
1801  }
1802 
1803  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1804  double minCenter = minHorizontalValue - horizontalDelta/2.;
1805  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1806  UniformOneDGrid<T> tmpGrid(m_env,
1807  "",
1808  bins.size(),
1809  minCenter,
1810  maxCenter);
1811  gridValues.clear();
1812  gridValues.resize(tmpGrid.size(),0.);
1813  for (unsigned int i = 0; i < tmpGrid.size(); ++i) {
1814  gridValues[i] = tmpGrid[i];
1815  }
1816 
1817  unsigned int dataSize = this->subSequenceSize();
1818  for (unsigned int j = 0; j < dataSize; ++j) {
1819  double value = m_seq[j];
1820  if (value < minHorizontalValue) {
1821  bins[0]++;
1822  }
1823  else if (value >= maxHorizontalValue) {
1824  bins[bins.size()-1]++;
1825  }
1826  else {
1827  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1828  bins[index]++;
1829  }
1830  }
1831 
1832  //std::cout << "In ScalarSequence<T>::subWeightHistogram():" << std::endl;
1833  //for (unsigned int j = 0; j < bins.size(); ++j) {
1834  // std::cout << "bins[" << j << "] = " << bins[j] << std::endl;
1835  //}
1836 
1837  return;
1838 }
1839 // --------------------------------------------------
1840 template <class T>
1841 void
1842 ScalarSequence<T>::subSort(
1843  unsigned int initialPos,
1844  ScalarSequence<T>& sortedSequence) const
1845 {
1846  unsigned int numPos = this->subSequenceSize() - initialPos;
1847  sortedSequence.resizeSequence(numPos);
1848  this->extractScalarSeq(initialPos,
1849  1,
1850  numPos,
1851  sortedSequence);
1852  sortedSequence.subSort();
1853 
1854  return;
1855 }
1856 // --------------------------------------------------
1857 template <class T>
1858 void
1859 ScalarSequence<T>::unifiedSort(
1860  bool useOnlyInter0Comm,
1861  unsigned int initialPos,
1862  ScalarSequence<T>& unifiedSortedSequence) const
1863 {
1864  if (m_env.numSubEnvironments() == 1) {
1865  return this->subSort(initialPos,unifiedSortedSequence);
1866  }
1867 
1868  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
1869 
1870  if (useOnlyInter0Comm) {
1871  if (m_env.inter0Rank() >= 0) {
1872  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1873 
1874  unsigned int localNumPos = this->subSequenceSize() - initialPos;
1875 
1876  std::vector<T> leafData(localNumPos,0.);
1877  this->extractRawData(0,
1878  1,
1879  localNumPos,
1880  leafData);
1881 
1882  if (m_env.inter0Rank() == 0) {
1883  int minus1NumTreeLevels = 0;
1884  int power2NumTreeNodes = 1;
1885 
1886  while (power2NumTreeNodes < m_env.inter0Comm().NumProc()) {
1887  power2NumTreeNodes += power2NumTreeNodes;
1888  minus1NumTreeLevels++;
1889  }
1890 
1891  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1892  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedSort()"
1893  << ": sorting tree has " << m_env.inter0Comm().NumProc()
1894  << " nodes and " << minus1NumTreeLevels+1
1895  << " levels"
1896  << std::endl;
1897  }
1898 
1899  this->parallelMerge(unifiedSortedSequence.rawData(),
1900  leafData,
1901  minus1NumTreeLevels);
1902  }
1903  else if (m_env.inter0Rank() > 0) { // KAUST
1904  unsigned int uintBuffer[1];
1905  RawType_MPI_Status status;
1906  m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, RawValue_MPI_ANY_SOURCE, SCALAR_SEQUENCE_INIT_MPI_MSG, &status,
1907  "ScalarSequence<T>::unifiedSort()",
1908  "failed MPI.Recv() for init");
1909  //if (status) {}; // just to remove compiler warning
1910 
1911  unsigned int treeLevel = uintBuffer[0];
1912  this->parallelMerge(unifiedSortedSequence.rawData(),
1913  leafData,
1914  treeLevel);
1915  }
1916 
1917  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), returned from parallelMerge()",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1918 
1919  // Broadcast
1920  unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
1921  m_env.inter0Comm().Bcast((void *) &unifiedDataSize, (int) 1, RawValue_MPI_UNSIGNED, 0,
1922  "ScalarSequence<T>::unifiedSort()",
1923  "failed MPI.Bcast() for unified data size");
1924 
1925  unsigned int sumOfNumPos = 0;
1926  m_env.inter0Comm().template Allreduce<unsigned int>(&localNumPos, &sumOfNumPos, (int) 1, RawValue_MPI_SUM,
1927  "ScalarSequence<T>::unifiedSort()",
1928  "failed MPI.Allreduce() for data size");
1929 
1930  queso_require_equal_to_msg(sumOfNumPos, unifiedDataSize, "incompatible unified sizes");
1931 
1932  unifiedSortedSequence.resizeSequence(unifiedDataSize);
1933  m_env.inter0Comm().Bcast((void *) &unifiedSortedSequence.rawData()[0], (int) unifiedDataSize, RawValue_MPI_DOUBLE, 0,
1934  "ScalarSequence<T>::unifiedSort()",
1935  "failed MPI.Bcast() for unified data");
1936 
1937  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1938  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
1939  << ": tree node " << m_env.inter0Rank()
1940  << ", unifiedSortedSequence[0] = " << unifiedSortedSequence[0]
1941  << ", unifiedSortedSequence[" << unifiedSortedSequence.subSequenceSize()-1 << "] = " << unifiedSortedSequence[unifiedSortedSequence.subSequenceSize()-1]
1942  << std::endl;
1943  }
1944 
1945  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1946  }
1947  else {
1948  // Node not in the 'inter0' communicator
1949  this->subSort(initialPos,unifiedSortedSequence);
1950  }
1951  }
1952  else {
1953  queso_error_msg("parallel vectors not supported yet");
1954  }
1955 
1956  //m_env.fullComm().Barrier();
1957 
1958  return;
1959 }
1960 // --------------------------------------------------
1961 template <class T>
1962 T
1963 ScalarSequence<T>::subInterQuantileRange(unsigned int initialPos) const
1964 {
1965  queso_require_less_msg(initialPos, this->subSequenceSize(), "'initialPos' is too big");
1966 
1967  ScalarSequence sortedSequence(m_env,0,"");
1968  this->subSort(initialPos,
1969  sortedSequence);
1970 
1971  // The test above guarantees that 'dataSize >= 1'
1972  unsigned int dataSize = this->subSequenceSize() - initialPos;
1973 
1974  queso_require_equal_to_msg(dataSize, sortedSequence.subSequenceSize(), "inconsistent size variables");
1975 
1976  bool everythingOk = true;
1977 
1978  // pos1 = (dataSize+1)/4 - 1
1979  // pos1 >= 0 <==> dataSize >= 3
1980  // pos1 < (dataSize-1) <==> 3*dataSize > 1
1981  unsigned int pos1 = (unsigned int) ( (((double) dataSize) + 1.)*1./4. - 1. );
1982  if (pos1 > (dataSize-1)) {
1983  pos1 = 0;
1984  everythingOk = false;
1985  }
1986  unsigned int pos1inc = pos1+1;
1987  if (pos1inc > (dataSize-1)) {
1988  pos1inc = dataSize-1;
1989  everythingOk = false;
1990  }
1991 
1992  // pos3 = (dataSize+1)*3/4 - 1
1993  // pos3 >= 0 <==> dataSize >= 1/3
1994  // pos3 < (dataSize-1) <==> dataSize > 3
1995  unsigned int pos3 = (unsigned int) ( (((double) dataSize) + 1.)*3./4. - 1. );
1996  if (pos3 > (dataSize-1)) {
1997  pos3 = 0;
1998  everythingOk = false;
1999  }
2000  unsigned int pos3inc = pos3+1;
2001  if (pos3inc > (dataSize-1)) {
2002  pos3inc = dataSize-1;
2003  everythingOk = false;
2004  }
2005 
2006  double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
2007  if (fraction1 < 0.) {
2008  fraction1 = 0.;
2009  everythingOk = false;
2010  }
2011  double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
2012  if (fraction3 < 0.) {
2013  fraction3 = 0.;
2014  everythingOk = false;
2015  }
2016 
2017  if (everythingOk == false) {
2018  std::cerr << "In ScalarSequence<T>::subInterQuantileRange()"
2019  << ", worldRank = " << m_env.worldRank()
2020  << ": at least one adjustment was necessary"
2021  << std::endl;
2022  }
2023 
2024  //if (m_env.subDisplayFile()) {
2025  // *m_env.subDisplayFile() << "In ScalarSequence::subInterQuantileRange()"
2026  // << ", initialPos = " << initialPos
2027  // << ", this->subSequenceSize() = " << this->subSequenceSize()
2028  // << ", dataSize = " << dataSize
2029  // << ", sortedSequence.size() = " << sortedSequence.size()
2030  // << ", pos1 = " << pos1
2031  // << ", pos3 = " << pos3
2032  // << std::endl;
2033  //}
2034 
2035  T value1 = (1.-fraction1) * sortedSequence[pos1] + fraction1 * sortedSequence[pos1inc];
2036  T value3 = (1.-fraction3) * sortedSequence[pos3] + fraction3 * sortedSequence[pos3inc];
2037  T iqrValue = value3 - value1;
2038 
2039  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2040  *m_env.subDisplayFile() << "In ScalarSequence<T>::subInterQuantileRange()"
2041  << ": iqrValue = " << iqrValue
2042  << ", dataSize = " << dataSize
2043  << ", pos1 = " << pos1
2044  << ", pos3 = " << pos3
2045  << ", value1 = " << value1
2046  << ", value3 = " << value3
2047  << std::endl;
2048 
2049  // Save data only once into a separate file
2050  //std::ofstream* ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
2051  //if ((ofsvar == NULL ) ||
2052  // (ofsvar->is_open() == false)) {
2053  // delete ofsvar;
2054  // ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
2055 
2056  // *ofsvar << "var_sort_sub" << m_env.subIdString() << " = zeros(" << 1
2057  // << "," << dataSize
2058  // << ");"
2059  // << std::endl;
2060  // for (unsigned int j = 0; j < dataSize; ++j) {
2061  // *ofsvar << "var_sort_sub" << m_env.subIdString() << "(" << 1
2062  // << "," << j+1
2063  // << ") = " << sortedSequence[j]
2064  // << ";"
2065  // << std::endl;
2066  // }
2067  //}
2068  //delete ofsvar;
2069  }
2070 
2071  return iqrValue;
2072 }
2073 // --------------------------------------------------
2074 template <class T>
2075 T
2076 ScalarSequence<T>::unifiedInterQuantileRange(
2077  bool useOnlyInter0Comm,
2078  unsigned int initialPos) const
2079 {
2080  T unifiedIqrValue = 0.;
2081 
2082  if (m_env.numSubEnvironments() == 1) {
2083  return this->subInterQuantileRange(initialPos);
2084  }
2085 
2086  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2087 
2088  if (useOnlyInter0Comm) {
2089  if (m_env.inter0Rank() >= 0) {
2090  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
2091 
2092  ScalarSequence unifiedSortedSequence(m_env,0,"");
2093  this->unifiedSort(useOnlyInter0Comm,
2094  initialPos,
2095  unifiedSortedSequence);
2096  unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
2097 
2098  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2099  unsigned int sumOfLocalSizes = 0;
2100  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &sumOfLocalSizes, (int) 1, RawValue_MPI_SUM,
2101  "ScalarSequence<T>::unifiedInterQuantileRange()",
2102  "failed MPI.Allreduce() for data size");
2103 
2104  queso_require_equal_to_msg(sumOfLocalSizes, unifiedDataSize, "incompatible unified sizes");
2105 
2106  unsigned int pos1 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*1./4. - 1. );
2107  unsigned int pos3 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*3./4. - 1. );
2108 
2109  double fraction1 = (((double) unifiedDataSize) + 1.)*1./4. - 1. - ((double) pos1);
2110  double fraction3 = (((double) unifiedDataSize) + 1.)*3./4. - 1. - ((double) pos3);
2111 
2112  T value1 = (1.-fraction1) * unifiedSortedSequence[pos1] + fraction1 * unifiedSortedSequence[pos1+1];
2113  T value3 = (1.-fraction3) * unifiedSortedSequence[pos3] + fraction3 * unifiedSortedSequence[pos3+1];
2114  unifiedIqrValue = value3 - value1;
2115 
2116  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2117  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedInterQuantileRange()"
2118  << ": unifiedIqrValue = " << unifiedIqrValue
2119  << ", localDataSize = " << localDataSize
2120  << ", unifiedDataSize = " << unifiedDataSize
2121  << ", pos1 = " << pos1
2122  << ", pos3 = " << pos3
2123  << ", value1 = " << value1
2124  << ", value3 = " << value3
2125  << std::endl;
2126 
2127  // Save data only once into a separate file
2128  //std::ofstream* ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
2129  //if ((ofsvar == NULL ) ||
2130  // (ofsvar->is_open() == false)) {
2131  // delete ofsvar;
2132  // ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
2133 
2134  // *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << " = zeros(" << 1
2135  // << "," << unifiedDataSize
2136  // << ");"
2137  // << std::endl;
2138  // for (unsigned int j = 0; j < unifiedDataSize; ++j) {
2139  // *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << "(" << 1
2140  // << "," << j+1
2141  // << ") = " << unifiedSortedSequence[j]
2142  // << ";"
2143  // << std::endl;
2144  // }
2145  //}
2146  //delete ofsvar;
2147  }
2148 
2149  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
2150  }
2151  else {
2152  // Node not in the 'inter0' communicator
2153  unifiedIqrValue = this->subInterQuantileRange(initialPos);
2154  }
2155  }
2156  else {
2157  queso_error_msg("parallel vectors not supported yet");
2158  }
2159 
2160  //m_env.fullComm().Barrier();
2161 
2162  return unifiedIqrValue;
2163 }
2164 // --------------------------------------------------
2165 template <class T>
2166 T
2167 ScalarSequence<T>::subScaleForKde(
2168  unsigned int initialPos,
2169  const T& iqrValue,
2170  unsigned int kdeDimension) const
2171 {
2172  bool bRC = (initialPos < this->subSequenceSize());
2173  queso_require_msg(bRC, "invalid input data");
2174 
2175  unsigned int dataSize = this->subSequenceSize() - initialPos;
2176 
2177  T meanValue = this->subMeanExtra(initialPos,
2178  dataSize);
2179 
2180  T samValue = this->subSampleVarianceExtra(initialPos,
2181  dataSize,
2182  meanValue);
2183 
2184  T scaleValue;
2185  if (iqrValue <= 0.) {
2186  scaleValue = 1.06*std::sqrt(samValue)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2187  }
2188  else {
2189  scaleValue = 1.06*std::min(std::sqrt(samValue),iqrValue/1.34)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2190  }
2191 
2192  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2193  *m_env.subDisplayFile() << "In ScalarSequence<T>::subScaleForKde()"
2194  << ": iqrValue = " << iqrValue
2195  << ", meanValue = " << meanValue
2196  << ", samValue = " << samValue
2197  << ", dataSize = " << dataSize
2198  << ", scaleValue = " << scaleValue
2199  << std::endl;
2200  }
2201 
2202  return scaleValue;
2203 }
2204 // --------------------------------------------------
2205 template <class T>
2206 T
2207 ScalarSequence<T>::unifiedScaleForKde(
2208  bool useOnlyInter0Comm,
2209  unsigned int initialPos,
2210  const T& unifiedIqrValue,
2211  unsigned int kdeDimension) const
2212 {
2213  if (m_env.numSubEnvironments() == 1) {
2214  return this->subScaleForKde(initialPos,
2215  unifiedIqrValue,
2216  kdeDimension);
2217  }
2218 
2219  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2220 
2221  T unifiedScaleValue = 0.;
2222  if (useOnlyInter0Comm) {
2223  if (m_env.inter0Rank() >= 0) {
2224  bool bRC = (initialPos < this->subSequenceSize());
2225  queso_require_msg(bRC, "invalid input data");
2226 
2227  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2228 
2229  T unifiedMeanValue = this->unifiedMeanExtra(useOnlyInter0Comm,
2230  initialPos,
2231  localDataSize);
2232 
2233  T unifiedSamValue = this->unifiedSampleVarianceExtra(useOnlyInter0Comm,
2234  initialPos,
2235  localDataSize,
2236  unifiedMeanValue);
2237 
2238  unsigned int unifiedDataSize = 0;
2239  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &unifiedDataSize, (int) 1, RawValue_MPI_SUM,
2240  "ScalarSequence<T>::unifiedScaleForKde()",
2241  "failed MPI.Allreduce() for data size");
2242 
2243  if (unifiedIqrValue <= 0.) {
2244  unifiedScaleValue = 1.06*std::sqrt(unifiedSamValue)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
2245  }
2246  else {
2247  unifiedScaleValue = 1.06*std::min(std::sqrt(unifiedSamValue),unifiedIqrValue/1.34)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
2248  }
2249 
2250  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2251  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedScaleForKde()"
2252  << ": unifiedIqrValue = " << unifiedIqrValue
2253  << ", unifiedMeanValue = " << unifiedMeanValue
2254  << ", unifiedSamValue = " << unifiedSamValue
2255  << ", unifiedDataSize = " << unifiedDataSize
2256  << ", unifiedScaleValue = " << unifiedScaleValue
2257  << std::endl;
2258  }
2259  }
2260  else {
2261  // Node not in the 'inter0' communicator
2262  unifiedScaleValue = this->subScaleForKde(initialPos,
2263  unifiedIqrValue,
2264  kdeDimension);
2265  }
2266  }
2267  else {
2268  queso_error_msg("parallel vectors not supported yet");
2269  }
2270 
2271  //m_env.fullComm().Barrier();
2272 
2273  return unifiedScaleValue;
2274 }
2275 // --------------------------------------------------
2276 template <class T>
2277 void
2278 ScalarSequence<T>::subGaussian1dKde(
2279  unsigned int initialPos,
2280  double scaleValue,
2281  const std::vector<T>& evaluationPositions,
2282  std::vector<double>& densityValues) const
2283 {
2284  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2285  (0 < evaluationPositions.size()) &&
2286  (evaluationPositions.size() == densityValues.size() ));
2287  queso_require_msg(bRC, "invalid input data");
2288 
2289  unsigned int dataSize = this->subSequenceSize() - initialPos;
2290  unsigned int numEvals = evaluationPositions.size();
2291 
2292  double scaleInv = 1./scaleValue;
2293  for (unsigned int j = 0; j < numEvals; ++j) {
2294  double x = evaluationPositions[j];
2295  double value = 0.;
2296  for (unsigned int k = 0; k < dataSize; ++k) {
2297  double xk = m_seq[initialPos+k];
2298  value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
2299  }
2300  densityValues[j] = scaleInv * (value/(double) dataSize);
2301  }
2302 
2303  return;
2304 }
2305 // --------------------------------------------------
2306 template <class T>
2307 void
2308 ScalarSequence<T>::unifiedGaussian1dKde(
2309  bool useOnlyInter0Comm,
2310  unsigned int initialPos,
2311  double unifiedScaleValue,
2312  const std::vector<T>& unifiedEvaluationPositions,
2313  std::vector<double>& unifiedDensityValues) const
2314 {
2315  if (m_env.numSubEnvironments() == 1) {
2316  return this->subGaussian1dKde(initialPos,
2317  unifiedScaleValue,
2318  unifiedEvaluationPositions,
2319  unifiedDensityValues);
2320  }
2321 
2322  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2323 
2324  if (useOnlyInter0Comm) {
2325  if (m_env.inter0Rank() >= 0) {
2326  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2327  (0 < unifiedEvaluationPositions.size()) &&
2328  (unifiedEvaluationPositions.size() == unifiedDensityValues.size() ));
2329  queso_require_msg(bRC, "invalid input data");
2330 
2331  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2332  unsigned int unifiedDataSize = 0;
2333  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &unifiedDataSize, (int) 1, RawValue_MPI_SUM,
2334  "ScalarSequence<T>::unifiedGaussian1dKde()",
2335  "failed MPI.Allreduce() for data size");
2336 
2337  unsigned int numEvals = unifiedEvaluationPositions.size();
2338 
2339  std::vector<double> densityValues(numEvals,0.);
2340  double unifiedScaleInv = 1./unifiedScaleValue;
2341  for (unsigned int j = 0; j < numEvals; ++j) {
2342  double x = unifiedEvaluationPositions[j];
2343  double value = 0.;
2344  for (unsigned int k = 0; k < localDataSize; ++k) {
2345  double xk = m_seq[initialPos+k];
2346  value += MiscGaussianDensity((x-xk)*unifiedScaleInv,0.,1.);
2347  }
2348  densityValues[j] = value;
2349  }
2350 
2351  for (unsigned int j = 0; j < numEvals; ++j) {
2352  unifiedDensityValues[j] = 0.;
2353  }
2354  m_env.inter0Comm().template Allreduce<double>(&densityValues[0], &unifiedDensityValues[0], (int) numEvals, RawValue_MPI_SUM,
2355  "ScalarSequence<T>::unifiedGaussian1dKde()",
2356  "failed MPI.Allreduce() for density values");
2357 
2358  for (unsigned int j = 0; j < numEvals; ++j) {
2359  unifiedDensityValues[j] *= unifiedScaleInv/((double) unifiedDataSize);
2360  }
2361 
2362  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2363  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedGaussian1dKde()"
2364  << ": unifiedDensityValues[0] = " << unifiedDensityValues[0]
2365  << ", unifiedDensityValues[" << unifiedDensityValues.size()-1 << "] = " << unifiedDensityValues[unifiedDensityValues.size()-1]
2366  << std::endl;
2367  }
2368  }
2369  else {
2370  // Node not in the 'inter0' communicator
2371  this->subGaussian1dKde(initialPos,
2372  unifiedScaleValue,
2373  unifiedEvaluationPositions,
2374  unifiedDensityValues);
2375  }
2376  }
2377  else {
2378  queso_error_msg("parallel vectors not supported yet");
2379  }
2380 
2381  //m_env.fullComm().Barrier();
2382 
2383  return;
2384 }
2385 // --------------------------------------------------
2386 template <class T>
2387 void
2388 ScalarSequence<T>::filter(
2389  unsigned int initialPos,
2390  unsigned int spacing)
2391 {
2392  if (m_env.subDisplayFile()) {
2393  *m_env.subDisplayFile() << "Entering ScalarSequence<V,M>::filter()"
2394  << ": initialPos = " << initialPos
2395  << ", spacing = " << spacing
2396  << ", subSequenceSize = " << this->subSequenceSize()
2397  << std::endl;
2398  }
2399 
2400  unsigned int i = 0;
2401  unsigned int j = initialPos;
2402  unsigned int originalSubSequenceSize = this->subSequenceSize();
2403  while (j < originalSubSequenceSize) {
2404  if (i != j) {
2405  //*m_env.subDisplayFile() << i << "--" << j << " ";
2406  m_seq[i] = m_seq[j];
2407  }
2408  i++;
2409  j += spacing;
2410  }
2411 
2412  this->resizeSequence(i);
2413 
2414  if (m_env.subDisplayFile()) {
2415  *m_env.subDisplayFile() << "Leaving ScalarSequence<V,M>::filter()"
2416  << ": initialPos = " << initialPos
2417  << ", spacing = " << spacing
2418  << ", subSequenceSize = " << this->subSequenceSize()
2419  << std::endl;
2420  }
2421 
2422  return;
2423 }
2424 // --------------------------------------------------
2425 template <class T>
2426 T
2427 ScalarSequence<T>::brooksGelmanConvMeasure(
2428  bool useOnlyInter0Comm,
2429  unsigned int initialPos,
2430  unsigned int spacing) const
2431 {
2432  double resultValue = 0.;
2433 
2434  // FIX ME: put logic if (numSubEnvs == 1) ...
2435 
2436  if (useOnlyInter0Comm) {
2437  if (m_env.inter0Rank() >= 0) {
2438  queso_not_implemented();
2439  }
2440  else {
2441  // Node not in the 'inter0' communicator
2442  // Do nothing
2443  }
2444  }
2445  else {
2446  queso_error_msg("parallel vectors not supported yet");
2447  }
2448 
2449  //m_env.fullComm().Barrier();
2450 
2451  return resultValue;
2452 }
2453 // --------------------------------------------------
2454 template <class T>
2455 void
2456 ScalarSequence<T>::append(
2457  const ScalarSequence<T>& src,
2458  unsigned int srcInitialPos,
2459  unsigned int srcNumPos)
2460 {
2461  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+1), "srcInitialPos is too big");
2462 
2463  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+srcNumPos), "srcNumPos is too big");
2464 
2465  deleteStoredScalars();
2466  unsigned int currentSize = this->subSequenceSize();
2467  m_seq.resize(currentSize+srcNumPos,0.);
2468  for (unsigned int i = 0; i < srcNumPos; ++i) {
2469  m_seq[currentSize+i] = src.m_seq[srcInitialPos+i];
2470  }
2471 
2472  return;
2473 }
2474 // --------------------------------------------------
2475 template <class T>
2476 T
2477 ScalarSequence<T>::subPositionsOfMaximum(
2478  const ScalarSequence<T>& subCorrespondingScalarValues,
2479  ScalarSequence<T>& subPositionsOfMaximum)
2480 {
2481  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2482 
2483  T subMaxValue = subCorrespondingScalarValues.subMaxPlain();
2484  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2485 
2486  unsigned int subNumPos = 0;
2487  for (unsigned int i = 0; i < iMax; ++i) {
2488  if (subCorrespondingScalarValues[i] == subMaxValue) {
2489  subNumPos++;
2490  }
2491  }
2492 
2493  subPositionsOfMaximum.resizeSequence(subNumPos);
2494  unsigned int j = 0;
2495  for (unsigned int i = 0; i < iMax; ++i) {
2496  if (subCorrespondingScalarValues[i] == subMaxValue) {
2497  subPositionsOfMaximum[j] = (*this)[i];
2498  j++;
2499  }
2500  }
2501 
2502  return subMaxValue;
2503 }
2504 // --------------------------------------------------
2505 template <class T>
2506 T
2507 ScalarSequence<T>::unifiedPositionsOfMaximum( // rr0
2508  const ScalarSequence<T>& subCorrespondingScalarValues,
2509  ScalarSequence<T>& unifiedPositionsOfMaximum)
2510 {
2511  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2512 
2513  T maxValue = subCorrespondingScalarValues.subMaxPlain();
2514  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2515 
2516  unsigned int numPos = 0;
2517  for (unsigned int i = 0; i < iMax; ++i) {
2518  if (subCorrespondingScalarValues[i] == maxValue) {
2519  numPos++;
2520  }
2521  }
2522 
2523  unifiedPositionsOfMaximum.resizeSequence(numPos);
2524  unsigned int j = 0;
2525  for (unsigned int i = 0; i < iMax; ++i) {
2526  if (subCorrespondingScalarValues[i] == maxValue) {
2527  unifiedPositionsOfMaximum[j] = (*this)[i];
2528  j++;
2529  }
2530  }
2531 
2532  return maxValue;
2533 }
2534 // --------------------------------------------------
2535 template <class T>
2536 void
2537 ScalarSequence<T>::subWriteContents(
2538  unsigned int initialPos,
2539  unsigned int numPos,
2540  const std::string& fileName,
2541  const std::string& fileType,
2542  const std::set<unsigned int>& allowedSubEnvIds) const
2543 {
2544  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
2545 
2546  FilePtrSetStruct filePtrSet;
2547  if (m_env.openOutputFile(fileName,
2548  fileType,
2549  allowedSubEnvIds,
2550  false, // A 'true' causes problems when the user chooses (via options
2551  // in the input file) to use just one file for all outputs.
2552  filePtrSet)) {
2553 
2554  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT ||
2555  fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2556  this->subWriteContents(initialPos,
2557  numPos,
2558  *filePtrSet.ofsVar,
2559  fileType);
2560  }
2561 #ifdef QUESO_HAS_HDF5
2562  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2563 
2564  // Create dataspace
2565  hsize_t dims[1] = { this->subSequenceSize() };
2566  hid_t dataspace_id = H5Screate_simple(1, dims, dims);
2567  queso_require_greater_equal_msg(
2568  dataspace_id,
2569  0,
2570  "error create dataspace with id: " << dataspace_id);
2571 
2572  // Create dataset
2573  hid_t dataset_id = H5Dcreate(filePtrSet.h5Var,
2574  "data",
2575  H5T_IEEE_F64LE,
2576  dataspace_id,
2577  H5P_DEFAULT,
2578  H5P_DEFAULT,
2579  H5P_DEFAULT);
2580 
2581  queso_require_greater_equal_msg(
2582  dataset_id,
2583  0,
2584  "error creating dataset with id: " << dataset_id);
2585 
2586  // Write the dataset
2587  herr_t status = H5Dwrite(
2588  dataset_id,
2589  H5T_NATIVE_DOUBLE, // The type in memory
2590  H5S_ALL, // The dataspace in memory
2591  dataspace_id, // The file dataspace
2592  H5P_DEFAULT, // Xfer property list
2593  &m_seq[0]);
2594 
2595  queso_require_greater_equal_msg(
2596  status,
2597  0,
2598  "error writing dataset to file with id: " << filePtrSet.h5Var);
2599 
2600  // Clean up
2601  H5Dclose(dataset_id);
2602  H5Sclose(dataspace_id);
2603  }
2604 #endif
2605 
2606  m_env.closeFile(filePtrSet,fileType);
2607  }
2608 }
2609 // --------------------------------------------------
2610 template <class T>
2611 void
2612 ScalarSequence<T>::subWriteContents( // rr0
2613  unsigned int /* initialPos */,
2614  unsigned int /* numPos */,
2615  std::ofstream& ofs,
2616  const std::string& fileType) const
2617 {
2618  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2619  this->writeSubMatlabHeader(ofs, this->subSequenceSize());
2620  }
2621  else if (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2622  this->writeTxtHeader(ofs, this->subSequenceSize());
2623  }
2624 
2625  unsigned int chainSize = this->subSequenceSize();
2626  for (unsigned int j = 0; j < chainSize; ++j) {
2627  ofs << m_seq[j]
2628  << std::endl;
2629  }
2630 
2631  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2632  ofs << "];\n";
2633  }
2634 
2635  return;
2636 }
2637 // --------------------------------------------------
2638 template <class T>
2639 void
2640 ScalarSequence<T>::unifiedWriteContents(
2641  const std::string& fileName,
2642  const std::string& inputFileType) const
2643 {
2644  std::string fileType(inputFileType);
2645 #ifdef QUESO_HAS_HDF5
2646  // Do nothing
2647 #else
2648  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2649  if (m_env.subDisplayFile()) {
2650  *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
2651  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2652  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2653  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2654  << "' instead..."
2655  << std::endl;
2656  }
2657  if (m_env.subRank() == 0) {
2658  std::cerr << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
2659  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2660  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2661  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2662  << "' instead..."
2663  << std::endl;
2664  }
2665  fileType = UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT;
2666  }
2667 #endif
2668 
2669  // All processors in 'fullComm' should call this routine...
2670 
2671  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2672  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
2673  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedWriteContents()"
2674  << ": worldRank " << m_env.worldRank()
2675  << ", subEnvironment " << m_env.subId()
2676  << ", subRank " << m_env.subRank()
2677  << ", inter0Rank " << m_env.inter0Rank()
2678  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
2679  << ", fileName = " << fileName
2680  << ", fileType = " << fileType
2681  << std::endl;
2682  }
2683 
2684  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
2685 
2686  if (m_env.inter0Rank() >= 0) {
2687  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
2688  if (m_env.inter0Rank() == (int) r) {
2689  // My turn
2690  FilePtrSetStruct unifiedFilePtrSet;
2691  // bool writeOver = (r == 0);
2692  bool writeOver = false; // A 'true' causes problems when the user chooses (via options
2693  // in the input file) to use just one file for all outputs.
2694  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 000 \n" << std::endl;
2695  if (m_env.openUnifiedOutputFile(fileName,
2696  fileType, // "m or hdf"
2697  writeOver,
2698  unifiedFilePtrSet)) {
2699  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 001 \n" << std::endl;
2700 
2701  unsigned int chainSize = this->subSequenceSize();
2702  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2703  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2704 
2705  // Write header info
2706  if (r == 0) {
2707  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2708  this->writeUnifiedMatlabHeader(*unifiedFilePtrSet.ofsVar,
2709  this->subSequenceSize()*m_env.inter0Comm().NumProc());
2710  }
2711  else { // It's definitely txt if we get in here
2712  this->writeTxtHeader(*unifiedFilePtrSet.ofsVar,
2713  this->subSequenceSize()*m_env.inter0Comm().NumProc());
2714  }
2715  }
2716 
2717  for (unsigned int j = 0; j < chainSize; ++j) {
2718  *unifiedFilePtrSet.ofsVar << m_seq[j]
2719  << std::endl;
2720  }
2721 
2722  m_env.closeFile(unifiedFilePtrSet,fileType);
2723  }
2724 #ifdef QUESO_HAS_HDF5
2725  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2726  unsigned int numParams = 1; // m_vectorSpace.dimLocal();
2727  if (r == 0) {
2728  hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
2729  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type created" << std::endl;
2730  hsize_t dimsf[1];
2731  dimsf[0] = chainSize;
2732  hid_t dataspace = H5Screate_simple(1, dimsf, NULL); // HDF5_rank = 2
2733  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space created" << std::endl;
2734  hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var,
2735  "data",
2736  datatype,
2737  dataspace,
2738  H5P_DEFAULT, // Link creation property list
2739  H5P_DEFAULT, // Dataset creation property list
2740  H5P_DEFAULT); // Dataset access property list
2741  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set created" << std::endl;
2742 
2743  struct timeval timevalBegin;
2744  int iRC = UQ_OK_RC;
2745  iRC = gettimeofday(&timevalBegin,NULL);
2746  if (iRC) {}; // just to remove compiler warning
2747 
2748  herr_t status;
2749  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 002 \n" << std::endl;
2750  status = H5Dwrite(dataset,
2751  H5T_NATIVE_DOUBLE,
2752  H5S_ALL,
2753  H5S_ALL,
2754  H5P_DEFAULT,
2755  &m_seq[0]);
2756  if (status) {}; // just to remove compiler warning
2757 
2758  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 003 \n" << std::endl;
2759  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data written" << std::endl;
2760 
2761  double writeTime = MiscGetEllapsedSeconds(&timevalBegin);
2762  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2763  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedWriteContents()"
2764  << ": worldRank " << m_env.worldRank()
2765  << ", fullRank " << m_env.fullRank()
2766  << ", subEnvironment " << m_env.subId()
2767  << ", subRank " << m_env.subRank()
2768  << ", inter0Rank " << m_env.inter0Rank()
2769  << ", fileName = " << fileName
2770  << ", numParams = " << numParams
2771  << ", chainSize = " << chainSize
2772  << ", writeTime = " << writeTime << " seconds"
2773  << std::endl;
2774  }
2775 
2776  H5Dclose(dataset);
2777  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set closed" << std::endl;
2778  H5Sclose(dataspace);
2779  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space closed" << std::endl;
2780  H5Tclose(datatype);
2781  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type closed" << std::endl;
2782  }
2783  else {
2784  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
2785  }
2786  }
2787 #endif
2788  } // if (m_env.openUnifiedOutputFile())
2789  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 004 \n" << std::endl;
2790  } // if (m_env.inter0Rank() == (int) r)
2791  m_env.inter0Comm().Barrier();
2792  } // for r
2793 
2794  if (m_env.inter0Rank() == 0) {
2795  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2796  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2797  FilePtrSetStruct unifiedFilePtrSet;
2798  if (m_env.openUnifiedOutputFile(fileName,
2799  fileType,
2800  false, // Yes, 'writeOver = false' in order to close the array for matlab
2801  unifiedFilePtrSet)) {
2802 
2803  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2804  *unifiedFilePtrSet.ofsVar << "];\n";
2805  }
2806 
2807  m_env.closeFile(unifiedFilePtrSet,fileType);
2808  }
2809  }
2810  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2811  // Do nothing
2812  }
2813  else {
2814  queso_error_msg("invalid file type");
2815  }
2816  }
2817  } // if (m_env.inter0Rank() >= 0)
2818 
2819  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
2820  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedWriteContents()"
2821  << ", fileName = " << fileName
2822  << ", fileType = " << fileType
2823  << std::endl;
2824  }
2825  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2826 
2827  return;
2828 }
2829 
2830 template <class T>
2831 void
2832 ScalarSequence<T>::writeUnifiedMatlabHeader(std::ofstream & ofs,
2833  double sequenceSize) const
2834 {
2835  ofs << m_name << "_unified" << " = zeros(" << sequenceSize
2836  << "," << 1
2837  << ");"
2838  << std::endl;
2839  ofs << m_name << "_unified" << " = [";
2840 }
2841 
2842 template <class T>
2843 void
2844 ScalarSequence<T>::writeSubMatlabHeader(std::ofstream & ofs,
2845  double sequenceSize) const
2846 {
2847  ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << sequenceSize
2848  << "," << 1
2849  << ");"
2850  << std::endl;
2851  ofs << m_name << "_sub" << m_env.subIdString() << " = [";
2852 }
2853 
2854 template <class T>
2855 void
2856 ScalarSequence<T>::writeTxtHeader(std::ofstream & ofs,
2857  double sequenceSize) const
2858 {
2859  ofs << sequenceSize << " " << 1 << std::endl;
2860 }
2861 
2862 // --------------------------------------------------
2863 template <class T>
2864 void
2865 ScalarSequence<T>::unifiedReadContents(
2866  const std::string& fileName,
2867  const std::string& inputFileType,
2868  const unsigned int subReadSize)
2869 {
2870  queso_require_not_equal_to_msg(inputFileType, std::string(UQ_FILE_EXTENSION_FOR_TXT_FORMAT), std::string("reading txt files is not yet supported"));
2871  std::string fileType(inputFileType);
2872 #ifdef QUESO_HAS_HDF5
2873  // Do nothing
2874 #else
2875  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2876  if (m_env.subDisplayFile()) {
2877  *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedReadContents()"
2878  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2879  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2880  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2881  << "' instead..."
2882  << std::endl;
2883  }
2884  if (m_env.subRank() == 0) {
2885  std::cerr << "WARNING in ScalarSequence<T>::unifiedReadContents()"
2886  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2887  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2888  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2889  << "' instead..."
2890  << std::endl;
2891  }
2892  fileType = UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT;
2893  }
2894 #endif
2895 
2896  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2897  if (m_env.subDisplayFile()) {
2898  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedReadContents()"
2899  << ": worldRank " << m_env.worldRank()
2900  << ", fullRank " << m_env.fullRank()
2901  << ", subEnvironment " << m_env.subId()
2902  << ", subRank " << m_env.subRank()
2903  << ", inter0Rank " << m_env.inter0Rank()
2904  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
2905  << ", fileName = " << fileName
2906  << ", subReadSize = " << subReadSize
2907  //<< ", unifiedReadSize = " << unifiedReadSize
2908  << std::endl;
2909  }
2910 
2911  this->resizeSequence(subReadSize);
2912 
2913  if (m_env.inter0Rank() >= 0) {
2914  double unifiedReadSize = subReadSize*m_env.inter0Comm().NumProc();
2915 
2916  // In the logic below, the id of a line' begins with value 0 (zero)
2917  unsigned int idOfMyFirstLine = 1 + m_env.inter0Rank()*subReadSize;
2918  unsigned int idOfMyLastLine = (1 + m_env.inter0Rank())*subReadSize;
2919  unsigned int numParams = 1; // this->vectorSizeLocal();
2920 
2921  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // "m or hdf"
2922  if (m_env.inter0Rank() == (int) r) {
2923  // My turn
2924  FilePtrSetStruct unifiedFilePtrSet;
2925  if (m_env.openUnifiedInputFile(fileName,
2926  fileType,
2927  unifiedFilePtrSet)) {
2928  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2929  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2930  if (r == 0) {
2931  // Read number of chain positions in the file by taking care of the first line,
2932  // which resembles something like 'variable_name = zeros(n_positions,m_params);'
2933  std::string tmpString;
2934 
2935  // Read 'variable name' string
2936  *unifiedFilePtrSet.ifsVar >> tmpString;
2937  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2938 
2939  // Read '=' sign
2940  *unifiedFilePtrSet.ifsVar >> tmpString;
2941  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2942  queso_require_equal_to_msg(tmpString, std::string("="), std::string("string should be the '=' sign"));
2943 
2944  // Read 'zeros(n_positions,n_params)' string
2945  *unifiedFilePtrSet.ifsVar >> tmpString;
2946  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2947  unsigned int posInTmpString = 6;
2948 
2949  // Isolate 'n_positions' in a string
2950  //char nPositionsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
2951  std::string nPositionsString((size_t) (tmpString.size()-posInTmpString+1),' ');
2952  unsigned int posInPositionsString = 0;
2953  do {
2954  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
2955  nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
2956  } while (tmpString[posInTmpString] != ',');
2957  nPositionsString[posInPositionsString] = '\0';
2958 
2959  // Isolate 'n_params' in a string
2960  posInTmpString++; // Avoid reading ',' char
2961  //char nParamsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
2962  std::string nParamsString((size_t) (tmpString.size()-posInTmpString+1),' ');
2963  unsigned int posInParamsString = 0;
2964  do {
2965  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
2966  nParamsString[posInParamsString++] = tmpString[posInTmpString++];
2967  } while (tmpString[posInTmpString] != ')');
2968  nParamsString[posInParamsString] = '\0';
2969 
2970  // Convert 'n_positions' and 'n_params' strings to numbers
2971  unsigned int sizeOfChainInFile = (unsigned int) strtod(nPositionsString.c_str(),NULL);
2972  unsigned int numParamsInFile = (unsigned int) strtod(nParamsString.c_str(), NULL);
2973  if (m_env.subDisplayFile()) {
2974  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
2975  << ": worldRank " << m_env.worldRank()
2976  << ", fullRank " << m_env.fullRank()
2977  << ", sizeOfChainInFile = " << sizeOfChainInFile
2978  << ", numParamsInFile = " << numParamsInFile
2979  << std::endl;
2980  }
2981 
2982  // Check if [size of chain in file] >= [requested unified sequence size]
2983  queso_require_greater_equal_msg(sizeOfChainInFile, unifiedReadSize, "size of chain in file is not big enough");
2984 
2985  // Check if [num params in file] == [num params in current chain]
2986  queso_require_equal_to_msg(numParamsInFile, numParams, "number of parameters of chain in file is different than number of parameters in this chain object");
2987  } // if (r == 0)
2988 
2989  // Code common to any core in 'inter0Comm', including core of rank 0
2990  unsigned int maxCharsPerLine = 64*numParams; // Up to about 60 characters to represent each parameter value
2991 
2992  unsigned int lineId = 0;
2993  while (lineId < idOfMyFirstLine) {
2994  unifiedFilePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
2995  lineId++;
2996  };
2997 
2998  if (r == 0) {
2999  // Take care of initial part of the first data line,
3000  // which resembles something like 'variable_name = [value1 value2 ...'
3001  std::string tmpString;
3002 
3003  // Read 'variable name' string
3004  *unifiedFilePtrSet.ifsVar >> tmpString;
3005  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
3006 
3007  // Read '=' sign
3008  *unifiedFilePtrSet.ifsVar >> tmpString;
3009  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
3010  queso_require_equal_to_msg(tmpString, std::string("="), std::string("in core 0, string should be the '=' sign"));
3011 
3012  // Take into account the ' [' portion
3013  std::streampos tmpPos = unifiedFilePtrSet.ifsVar->tellg();
3014  unifiedFilePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
3015  }
3016 
3017  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3018  while (lineId <= idOfMyLastLine) {
3019  for (unsigned int i = 0; i < numParams; ++i) {
3020  *unifiedFilePtrSet.ifsVar >> tmpScalar;
3021  }
3022  m_seq[lineId - idOfMyFirstLine] = tmpScalar;
3023  lineId++;
3024  };
3025  }
3026 #ifdef QUESO_HAS_HDF5
3027  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
3028  if (r == 0) {
3029  hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var,
3030  "data",
3031  H5P_DEFAULT); // Dataset access property list
3032  hid_t datatype = H5Dget_type(dataset);
3033  H5T_class_t t_class = H5Tget_class(datatype);
3034  queso_require_equal_to_msg(t_class, H5T_FLOAT, "t_class is not H5T_DOUBLE");
3035  hid_t dataspace = H5Dget_space(dataset);
3036  int rank = H5Sget_simple_extent_ndims(dataspace);
3037  queso_require_equal_to_msg(rank, 2, "hdf rank is not 2");
3038  hsize_t dims_in[2];
3039  int status_n;
3040  status_n = H5Sget_simple_extent_dims(dataspace, dims_in, NULL);
3041  if (status_n) {}; // just to remove compiler warning
3042  //std::cout << "In ScalarSequence<T>::unifiedReadContents()"
3043  // << ": dims_in[0] = " << dims_in[0]
3044  // << ", dims_in[1] = " << dims_in[1]
3045  // << std::endl;
3046  queso_require_equal_to_msg(dims_in[0], numParams, "dims_in[0] is not equal to 'numParams'");
3047  queso_require_greater_equal_msg(dims_in[1], subReadSize, "dims_in[1] is smaller that requested 'subReadSize'");
3048 
3049  struct timeval timevalBegin;
3050  int iRC = UQ_OK_RC;
3051  iRC = gettimeofday(&timevalBegin,NULL);
3052  if (iRC) {}; // just to remove compiler warning
3053 
3054  unsigned int chainSizeIn = (unsigned int) dims_in[1];
3055  //double* dataIn[numParams]; // avoid compiler warning
3056  std::vector<double*> dataIn((size_t) numParams,NULL);
3057  dataIn[0] = (double*) malloc(numParams*chainSizeIn*sizeof(double));
3058  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
3059  dataIn[i] = dataIn[i-1] + chainSizeIn; // Yes, just 'chainSizeIn', not 'chainSizeIn*sizeof(double)'
3060  }
3061  //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, memory allocated" << std::endl;
3062  herr_t status;
3063  status = H5Dread(dataset,
3064  H5T_NATIVE_DOUBLE,
3065  H5S_ALL,
3066  dataspace,
3067  H5P_DEFAULT,
3068  dataIn[0]);
3069  if (status) {}; // just to remove compiler warning
3070  //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, data read" << std::endl;
3071  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3072  for (unsigned int j = 0; j < subReadSize; ++j) { // Yes, 'subReadSize', not 'chainSizeIn'
3073  for (unsigned int i = 0; i < numParams; ++i) {
3074  tmpScalar = dataIn[i][j];
3075  }
3076  m_seq[j] = tmpScalar;
3077  }
3078 
3079  double readTime = MiscGetEllapsedSeconds(&timevalBegin);
3080  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
3081  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
3082  << ": worldRank " << m_env.worldRank()
3083  << ", fullRank " << m_env.fullRank()
3084  << ", subEnvironment " << m_env.subId()
3085  << ", subRank " << m_env.subRank()
3086  << ", inter0Rank " << m_env.inter0Rank()
3087  << ", fileName = " << fileName
3088  << ", numParams = " << numParams
3089  << ", chainSizeIn = " << chainSizeIn
3090  << ", subReadSize = " << subReadSize
3091  << ", readTime = " << readTime << " seconds"
3092  << std::endl;
3093  }
3094 
3095  H5Sclose(dataspace);
3096  H5Tclose(datatype);
3097  H5Dclose(dataset);
3098  //free(dataIn[0]); // related to the changes above for compiler warning
3099  for (unsigned int tmpIndex = 0; tmpIndex < dataIn.size(); tmpIndex++) {
3100  free (dataIn[tmpIndex]);
3101  }
3102  }
3103  else {
3104  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
3105  }
3106  }
3107 #endif
3108  else {
3109  queso_error_msg("invalid file type");
3110  }
3111  m_env.closeFile(unifiedFilePtrSet,fileType);
3112  } // if (m_env.openUnifiedInputFile())
3113  } // if (m_env.inter0Rank() == (int) r)
3114  m_env.inter0Comm().Barrier();
3115  } // for r
3116  } // if (m_env.inter0Rank() >= 0)
3117  else {
3118  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3119  for (unsigned int i = 1; i < subReadSize; ++i) {
3120  m_seq[i] = tmpScalar;
3121  }
3122  }
3123 
3124  if (m_env.subDisplayFile()) {
3125  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedReadContents()"
3126  << ", fileName = " << fileName
3127  << std::endl;
3128  }
3129  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
3130 
3131  return;
3132 }
3133 // Private methods-----------------------------------
3134 template <class T>
3135 void
3136 ScalarSequence<T>::copy(const ScalarSequence<T>& src)
3137 {
3138  m_name = src.m_name;
3139  m_seq.clear();
3140  m_seq.resize(src.subSequenceSize(),0.);
3141  for (unsigned int i = 0; i < m_seq.size(); ++i) {
3142  m_seq[i] = src.m_seq[i];
3143  }
3144  deleteStoredScalars();
3145 
3146  return;
3147 }
3148 // --------------------------------------------------
3149 template <class T>
3150 void
3151 ScalarSequence<T>::extractScalarSeq(
3152  unsigned int initialPos,
3153  unsigned int spacing,
3154  unsigned int numPos,
3155  ScalarSequence<T>& scalarSeq) const
3156 {
3157  scalarSeq.resizeSequence(numPos);
3158  if (spacing == 1) {
3159  for (unsigned int j = 0; j < numPos; ++j) {
3160  //if ((initialPos+j*spacing) > m_seq.size()) {
3161  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3162  // << ": initialPos = " << initialPos
3163  // << ", spacing = " << spacing
3164  // << ", numPos = " << numPos
3165  // << ", j = " << j
3166  // << ", position got too large"
3167  // << std::endl;
3168  //}
3169  scalarSeq[j] = m_seq[initialPos+j ];
3170  }
3171  }
3172  else {
3173  for (unsigned int j = 0; j < numPos; ++j) {
3174  //if ((initialPos+j*spacing) > m_seq.size()) {
3175  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3176  // << ": initialPos = " << initialPos
3177  // << ", spacing = " << spacing
3178  // << ", numPos = " << numPos
3179  // << ", j = " << j
3180  // << ", position got too large"
3181  // << std::endl;
3182  //}
3183  scalarSeq[j] = m_seq[initialPos+j*spacing];
3184  }
3185  }
3186 
3187  return;
3188 }
3189 // --------------------------------------------------
3190 template <class T>
3191 void
3192 ScalarSequence<T>::extractRawData(
3193  unsigned int initialPos,
3194  unsigned int spacing,
3195  unsigned int numPos,
3196  std::vector<double>& rawDataVec) const
3197 {
3198  rawDataVec.resize(numPos);
3199  if (spacing == 1) {
3200  for (unsigned int j = 0; j < numPos; ++j) {
3201  rawDataVec[j] = m_seq[initialPos+j ];
3202  }
3203  }
3204  else {
3205  for (unsigned int j = 0; j < numPos; ++j) {
3206  rawDataVec[j] = m_seq[initialPos+j*spacing];
3207  }
3208  }
3209 
3210  return;
3211 }
3212 // --------------------------------------------------
3213 template <class T>
3214 std::vector<T>&
3215 ScalarSequence<T>::rawData()
3216 {
3217  return m_seq;
3218 }
3219 
3220 // --------------------------------------------------
3221 template <class T>
3222 void
3223 ScalarSequence<T>::subSort()
3224 {
3225  std::sort(m_seq.begin(), m_seq.end());
3226  return;
3227 }
3228 // --------------------------------------------------
3229 // Acknowledgement: the tree idea in this routine was taken from
3230 // 'http://penguin.ewu.edu/~trolfe/ParallelMerge/ParallelMerge.html', as of March 08, 2009
3231 template <class T>
3232 void
3233 ScalarSequence<T>::parallelMerge(
3234  std::vector<T>& sortedBuffer,
3235  const std::vector<T>& leafData,
3236  unsigned int currentTreeLevel) const
3237 {
3238  int parentNode = m_env.inter0Rank() & ~(1 << currentTreeLevel);
3239 
3240  if (m_env.inter0Rank() >= 0) { // KAUST
3241  if (currentTreeLevel == 0) {
3242  // Leaf node: sort own local data.
3243  unsigned int leafDataSize = leafData.size();
3244  sortedBuffer.resize(leafDataSize,0.);
3245  for (unsigned int i = 0; i < leafDataSize; ++i) {
3246  sortedBuffer[i] = leafData[i];
3247  }
3248  std::sort(sortedBuffer.begin(), sortedBuffer.end());
3249  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3250  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3251  << ": tree node " << m_env.inter0Rank()
3252  << ", leaf sortedBuffer[0] = " << sortedBuffer[0]
3253  << ", leaf sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3254  << std::endl;
3255  }
3256  }
3257  else {
3258  int nextTreeLevel = currentTreeLevel - 1;
3259  int rightChildNode = m_env.inter0Rank() | (1 << nextTreeLevel);
3260 
3261  if (rightChildNode >= m_env.inter0Comm().NumProc()) { // No right child. Move down one level.
3262  this->parallelMerge(sortedBuffer,
3263  leafData,
3264  nextTreeLevel);
3265  }
3266  else {
3267  unsigned int uintBuffer[1];
3268  uintBuffer[0] = nextTreeLevel;
3269  m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_INIT_MPI_MSG,
3270  "ScalarSequence<T>::parallelMerge()",
3271  "failed MPI.Send() for init");
3272 
3273  this->parallelMerge(sortedBuffer,
3274  leafData,
3275  nextTreeLevel);
3276 
3277  // Prepare variable 'leftSortedBuffer': just copy own current sorted data.
3278  unsigned int leftSize = sortedBuffer.size();
3279  std::vector<T> leftSortedBuffer(leftSize,0.);
3280  for (unsigned int i = 0; i < leftSize; ++i) {
3281  leftSortedBuffer[i] = sortedBuffer[i];
3282  }
3283 
3284  // Prepare variable 'rightSortedBuffer': receive data from right child node.
3285  RawType_MPI_Status status;
3286  m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_SIZE_MPI_MSG, &status,
3287  "ScalarSequence<T>::parallelMerge()",
3288  "failed MPI.Recv() for size");
3289  //if (status) {}; // just to remove compiler warning
3290 
3291  unsigned int rightSize = uintBuffer[0];
3292  std::vector<T> rightSortedBuffer(rightSize,0.);
3293  m_env.inter0Comm().Recv((void *) &rightSortedBuffer[0], (int) rightSize, RawValue_MPI_DOUBLE, rightChildNode, SCALAR_SEQUENCE_DATA_MPI_MSG, &status,
3294  "ScalarSequence<T>::parallelMerge()",
3295  "failed MPI.Recv() for data");
3296 
3297  // Merge the two results into 'sortedBuffer'.
3298  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3299  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3300  << ": tree node " << m_env.inter0Rank()
3301  << " is combining " << leftSortedBuffer.size()
3302  << " left doubles with " << rightSortedBuffer.size()
3303  << " right doubles"
3304  << std::endl;
3305  }
3306 
3307  sortedBuffer.clear();
3308  sortedBuffer.resize(leftSortedBuffer.size()+rightSortedBuffer.size(),0.);
3309  unsigned int i = 0;
3310  unsigned int j = 0;
3311  unsigned int k = 0;
3312  while ((i < leftSize ) &&
3313  (j < rightSize)) {
3314  if (leftSortedBuffer[i] > rightSortedBuffer[j]) sortedBuffer[k++] = rightSortedBuffer[j++];
3315  else sortedBuffer[k++] = leftSortedBuffer [i++];
3316  }
3317  while (i < leftSize ) sortedBuffer[k++] = leftSortedBuffer [i++];
3318  while (j < rightSize) sortedBuffer[k++] = rightSortedBuffer[j++];
3319  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3320  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3321  << ": tree node " << m_env.inter0Rank()
3322  << ", merged sortedBuffer[0] = " << sortedBuffer[0]
3323  << ", merged sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3324  << std::endl;
3325  }
3326  }
3327  }
3328 
3329  if (parentNode != m_env.inter0Rank()) {
3330  // Transmit data to parent node.
3331  unsigned int uintBuffer[1];
3332  uintBuffer[0] = sortedBuffer.size();
3333  m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, parentNode, SCALAR_SEQUENCE_SIZE_MPI_MSG,
3334  "ScalarSequence<T>::parallelMerge()",
3335  "failed MPI.Send() for size");
3336 
3337  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
3338  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3339  << ": tree node " << m_env.inter0Rank()
3340  << " is sending " << sortedBuffer.size()
3341  << " doubles to tree node " << parentNode
3342  << ", with sortedBuffer[0] = " << sortedBuffer[0]
3343  << " and sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3344  << std::endl;
3345  }
3346 
3347  m_env.inter0Comm().Send((void *) &sortedBuffer[0], (int) sortedBuffer.size(), RawValue_MPI_DOUBLE, parentNode, SCALAR_SEQUENCE_DATA_MPI_MSG,
3348  "ScalarSequence<T>::parallelMerge()",
3349  "failed MPI.Send() for data");
3350  }
3351  } // KAUST
3352 
3353  return;
3354 }
3355 
3356 // --------------------------------------------------
3357 // Methods conditionally available ------------------
3358 // --------------------------------------------------
3359 // --------------------------------------------------
3360 
3361 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
3362 template <class T>
3363 void
3364 ScalarSequence<T>::subUniformlySampledMdf(
3365  unsigned int numEvaluationPoints,
3366  T& minDomainValue,
3367  T& maxDomainValue,
3368  std::vector<T>& mdfValues) const
3369 {
3370  T tmpMinValue;
3371  T tmpMaxValue;
3372  std::vector<T> centers(numEvaluationPoints,0.);
3373  std::vector<unsigned int> bins (numEvaluationPoints,0);
3374 
3375  subMinMaxExtra(0, // initialPos
3376  this->subSequenceSize(),
3377  tmpMinValue,
3378  tmpMaxValue);
3379  subHistogram(0, // initialPos,
3380  tmpMinValue,
3381  tmpMaxValue,
3382  centers,
3383  bins);
3384 
3385  minDomainValue = centers[0];
3386  maxDomainValue = centers[centers.size()-1];
3387  T binSize = (maxDomainValue - minDomainValue)/((double)(centers.size() - 1));
3388 
3389  unsigned int totalCount = 0;
3390  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
3391  totalCount += bins[i];
3392  }
3393 
3394  mdfValues.clear();
3395  mdfValues.resize(numEvaluationPoints);
3396  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
3397  mdfValues[i] = ((T) bins[i])/((T) totalCount)/binSize;
3398  }
3399 
3400  return;
3401 }
3402 // --------------------------------------------------
3403 template <class T>
3404 T
3405 ScalarSequence<T>::subMeanCltStd(
3406  unsigned int initialPos,
3407  unsigned int numPos,
3408  const T& meanValue) const
3409 {
3410  if (this->subSequenceSize() == 0) return 0.;
3411 
3412  bool bRC = ((initialPos < this->subSequenceSize()) &&
3413  (0 < numPos ) &&
3414  ((initialPos+numPos) <= this->subSequenceSize()));
3415  queso_require_msg(bRC, "invalid input data");
3416 
3417  unsigned int finalPosPlus1 = initialPos + numPos;
3418  T diff;
3419  T stdValue = 0.;
3420  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
3421  diff = m_seq[j] - meanValue;
3422  stdValue += diff*diff;
3423  }
3424 
3425  stdValue /= (((T) numPos) - 1.);
3426  stdValue /= (((T) numPos) - 1.);
3427  stdValue = sqrt(stdValue);
3428 
3429  return stdValue;
3430 }
3431 // --------------------------------------------------
3432 template <class T>
3433 T
3434 ScalarSequence<T>::unifiedMeanCltStd(
3435  bool useOnlyInter0Comm,
3436  unsigned int initialPos,
3437  unsigned int numPos,
3438  const T& unifiedMeanValue) const
3439 {
3440  if (m_env.numSubEnvironments() == 1) {
3441  return this->subMeanCltStd(initialPos,
3442  numPos,
3443  unifiedMeanValue);
3444  }
3445 
3446  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
3447 
3448  T unifiedStdValue = 0.;
3449  if (useOnlyInter0Comm) {
3450  if (m_env.inter0Rank() >= 0) {
3451  bool bRC = ((initialPos < this->subSequenceSize()) &&
3452  (0 < numPos ) &&
3453  ((initialPos+numPos) <= this->subSequenceSize()));
3454  queso_require_msg(bRC, "invalid input data");
3455 
3456  unsigned int finalPosPlus1 = initialPos + numPos;
3457  T diff;
3458  T localStdValue = 0.;
3459  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
3460  diff = m_seq[j] - unifiedMeanValue;
3461  localStdValue += diff*diff;
3462  }
3463 
3464  unsigned int unifiedNumPos = 0;
3465  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
3466  "ScalarSequence<T>::unifiedMeanCltStd()",
3467  "failed MPI.Allreduce() for numPos");
3468 
3469  m_env.inter0Comm().template Allreduce<double>(&localStdValue, &unifiedStdValue, (int) 1, RawValue_MPI_SUM,
3470  "ScalarSequence<T>::unifiedMeanCltStd()",
3471  "failed MPI.Allreduce() for stdValue");
3472 
3473  unifiedStdValue /= (((T) unifiedNumPos) - 1.);
3474  unifiedStdValue /= (((T) unifiedNumPos) - 1.);
3475  unifiedStdValue /= sqrt(unifiedStdValue);
3476  }
3477  else {
3478  // Node not in the 'inter0' communicator
3479  this->subMeanCltStd(initialPos,
3480  numPos,
3481  unifiedMeanValue);
3482  }
3483  }
3484  else {
3485  queso_error_msg("parallel vectors not supported yet");
3486  }
3487  //m_env.fullComm().Barrier();
3488 
3489  return unifiedStdValue;
3490 }
3491 //---------------------------------------------------
3492 template <class T>
3493 T
3494 ScalarSequence<T>::bmm(
3495  unsigned int initialPos,
3496  unsigned int batchLength) const
3497 {
3498  bool bRC = ((initialPos < this->subSequenceSize() ) &&
3499  (batchLength < (this->subSequenceSize()-initialPos)));
3500  queso_require_msg(bRC, "invalid input data");
3501 
3502  unsigned int numberOfBatches = (this->subSequenceSize() - initialPos)/batchLength;
3503  ScalarSequence<T> batchMeans(m_env,numberOfBatches,"");
3504 
3505  for (unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
3506  batchMeans[batchId] = this->subMeanExtra(initialPos + batchId*batchLength,
3507  batchLength);
3508  }
3509 
3510  T meanOfBatchMeans = batchMeans.subMeanExtra(0,
3511  batchMeans.subSequenceSize());
3512 
3513  //T covLag0OfBatchMeans = batchMeans.autoCovariance(0,
3514  // batchMeans.subSequenceSize(),
3515  // meanOfBatchMeans,
3516  // 0); // lag
3517 
3518  T bmmValue = batchMeans.subSampleVarianceExtra(0,
3519  batchMeans.subSequenceSize(),
3520  meanOfBatchMeans);
3521 
3522  bmmValue /= (T) batchMeans.subSequenceSize(); // CHECK
3523 //bmmValue *= (T) (this->subSequenceSize() - initialPos); // CHECK
3524 
3525  return bmmValue;
3526 }
3527 //---------------------------------------------------
3528 template <class T>
3529 void
3530 ScalarSequence<T>::psd(
3531  unsigned int initialPos,
3532  unsigned int numBlocks,
3533  double hopSizeRatio,
3534  std::vector<double>& psdResult) const
3535 {
3536  bool bRC = ((initialPos < this->subSequenceSize() ) &&
3537  (hopSizeRatio != 0. ) &&
3538  (numBlocks < (this->subSequenceSize() - initialPos)) &&
3539  (fabs(hopSizeRatio) < (double) (this->subSequenceSize() - initialPos)));
3540  queso_require_msg(bRC, "invalid input data");
3541 
3542  unsigned int dataSize = this->subSequenceSize() - initialPos;
3543 
3544  T meanValue = this->subMeanExtra(initialPos,
3545  dataSize);
3546 
3547  // Determine hopSize and blockSize
3548  unsigned int hopSize = 0;
3549  unsigned int blockSize = 0;
3550  if (hopSizeRatio <= -1.) {
3551  double overlapSize = -hopSizeRatio;
3552  double tmp = ((double) (dataSize + (numBlocks - 1)*overlapSize))/((double) numBlocks);
3553  blockSize = (unsigned int) tmp;
3554  }
3555  else if (hopSizeRatio < 0.) {
3556  double tmp = ((double) dataSize)/(( ((double) numBlocks) - 1. )*(-hopSizeRatio) + 1.);
3557  blockSize = (unsigned int) tmp;
3558  hopSize = (unsigned int) ( ((double) blockSize) * (-hopSizeRatio) );
3559  }
3560  else if (hopSizeRatio == 0.) {
3561  queso_error_msg("hopSizeRatio == 0");
3562  }
3563  else if (hopSizeRatio < 1.) {
3564  double tmp = ((double) dataSize)/(( ((double) numBlocks) - 1. )*hopSizeRatio + 1.);
3565  blockSize = (unsigned int) tmp;
3566  hopSize = (unsigned int) ( ((double) blockSize) * hopSizeRatio );
3567  }
3568  else {
3569  hopSize = (unsigned int) hopSizeRatio;
3570  blockSize = dataSize - (numBlocks - 1)*hopSize;
3571  }
3572 
3573  int numberOfDiscardedDataElements = dataSize - (numBlocks-1)*hopSize - blockSize;
3574  //if (m_env.subDisplayFile()) {
3575  // *m_env.subDisplayFile() << "initialPos = " << initialPos
3576  // << ", N = " << dataSize
3577  // << ", #Blocks = " << numBlocks
3578  // << ", R (hop size) = " << hopSize
3579  // << ", B (block size) = " << blockSize
3580  // << ", overlap = " << blockSize - hopSize
3581  // << ", [(#Blocks - 1) * R + B] = " << (numBlocks-1)*hopSize + blockSize
3582  // << ", numberOfDiscardedDataElements = " << numberOfDiscardedDataElements
3583  // << std::endl;
3584  //}
3585  queso_require_greater_equal_msg(numberOfDiscardedDataElements, 0., "eventual extra space for last block should not be negative");
3586 
3587  double tmp = log((double) blockSize)/log(2.);
3588  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
3589  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
3590  unsigned int fftSize = (unsigned int) std::pow(2.,tmp);
3591  //if (m_env.subDisplayFile()) {
3592  // *m_env.subDisplayFile() << "fractionalPart = " << fractionalPart
3593  // << ", B = " << blockSize
3594  // << ", fftSize = " << fftSize
3595  // << std::endl;
3596  //}
3597 
3598  double modificationScale = 0.;
3599  for (unsigned int j = 0; j < blockSize; ++j) {
3600  double tmpValue = MiscHammingWindow(blockSize-1,j);
3601  modificationScale += tmpValue*tmpValue;
3602  }
3603  modificationScale = 1./modificationScale;
3604 
3605  std::vector<double> blockData(blockSize,0.);
3606  Fft<T> fftObj(m_env);
3607  std::vector<std::complex<double> > fftResult(0,std::complex<double>(0.,0.));
3608 
3609 #if 0
3610  unsigned int halfFFTSize = fftSize/2;
3611  psdResult.clear();
3612  psdResult.resize(1+halfFFTSize,0.);
3613  double factor = 1./M_PI/((double) numBlocks); // /((double) blockSize);
3614 #else
3615  psdResult.clear();
3616  psdResult.resize(fftSize,0.);
3617  double factor = 1./2./M_PI/((double) numBlocks); // /((double) blockSize);
3618 #endif
3619 
3620  for (unsigned int blockId = 0; blockId < numBlocks; blockId++) {
3621  // Fill block using Hamming window
3622  unsigned int initialDataPos = initialPos + blockId*hopSize;
3623  for (unsigned int j = 0; j < blockSize; ++j) {
3624  unsigned int dataPos = initialDataPos + j;
3625  queso_require_less_msg(dataPos, dataSize, "too large position to be accessed in data");
3626  blockData[j] = MiscHammingWindow(blockSize-1,j) * ( m_seq[dataPos] - meanValue ); // IMPORTANT
3627  }
3628 
3629  fftObj.forward(blockData,fftSize,fftResult);
3630 
3631  //if (m_env.subDisplayFile()) {
3632  // *m_env.subDisplayFile() << "blockData.size() = " << blockData.size()
3633  // << ", fftSize = " << fftSize
3634  // << ", fftResult.size() = " << fftResult.size()
3635  // << ", psdResult.size() = " << psdResult.size()
3636  // << std::endl;
3637  //}
3638  // Normalized spectral density: power per radians per sample
3639  for (unsigned int j = 0; j < psdResult.size(); ++j) {
3640  psdResult[j] += norm(fftResult[j])*modificationScale;
3641  }
3642  }
3643 
3644  for (unsigned int j = 0; j < psdResult.size(); ++j) {
3645  psdResult[j] *= factor;
3646  }
3647 
3648  return;
3649 }
3650 //---------------------------------------------------
3651 template <class T>
3652 T
3653 ScalarSequence<T>::geweke(
3654  unsigned int initialPos,
3655  double ratioNa,
3656  double ratioNb) const
3657 {
3658  double doubleFullDataSize = (double) (this->subSequenceSize() - initialPos);
3659  ScalarSequence<T> tmpSeq(m_env,0,"");
3660  std::vector<double> psdResult(0,0.);
3661 
3662  unsigned int dataSizeA = (unsigned int) (doubleFullDataSize * ratioNa);
3663  double doubleDataSizeA = (double) dataSizeA;
3664  unsigned int initialPosA = initialPos;
3665  this->extractScalarSeq(initialPosA,
3666  1,
3667  dataSizeA,
3668  tmpSeq);
3669  double meanA = tmpSeq.subMeanExtra(0,
3670  dataSizeA);
3671  tmpSeq.psd(0,
3672  (unsigned int) std::sqrt((double) dataSizeA), // numBlocks
3673  .5, // hopSizeRatio
3674  psdResult);
3675  double psdA = psdResult[0];
3676  double varOfMeanA = 2.*M_PI*psdA/doubleDataSizeA;
3677 
3678  unsigned int dataSizeB = (unsigned int) (doubleFullDataSize * ratioNb);
3679  double doubleDataSizeB = (double) dataSizeB;
3680  unsigned int initialPosB = this->subSequenceSize() - dataSizeB;
3681  this->extractScalarSeq(initialPosB,
3682  1,
3683  dataSizeB,
3684  tmpSeq);
3685  double meanB = tmpSeq.subMeanExtra(0,
3686  dataSizeB);
3687  tmpSeq.psd(0,
3688  (unsigned int) std::sqrt((double) dataSizeB), // numBlocks
3689  .5, // hopSizeRatio
3690  psdResult);
3691  double psdB = psdResult[0];
3692  double varOfMeanB = 2.*M_PI*psdB/doubleDataSizeB;
3693 
3694  if (m_env.subDisplayFile()) {
3695  *m_env.subDisplayFile() << "In ScalarSequence<T>::geweke()"
3696  << ", before computation of gewCoef"
3697  << ":\n"
3698  << ", dataSizeA = " << dataSizeA
3699  << ", numBlocksA = " << (unsigned int) std::sqrt((double) dataSizeA)
3700  << ", meanA = " << meanA
3701  << ", psdA = " << psdA
3702  << ", varOfMeanA = " << varOfMeanA
3703  << "\n"
3704  << ", dataSizeB = " << dataSizeB
3705  << ", numBlocksB = " << (unsigned int) std::sqrt((double) dataSizeB)
3706  << ", meanB = " << meanB
3707  << ", psdB = " << psdB
3708  << ", varOfMeanB = " << varOfMeanB
3709  << std::endl;
3710  }
3711  double gewCoef = (meanA - meanB)/std::sqrt(varOfMeanA + varOfMeanB);
3712 
3713  return gewCoef;
3714 }
3715 //---------------------------------------------------
3716 template <class T>
3717 T
3718 ScalarSequence<T>::meanStacc(
3719  unsigned int initialPos) const
3720 {
3721  queso_not_implemented(); // ERNESTO
3722 
3723  double value = 0.;
3724 
3725  return value;
3726 }
3727 //---------------------------------------------------
3728 template <class T>
3729 void
3730 ScalarSequence<T>::subCdfPercentageRange(
3731  unsigned int initialPos,
3732  unsigned int numPos,
3733  double range, // \in [0,1]
3734  T& lowerValue,
3735  T& upperValue) const
3736 {
3737  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
3738 
3739  queso_require_msg(!((range < 0) || (range > 1.)), "invalid 'range' value");
3740 
3741  ScalarSequence<T> sortedSequence(m_env,0,"");;
3742  sortedSequence.resizeSequence(numPos);
3743  this->extractScalarSeq(initialPos,
3744  1,
3745  numPos,
3746  sortedSequence);
3747  sortedSequence.subSort();
3748 
3749  unsigned int lowerId = (unsigned int) round( 0.5*(1.-range)*((double) numPos) );
3750  lowerValue = sortedSequence[lowerId];
3751 
3752  unsigned int upperId = (unsigned int) round( 0.5*(1.+range)*((double) numPos) );
3753  if (upperId == numPos) {
3754  upperId = upperId-1;
3755  }
3756  queso_require_less_msg(upperId, numPos, "'upperId' got too big");
3757  upperValue = sortedSequence[upperId];
3758 
3759  return;
3760 }
3761 //---------------------------------------------------
3762 template <class T>
3763 void
3764 ScalarSequence<T>::unifiedCdfPercentageRange(
3765  bool useOnlyInter0Comm,
3766  unsigned int initialPos,
3767  unsigned int numPos,
3768  double range, // \in [0,1]
3769  T& unifiedLowerValue,
3770  T& unifiedUpperValue) const
3771 {
3772  if (m_env.numSubEnvironments() == 1) {
3773  return this->subCdfPercentageRange(initialPos,
3774  numPos,
3775  range,
3776  unifiedLowerValue,
3777  unifiedUpperValue);
3778  }
3779 
3780  // As of 07/Feb/2011, this routine does *not* require sub sequences to have equal size. Good.
3781 
3782  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
3783 
3784  queso_require_msg(!((range < 0) || (range > 1.)), "invalid 'range' value");
3785 
3786  if (useOnlyInter0Comm) {
3787  if (m_env.inter0Rank() >= 0) {
3788  queso_error_msg("code not implemented yet");
3789  }
3790  else {
3791  // Node not in the 'inter0' communicator
3792  this->subCdfPercentageRange(initialPos,
3793  numPos,
3794  unifiedLowerValue,
3795  unifiedUpperValue);
3796  }
3797  }
3798  else {
3799  queso_error_msg("parallel vectors not supported yet");
3800  }
3801  //m_env.fullComm().Barrier();
3802 
3803  return;
3804 }
3805 //---------------------------------------------------
3806 template <class T>
3807 void
3808 ScalarSequence<T>::subCdfStacc(
3809  unsigned int initialPos,
3810  std::vector<double>& cdfStaccValues,
3811  std::vector<double>& cdfStaccValuesUp,
3812  std::vector<double>& cdfStaccValuesLow,
3813  const ScalarSequence<T>& sortedDataValues) const
3814 {
3815  queso_require_msg(!(false), "not implemented yet"); // Joseph
3816  bool bRC = (initialPos < this->subSequenceSize() );
3817  queso_require_msg(bRC, "invalid input data");
3818 
3819  unsigned int numPoints = subSequenceSize()-initialPos;
3820  double auxNumPoints = numPoints;
3821  double maxLamb = 0.;
3822  std::vector<double> ro (numPoints,0.);
3823  std::vector<double> Isam_mat(numPoints,0.);
3824 
3825  for (unsigned int pointId = 0; pointId < numPoints; pointId++) {
3826  double p = ( ((double) pointId) + 1.0 )/auxNumPoints;
3827  double ro0 = p*(1.0-p);
3828  cdfStaccValues[pointId] = p;
3829 
3830  //std::cout << "x-data" << data[pointId]
3831  // << std::endl;
3832 
3833  for (unsigned int k = 0; k < numPoints; k++) {
3834  if (m_seq[k] <= sortedDataValues[pointId]) {
3835  Isam_mat[k] = 1;
3836  }
3837  else {
3838  Isam_mat[k] = 0;
3839  }
3840  }
3841 
3842  for (unsigned int tau = 0; tau < (numPoints-1); tau++) {
3843  ro[tau] = 0.;
3844  for (unsigned int kk = 0; kk < numPoints-(tau+1); kk++) {
3845  ro[tau] += (Isam_mat[kk+tau+1]-p)*(Isam_mat[kk]-p);
3846  }
3847  ro[tau] *= 1.0/auxNumPoints;
3848  }
3849  double lamb = 0.;
3850  for (unsigned int tau = 0; tau < (numPoints-1); tau++) {
3851  double auxTau = tau;
3852  lamb += (1.-(auxTau+1.)/auxNumPoints)*ro[tau]/ro0;
3853  if (lamb > maxLamb) maxLamb = lamb;
3854  }
3855  lamb = 2.*maxLamb;
3856 
3857  //double ll=gsl_cdf_gaussian_Pinv(-0.05);
3858  cdfStaccValuesUp [pointId] = cdfStaccValues[pointId]+1.96*sqrt(ro0/auxNumPoints*(1.+lamb));
3859  cdfStaccValuesLow[pointId] = cdfStaccValues[pointId]-1.96*sqrt(ro0/auxNumPoints*(1.+lamb));
3860  if (cdfStaccValuesLow[pointId] < 0.) cdfStaccValuesLow[pointId] = 0.;
3861  if (cdfStaccValuesUp [pointId] > 1.) cdfStaccValuesUp [pointId] = 1.;
3862  //if (pointId==1 | pointId==100 | pointId==900) {
3863  // std::cout<<lamb<<ll<<std::endl;
3864  // std::cout<<lamb<<std::endl;
3865  //}
3866  }
3867 
3868 #if 0
3869  unsigned int dataSize = this->subSequenceSize() - initialPos;
3870  unsigned int numEvals = evaluationPositions.size();
3871 
3872  for (unsigned int j = 0; j < numEvals; ++j) {
3873  double x = evaluationPositions[j];
3874  double value = 0.;
3875  for (unsigned int k = 0; k < dataSize; ++k) {
3876  double xk = m_seq[initialPos+k];
3877  value += 0.;//MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
3878  }
3879  cdfStaccValues[j] = value/(double) dataSize;
3880  }
3881 #endif
3882 
3883  return;
3884 }
3885 //---------------------------------------------------
3886 template <class T>
3887 void
3888 ScalarSequence<T>::subCdfStacc(
3889  unsigned int initialPos,
3890  const std::vector<T>& evaluationPositions,
3891  std::vector<double>& cdfStaccValues) const
3892 {
3893  queso_not_implemented();
3894 
3895  bool bRC = ((initialPos < this->subSequenceSize() ) &&
3896  (0 < evaluationPositions.size()) &&
3897  (evaluationPositions.size() == cdfStaccValues.size() ));
3898  queso_require_msg(bRC, "invalid input data");
3899 
3900  // For Joseph:
3901  // Maybe sort first
3902  // For each of the evaluation positions:
3903  // --> 1) form temporary scalar seq that contains only 0 and 1
3904  // --> 2) call "temporary scalar seq"->meanStacc()
3905 
3906 #if 0
3907  unsigned int dataSize = this->subSequenceSize() - initialPos;
3908  unsigned int numEvals = evaluationPositions.size();
3909 
3910  for (unsigned int j = 0; j < numEvals; ++j) {
3911  //double x = evaluationPositions[j];
3912  double value = 0.;
3913  for (unsigned int k = 0; k < dataSize; ++k) {
3914  //double xk = m_seq[initialPos+k];
3915  value += 0.;//MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
3916  }
3917  cdfStaccValues[j] = value/(double) dataSize;
3918  }
3919 #endif
3920 
3921  return;
3922 }
3923 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
3924 
3925 // --------------------------------------------------
3926 // Additional methods -------------------------------
3927 // (outside class declaration) ----------------------
3928 //---------------------------------------------------
3929 template <class T>
3930 void
3931 ComputeSubGaussian2dKde(const ScalarSequence<T>& scalarSeq1,
3932  const ScalarSequence<T>& scalarSeq2,
3933  unsigned int initialPos,
3934  double scaleValue1,
3935  double scaleValue2,
3936  const std::vector<T>& evaluationPositions1,
3937  const std::vector<T>& evaluationPositions2,
3938  std::vector<double>& densityValues)
3939 {
3940  queso_require_equal_to_msg(initialPos, 0, "not implemented yet for initialPos != 0");
3941 
3942  double covValue = 0.;
3943  double corrValue = 0.;
3944  ComputeCovCorrBetweenScalarSequences(scalarSeq1,
3945  scalarSeq2,
3946  scalarSeq1.subSequenceSize(),
3947  covValue,
3948  corrValue);
3949 
3950  bool bRC = ((initialPos < scalarSeq1.subSequenceSize()) &&
3951  (scalarSeq1.subSequenceSize() == scalarSeq2.subSequenceSize()) &&
3952  (0 < evaluationPositions1.size() ) &&
3953  (evaluationPositions1.size() == evaluationPositions2.size() ) &&
3954  (evaluationPositions1.size() == densityValues.size() ));
3955  queso_require_msg(bRC, "invalid input data");
3956 
3957  unsigned int dataSize = scalarSeq1.subSequenceSize() - initialPos;
3958  unsigned int numEvals = evaluationPositions1.size();
3959 
3960  double scale1Inv = 1./scaleValue1;
3961  double scale2Inv = 1./scaleValue2;
3962  //corrValue = 0.;
3963  double r = 1. - corrValue*corrValue;
3964  if (r <= 0.) { // prudencio 2010-07-23
3965  std::cerr << "In ComputeSubGaussian2dKde()"
3966  << ": WARNING, r = " << r
3967  << std::endl;
3968  }
3969 
3970  // scalarSeq1.env().worldRank(),
3971  // "ComputeSubGaussian2dKde()",
3972  // "negative r");
3973  for (unsigned int j = 0; j < numEvals; ++j) {
3974  double x1 = evaluationPositions1[j];
3975  double x2 = evaluationPositions2[j];
3976  double value = 0.;
3977  for (unsigned int k = 0; k < dataSize; ++k) {
3978  double d1k = scale1Inv*(x1 - scalarSeq1[initialPos+k]);
3979  double d2k = scale2Inv*(x2 - scalarSeq2[initialPos+k]);
3980  value += exp(-.5*( d1k*d1k + 2*corrValue*d1k*d2k + d2k*d2k )/r);
3981  }
3982  densityValues[j] = scale1Inv * scale2Inv * (value/(double) dataSize) / 2. / M_PI / sqrt(r);
3983  }
3984 
3985  return;
3986 }
3987 //---------------------------------------------------
3988 template <class T>
3989 void
3990 ComputeUnifiedGaussian2dKde(bool useOnlyInter0Comm, // INPUT
3991  const ScalarSequence<T>& scalarSeq1, // INPUT
3992  const ScalarSequence<T>& scalarSeq2, // INPUT
3993  unsigned int initialPos, // INPUT
3994  double unifiedScaleValue1, // INPUT
3995  double unifiedScaleValue2, // INPUT
3996  const std::vector<T>& unifiedEvaluationPositions1, // INPUT
3997  const std::vector<T>& unifiedEvaluationPositions2, // INPUT
3998  std::vector<double>& unifiedDensityValues) // OUTPUT
3999 {
4000  if (scalarSeq1.env().numSubEnvironments() == 1) {
4001  return ComputeSubGaussian2dKde(scalarSeq1,
4002  scalarSeq2,
4003  initialPos,
4004  unifiedScaleValue1,
4005  unifiedScaleValue2,
4006  unifiedEvaluationPositions1,
4007  unifiedEvaluationPositions2,
4008  unifiedDensityValues);
4009  }
4010 
4011  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
4012 
4013  if (useOnlyInter0Comm) {
4014  if (scalarSeq1.env().inter0Rank() >= 0) {
4015  queso_error_msg("inter0 case not supported yet");
4016  }
4017  else {
4018  // Node not in the 'inter0' communicator
4019  ComputeSubGaussian2dKde(scalarSeq1,
4020  scalarSeq2,
4021  initialPos,
4022  unifiedScaleValue1,
4023  unifiedScaleValue2,
4024  unifiedEvaluationPositions1,
4025  unifiedEvaluationPositions2,
4026  unifiedDensityValues);
4027  }
4028  }
4029  else {
4030  queso_error_msg("parallel vectors not supported yet");
4031  }
4032 
4033  //scalarSeq1.env().fullComm().Barrier();
4034 
4035  return;
4036 }
4037 // --------------------------------------------------
4038 template <class T>
4039 void
4040 ComputeCovCorrBetweenScalarSequences(
4041  const ScalarSequence<T>& subPSeq,
4042  const ScalarSequence<T>& subQSeq,
4043  unsigned int subNumSamples,
4044  T& covValue,
4045  T& corrValue)
4046 {
4047  // Check input data consistency
4048  const BaseEnvironment& env = subPSeq.env();
4049 
4050  queso_require_msg(!((subNumSamples > subPSeq.subSequenceSize()) || (subNumSamples > subQSeq.subSequenceSize())), "subNumSamples is too large");
4051 
4052  // For both P and Q vector sequences: fill them
4053  T tmpP = 0.;
4054  T tmpQ = 0.;
4055 
4056  // For both P and Q vector sequences: compute the unified mean
4057  T unifiedMeanP = subPSeq.unifiedMeanExtra(true,0,subNumSamples);
4058  T unifiedMeanQ = subQSeq.unifiedMeanExtra(true,0,subNumSamples);
4059 
4060  // Compute "sub" covariance matrix
4061  covValue = 0.;
4062  for (unsigned k = 0; k < subNumSamples; ++k) {
4063  // For both P and Q vector sequences: get the difference (wrt the unified mean) in them
4064  tmpP = subPSeq[k] - unifiedMeanP;
4065  tmpQ = subQSeq[k] - unifiedMeanQ;
4066  covValue += tmpP*tmpQ;
4067  }
4068 
4069  // For both P and Q vector sequences: compute the unified variance
4070  T unifiedSampleVarianceP = subPSeq.unifiedSampleVarianceExtra(true,
4071  0,
4072  subNumSamples,
4073  unifiedMeanP);
4074 
4075  T unifiedSampleVarianceQ = subQSeq.unifiedSampleVarianceExtra(true,
4076  0,
4077  subNumSamples,
4078  unifiedMeanQ);
4079 
4080  // Compute unified covariance
4081  if (env.inter0Rank() >= 0) {
4082  unsigned int unifiedNumSamples = 0;
4083  env.inter0Comm().template Allreduce<unsigned int>(&subNumSamples, &unifiedNumSamples, (int) 1, RawValue_MPI_SUM,
4084  "ComputeCovCorrBetweenScalarSequences()",
4085  "failed MPI.Allreduce() for subNumSamples");
4086 
4087  double aux = 0.;
4088  env.inter0Comm().template Allreduce<double>(&covValue, &aux, (int) 1, RawValue_MPI_SUM,
4089  "ComputeCovCorrBetweenScalarSequences()",
4090  "failed MPI.Allreduce() for a matrix position");
4091  covValue = aux/((double) (unifiedNumSamples-1)); // Yes, '-1' in order to compensate for the 'N-1' denominator factor in the calculations of sample variances above (whose square roots will be used below)
4092 
4093  corrValue = covValue/std::sqrt(unifiedSampleVarianceP)/std::sqrt(unifiedSampleVarianceQ);
4094 
4095  if ((corrValue < -1.) || (corrValue > 1.)) { // prudencio 2010-07-23
4096  std::cerr << "In ComputeCovCorrBetweenScalarSequences()"
4097  << ": computed correlation is out of range, corrValue = " << corrValue
4098  << std::endl;
4099  }
4100 
4101  // env.worldRank(),
4102  // "ComputeCovCorrBetweenScalarSequences()",
4103  // "computed correlation is out of range");
4104  }
4105  else {
4106  // Node not in the 'inter0' communicator: do nothing extra
4107  }
4108 
4109  return;
4110 }
4111 
4112 } // End namespace QUESO
4113 
4114 template class QUESO::ScalarSequence<double>;
QUESO::ScalarSequence::unifiedScaleForKde
T unifiedScaleForKde(bool useOnlyInter0Comm, unsigned int initialPos, const T &unifiedIqrValue, unsigned int kdeDimension) const
Selects the scales (bandwidth) for the kernel density estimation, considering the unified sequence...
Definition: ScalarSequence.C:2207
QUESO::Fft::forward
void forward(const std::vector< T > &data, unsigned int fftSize, std::vector< std::complex< double > > &result)
Calculates the forward Fourier transform (for real data. TODO: complex data).
QUESO::ScalarSequence::subMaxPlain
const T & subMaxPlain() const
Finds the maximum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:369
QUESO::ScalarSequence::subInterQuantileRange
T subInterQuantileRange(unsigned int initialPos) const
Returns the interquartile range of the values in the sub-sequence.
Definition: ScalarSequence.C:1963
QUESO::ScalarSequence::subMeanPlain
const T & subMeanPlain() const
Finds the mean value of the sub-sequence of scalars.
Definition: ScalarSequence.C:395
QUESO::ScalarSequence::ScalarSequence
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
Definition: ScalarSequence.C:32
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
QUESO::ScalarSequence::unifiedMinPlain
const T & unifiedMinPlain(bool useOnlyInter0Comm) const
Finds the minimum value of the unified sequence of scalars.
Definition: ScalarSequence.C:356
QUESO::ScalarSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
Definition: ScalarSequence.C:143
QUESO::ScalarSequence::~ScalarSequence
~ScalarSequence()
Destructor.
Definition: ScalarSequence.C:54
QUESO::ScalarSequence::subMinPlain
const T & subMinPlain() const
Finds the minimum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:343
QUESO::ScalarSequence::operator=
ScalarSequence< T > & operator=(const ScalarSequence< T > &rhs)
Assignment operator; it copies rhs to this.
Definition: ScalarSequence.C:61
QUESO::ScalarSequence::unifiedMeanExtra
T unifiedMeanExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos) const
Finds the mean value of the unified sequence of numPos positions starting at position initialPos...
Definition: ScalarSequence.C:844
QUESO::BaseEnvironment
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
QUESO::ScalarSequence::erasePositions
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values of the sequence, starting at position initialPos.
Definition: ScalarSequence.C:206
QUESO::ScalarSequence::subMeanExtra
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
Definition: ScalarSequence.C:806
QUESO::ScalarSequence::unifiedSampleVarianceExtra
T unifiedSampleVarianceExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos, const T &unifiedMeanValue) const
Finds the sample variance of the unified sequence, considering numPos positions starting at position ...
Definition: ScalarSequence.C:1039
QUESO::ScalarSequence::unifiedGaussian1dKde
void unifiedGaussian1dKde(bool useOnlyInter0Comm, unsigned int initialPos, double unifiedScaleValue, const std::vector< T > &unifiedEvaluationPositions, std::vector< double > &unifiedDensityValues) const
Gaussian kernel for the KDE estimate of the unified sequence.
Definition: ScalarSequence.C:2308
QUESO::ScalarSequence::brooksGelmanConvMeasure
T brooksGelmanConvMeasure(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int spacing) const
Estimates convergence rate using Brooks &amp; Gelman method.
Definition: ScalarSequence.C:2427
QUESO::ScalarSequence::writeUnifiedMatlabHeader
void writeUnifiedMatlabHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write header info for matlab files from all chains.
Definition: ScalarSequence.C:2832
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:92
QUESO::ScalarSequence::unifiedMeanPlain
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Definition: ScalarSequence.C:407
QUESO::MiscHammingWindow
double MiscHammingWindow(unsigned int N, unsigned int j)
Definition: Miscellaneous.C:298
QUESO::ScalarSequence::unifiedCdfPercentageRange
void unifiedCdfPercentageRange(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, double range, T &lowerValue, T &upperValue) const
Definition: ScalarSequence.C:3764
QUESO::queso_require_not_equal_to_msg
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
Definition: MonteCarloSGOptions.C:417
QUESO::ScalarSequence::unifiedSampleStd
T unifiedSampleStd(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos, const T &unifiedMeanValue) const
Finds the sample standard deviation of the unified sequence, considering localnumPos positions starti...
Definition: ScalarSequence.C:1126
QUESO::BaseEnvironment::inter0Comm
const MpiComm & inter0Comm() const
Access function for MpiComm communicator for processes with subRank() 0.
Definition: Environment.C:313
QUESO::ScalarSequence::subScaleForKde
T subScaleForKde(unsigned int initialPos, const T &iqrValue, unsigned int kdeDimension) const
Selects the scales (output value) for the kernel density estimation, considering only the sub-sequenc...
Definition: ScalarSequence.C:2167
QUESO::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2640
QUESO::ScalarSequence::seqScalarPositionConstIteratorTypedef
std::vector< T >::const_iterator seqScalarPositionConstIteratorTypedef
Definition: ScalarSequence.h:60
QUESO::ScalarSequence::unifiedPopulationVariance
T unifiedPopulationVariance(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, const T &unifiedMeanValue) const
Finds the population variance of the unified sequence, considering numPos positions starting at posit...
Definition: ScalarSequence.C:1213
QUESO::ScalarSequence::unifiedMeanCltStd
T unifiedMeanCltStd(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos, const T &unifiedMeanValue) const
Definition: ScalarSequence.C:3434
QUESO::ScalarSequence::unifiedMedianExtra
T unifiedMedianExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos) const
Finds the median value of the unified sequence, considering numPos positions starting at position ini...
Definition: ScalarSequence.C:965
QUESO::ScalarSequence::subHistogram
void subHistogram(unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, std::vector< T > &centers, std::vector< unsigned int > &bins) const
Calculates the histogram of the sub-sequence.
Definition: ScalarSequence.C:1569
QUESO::ScalarSequence::autoCorrViaFft
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
Definition: ScalarSequence.C:1331
QUESO::ScalarSequence::subBasicHistogram
void subBasicHistogram(unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, UniformOneDGrid< T > *&gridValues, std::vector< unsigned int > &bins) const
Calculates the histogram of the sub-sequence.
Definition: ScalarSequence.C:1706
QUESO::ScalarSequence::extractScalarSeq
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, ScalarSequence< T > &scalarSeq) const
Extracts a sequence of scalars.
Definition: ScalarSequence.C:3151
QUESO::ScalarSequence::setGaussian
void setGaussian(const T &mean, const T &stdDev)
Sets the values of the sequence as a Gaussian distribution of mean given by meanVec and standard devi...
Definition: ScalarSequence.C:515
QUESO::ScalarSequence::rawData
std::vector< T > & rawData()
The sequence of scalars. Access to private attribute m_seq.
Definition: ScalarSequence.C:3215
QUESO::ScalarSequence::psd
void psd(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, std::vector< double > &psdSequence) const
Definition: ScalarSequence.C:3530
QUESO::ScalarSequence::unifiedPositionsOfMaximum
T unifiedPositionsOfMaximum(const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &unifiedPositionsOfMaximum)
Finds the positions where the maximum element occurs in the unified sequence.
Definition: ScalarSequence.C:2507
QUESO::ScalarSequence::subCdfStacc
void subCdfStacc(unsigned int initialPos, std::vector< double > &cdfStaccValues, std::vector< double > &cdfStaccValuesup, std::vector< double > &cdfStaccValueslow, const ScalarSequence< T > &sortedDataValues) const
Definition: ScalarSequence.C:3808
QUESO::ScalarSequence::unifiedSort
void unifiedSort(bool useOnlyInter0Comm, unsigned int initialPos, ScalarSequence< T > &unifiedSortedSequence) const
Sorts the unified sequence of scalars.
Definition: ScalarSequence.C:1859
QUESO::ScalarSequence::setUniform
void setUniform(const T &a, const T &b)
Sets the values of the sequence as a uniform distribution between the values given by vectors aVec an...
Definition: ScalarSequence.C:527
QUESO::ScalarSequence::subUniformlySampledMdf
void subUniformlySampledMdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &mdfValues) const
Definition: ScalarSequence.C:3364
QUESO::ScalarSequence::subSort
void subSort()
Sorts the sequence of scalars in the private attribute m_seq.
Definition: ScalarSequence.C:3223
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
QUESO::MiscGaussianDensity
double MiscGaussianDensity(double x, double mu, double sigma)
Definition: Miscellaneous.C:306
QUESO::ScalarSequence::seqScalarPositionIteratorTypedef
std::vector< T >::iterator seqScalarPositionIteratorTypedef
Definition: ScalarSequence.h:59
QUESO::ScalarSequence::writeSubMatlabHeader
void writeSubMatlabHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write header info for matlab files from one chain.
Definition: ScalarSequence.C:2844
QUESO::ScalarSequence::subPopulationVariance
T subPopulationVariance(unsigned int initialPos, unsigned int numPos, const T &meanValue) const
Finds the population variance of the sub-sequence, considering numPos positions starting at position ...
Definition: ScalarSequence.C:1186
QUESO::ScalarSequence::clear
void clear()
Clears the sequence of scalars.
Definition: ScalarSequence.C:123
QUESO::ComputeUnifiedGaussian2dKde
void ComputeUnifiedGaussian2dKde(bool useOnlyInter0Comm, const ScalarSequence< T > &scalarSeq1, const ScalarSequence< T > &scalarSeq2, unsigned int initialPos, double unifiedScaleValue1, double unifiedScaleValue2, const std::vector< T > &unifiedEvaluationPositions1, const std::vector< T > &unifiedEvaluationPositions2, std::vector< double > &unifiedDensityValues)
Definition: ScalarSequence.C:3990
QUESO::ScalarSequence::deleteStoredScalars
void deleteStoredScalars()
Deletes all stored scalars.
Definition: ScalarSequence.C:467
QUESO::ScalarSequence::writeTxtHeader
void writeTxtHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write txt info for matlab files.
Definition: ScalarSequence.C:2856
QUESO::ScalarSequence::unifiedHistogram
void unifiedHistogram(bool useOnlyInter0Comm, unsigned int initialPos, const T &unifiedMinHorizontalValue, const T &unifiedMaxHorizontalValue, std::vector< T > &unifiedCenters, std::vector< unsigned int > &unifiedBins) const
Calculates the histogram of the unified sequence.
Definition: ScalarSequence.C:1616
QUESO::ScalarSequence::operator[]
const T & operator[](unsigned int posId) const
Access position posId of the sequence of scalars (const).
Definition: ScalarSequence.C:69
QUESO::ScalarSequence::subSampleVariancePlain
const T & subSampleVariancePlain() const
Finds the variance of a sample of the sub-sequence of scalars.
Definition: ScalarSequence.C:443
QUESO::ScalarSequence
Class for handling scalar samples.
Definition: ScalarSequence.h:54
QUESO::ScalarSequence::subMinMaxExtra
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
Definition: ScalarSequence.C:1469
QUESO::ScalarSequence::m_seq
std::vector< T > m_seq
Definition: ScalarSequence.h:527
QUESO::ScalarSequence::resetValues
void resetValues(unsigned int initialPos, unsigned int)
Sets numPos values of the sequence to zero, starting at position initialPos.
Definition: ScalarSequence.C:186
QUESO::ScalarSequence::subWeightHistogram
void subWeightHistogram(unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, UniformOneDGrid< T > *&gridValues, std::vector< unsigned int > &bins) const
Calculates the weighted histogram of the sub-sequence.
Definition: ScalarSequence.C:1748
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2388
QUESO::ScalarSequence::extractRawData
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, std::vector< double > &rawData) const
Extracts the raw data.
Definition: ScalarSequence.C:3192
QUESO::ScalarSequence::autoCovariance
T autoCovariance(unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
Calculates the autocovariance.
Definition: ScalarSequence.C:1272
QUESO::Fft
Class for a Fast Fourier Transform (FFT) algorithm.
Definition: Fft.h:66
QUESO::ComputeCovCorrBetweenScalarSequences
void ComputeCovCorrBetweenScalarSequences(const ScalarSequence< T > &subPSeq, const ScalarSequence< T > &subQSeq, unsigned int subNumSamples, T &covValue, T &corrValue)
Definition: ScalarSequence.C:4040
QUESO::ScalarSequence::unifiedMinMaxExtra
void unifiedMinMaxExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, T &unifiedMinValue, T &unifiedMaxValue) const
Finds the minimum and the maximum values of the unified sequence, considering numPos positions starti...
Definition: ScalarSequence.C:1498
QUESO::FilePtrSetStruct::ifsVar
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:89
QUESO::ScalarSequence::unifiedInterQuantileRange
T unifiedInterQuantileRange(bool useOnlyInter0Comm, unsigned int initialPos) const
Returns the interquartile range of the values in the unified sequence.
Definition: ScalarSequence.C:2076
QUESO::ScalarSequence::bmm
T bmm(unsigned int initialPos, unsigned int batchLength) const
Definition: ScalarSequence.C:3494
QUESO::ScalarSequence::subWeightCdf
void subWeightCdf(unsigned int numIntervals, std::vector< T > &gridValues, std::vector< T > &cdfValues) const
Finds the Weighted Cumulative Distribution Function (CDF) of the sub-sequence of scalars.
Definition: ScalarSequence.C:717
QUESO::ScalarSequence::subGaussian1dKde
void subGaussian1dKde(unsigned int initialPos, double scaleValue, const std::vector< T > &evaluationPositions, std::vector< double > &densityValues) const
Gaussian kernel for the KDE estimate of the sub-sequence.
Definition: ScalarSequence.C:2278
QUESO::ScalarSequence::unifiedUniformlySampledCdf
void unifiedUniformlySampledCdf(bool useOnlyInter0Comm, unsigned int numIntervals, T &unifiedMinDomainValue, T &unifiedMaxDomainValue, std::vector< T > &unifiedCdfValues) const
Uniformly samples from the CDF from the unified sequence.
Definition: ScalarSequence.C:581
QUESO::ScalarSequence::unifiedSampleVariancePlain
const T & unifiedSampleVariancePlain(bool useOnlyInter0Comm) const
Finds the variance of a sample of the unified sequence of scalars.
Definition: ScalarSequence.C:455
QUESO::queso_require_equal_to_msg
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
Definition: MonteCarloSGOptions.C:460
QUESO::ScalarSequence::subMeanCltStd
T subMeanCltStd(unsigned int initialPos, unsigned int numPos, const T &meanValue) const
Definition: ScalarSequence.C:3405
QUESO::ScalarSequence::unifiedMaxPlain
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
Definition: ScalarSequence.C:382
QUESO::ComputeSubGaussian2dKde
void ComputeSubGaussian2dKde(const ScalarSequence< T > &scalarSeq1, const ScalarSequence< T > &scalarSeq2, unsigned int initialPos, double scaleValue1, double scaleValue2, const std::vector< T > &evaluationPositions1, const std::vector< T > &evaluationPositions2, std::vector< double > &densityValues)
Definition: ScalarSequence.C:3931
QUESO::ScalarSequence::subSort
void subSort(unsigned int initialPos, ScalarSequence< T > &sortedSequence) const
Sorts the sub-sequence of scalars.
Definition: ScalarSequence.C:1842
QUESO::ScalarSequence::parallelMerge
void parallelMerge(std::vector< T > &sortedBuffer, const std::vector< T > &leafData, unsigned int treeLevel) const
Sorts/merges data in parallel using MPI.
Definition: ScalarSequence.C:3233
QUESO::ScalarSequence::subMedianExtra
T subMedianExtra(unsigned int initialPos, unsigned int numPos) const
Finds the median value of the sub-sequence, considering numPos positions starting at position initial...
Definition: ScalarSequence.C:922
QUESO::ScalarSequence::append
void append(const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
Appends the scalar sequence src to this sequence.
Definition: ScalarSequence.C:2456
QUESO::ScalarSequence::meanStacc
T meanStacc(unsigned int initialPos) const
Definition: ScalarSequence.C:3718
QUESO::UniformOneDGrid
Class for accommodating uniform one-dimensional grids.
Definition: UniformOneDGrid.h:41
QUESO::ScalarSequence::autoCorrViaDef
T autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag) const
Calculates the autocorrelation via definition.
Definition: ScalarSequence.C:1302
QUESO::ScalarSequence::subMedianPlain
const T & subMedianPlain() const
Finds the median value of the sub-sequence of scalars.
Definition: ScalarSequence.C:419
QUESO::ScalarSequence::subSampleStd
T subSampleStd(unsigned int initialPos, unsigned int numPos, const T &meanValue) const
Finds the sample standard deviation of the unified sequence, considering numPos positions starting at...
Definition: ScalarSequence.C:1098
QUESO::ScalarSequence::subWriteContents
void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Writes the sub-sequence to a file.
Definition: ScalarSequence.C:2537
QUESO::ScalarSequence::subCdfPercentageRange
void subCdfPercentageRange(unsigned int initialPos, unsigned int numPos, double range, T &lowerValue, T &upperValue) const
Definition: ScalarSequence.C:3730
QUESO::ScalarSequence::setName
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
Definition: ScalarSequence.C:115
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:86
QUESO::ScalarSequence::env
const BaseEnvironment & env() const
Access to QUESO environment.
Definition: ScalarSequence.C:101
QUESO::Fft::inverse
void inverse(const std::vector< T > &data, unsigned int fftSize, std::vector< std::complex< double > > &result)
Calculates the inverse Fourier transform for real and complex data.
QUESO::ScalarSequence::getUnifiedContentsAtProc0Only
void getUnifiedContentsAtProc0Only(bool useOnlyInter0Comm, std::vector< T > &outputVec) const
Gets the unified contents of processor of rank equals to 0.
Definition: ScalarSequence.C:235
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:77
QUESO::ScalarSequence::subUniformlySampledCdf
void subUniformlySampledCdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
Definition: ScalarSequence.C:539
QUESO::ScalarSequence::copy
void copy(const ScalarSequence< T > &src)
Copies the scalar sequence src to this.
Definition: ScalarSequence.C:3136
QUESO::ScalarSequence::subBasicCdf
void subBasicCdf(unsigned int numIntervals, UniformOneDGrid< T > *&gridValues, std::vector< T > &cdfValues) const
Finds the Cumulative Distribution Function (CDF) of the sub-sequence of scalars.
Definition: ScalarSequence.C:680
QUESO::BaseEnvironment::numSubEnvironments
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:335
QUESO::ScalarSequence::unifiedMedianPlain
const T & unifiedMedianPlain(bool useOnlyInter0Comm) const
Finds the median value of the unified sequence of scalars.
Definition: ScalarSequence.C:431
QUESO::ScalarSequence::geweke
T geweke(unsigned int initialPos, double ratioNa, double ratioNb) const
Definition: ScalarSequence.C:3653
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:282
QUESO::ScalarSequence::unifiedReadContents
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
Definition: ScalarSequence.C:2865
QUESO::RawType_MPI_Status
MPI_Status RawType_MPI_Status
Definition: MpiComm.h:46
QUESO::BaseEnvironment::inter0Rank
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
QUESO::ScalarSequence::subPositionsOfMaximum
T subPositionsOfMaximum(const ScalarSequence< T > &subCorrespondingScalarValues, ScalarSequence< T > &subPositionsOfMaximum)
Finds the positions where the maximum element occurs in the sub-sequence.
Definition: ScalarSequence.C:2477
QUESO::ScalarSequence::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:136
QUESO::ScalarSequence::subSampleVarianceExtra
T subSampleVarianceExtra(unsigned int initialPos, unsigned int numPos, const T &meanValue) const
Finds the sample variance of the sub-sequence, considering numPos positions starting at position init...
Definition: ScalarSequence.C:1012
QUESO::ScalarSequence::name
const std::string & name() const
Access to the name of the sequence of scalars.
Definition: ScalarSequence.C:108
Generated on Tue Jun 5 2018 19:48:54 for queso-0.57.1 by  doxygen 1.8.5