queso-0.57.0
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::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::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::subSampleVariancePlain
const T & subSampleVariancePlain() const
Finds the variance of a sample of the sub-sequence of scalars.
Definition: ScalarSequence.C:443
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::ScalarSequence::deleteStoredScalars
void deleteStoredScalars()
Deletes all stored scalars.
Definition: ScalarSequence.C:467
QUESO::ScalarSequence::subSort
void subSort(unsigned int initialPos, ScalarSequence< T > &sortedSequence) const
Sorts the sub-sequence of scalars.
Definition: ScalarSequence.C:1842
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::FilePtrSetStruct::ifsVar
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:89
QUESO::ScalarSequence::m_seq
std::vector< T > m_seq
Definition: ScalarSequence.h:527
QUESO::ScalarSequence::subMeanPlain
const T & subMeanPlain() const
Finds the mean value of the sub-sequence of scalars.
Definition: ScalarSequence.C:395
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::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::seqScalarPositionIteratorTypedef
std::vector< T >::iterator seqScalarPositionIteratorTypedef
Definition: ScalarSequence.h:59
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::ScalarSequence
Class for handling scalar samples.
Definition: ScalarSequence.h:54
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::ScalarSequence::unifiedMaxPlain
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
Definition: ScalarSequence.C:382
QUESO::ScalarSequence::~ScalarSequence
~ScalarSequence()
Destructor.
Definition: ScalarSequence.C:54
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::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::subSort
void subSort()
Sorts the sequence of scalars in the private attribute m_seq.
Definition: ScalarSequence.C:3223
dataIn
istream * dataIn
Definition: ann_sample.cpp:58
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:92
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::unifiedMedianPlain
const T & unifiedMedianPlain(bool useOnlyInter0Comm) const
Finds the median value of the unified sequence of scalars.
Definition: ScalarSequence.C:431
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::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
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::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::ScalarSequence::unifiedCdfPercentageRange
void unifiedCdfPercentageRange(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, double range, T &lowerValue, T &upperValue) const
Definition: ScalarSequence.C:3764
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::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::MiscHammingWindow
double MiscHammingWindow(unsigned int N, unsigned int j)
Definition: Miscellaneous.C:298
QUESO::Fft
Class for a Fast Fourier Transform (FFT) algorithm.
Definition: Fft.h:66
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::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:136
k
int k
Definition: ann_sample.cpp:53
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::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::subMedianPlain
const T & subMedianPlain() const
Finds the median value of the sub-sequence of scalars.
Definition: ScalarSequence.C:419
QUESO::UniformOneDGrid
Class for accommodating uniform one-dimensional grids.
Definition: UniformOneDGrid.h:41
QUESO::BaseEnvironment::inter0Comm
const MpiComm & inter0Comm() const
Access function for MpiComm communicator for processes with subRank() 0.
Definition: Environment.C:313
QUESO::ScalarSequence::unifiedMeanCltStd
T unifiedMeanCltStd(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int localNumPos, const T &unifiedMeanValue) const
Definition: ScalarSequence.C:3434
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:86
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::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::autoCovariance
T autoCovariance(unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
Calculates the autocovariance.
Definition: ScalarSequence.C:1272
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::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::env
const BaseEnvironment & env() const
Access to QUESO environment.
Definition: ScalarSequence.C:101
QUESO::BaseEnvironment::inter0Rank
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
QUESO::BaseEnvironment::worldRank
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
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::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::unifiedSort
void unifiedSort(bool useOnlyInter0Comm, unsigned int initialPos, ScalarSequence< T > &unifiedSortedSequence) const
Sorts the unified sequence of scalars.
Definition: ScalarSequence.C:1859
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::ScalarSequence::unifiedMeanPlain
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Definition: ScalarSequence.C:407
QUESO::ScalarSequence::subMaxPlain
const T & subMaxPlain() const
Finds the maximum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:369
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::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::unifiedMinPlain
const T & unifiedMinPlain(bool useOnlyInter0Comm) const
Finds the minimum value of the unified sequence of scalars.
Definition: ScalarSequence.C:356
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::MiscGaussianDensity
double MiscGaussianDensity(double x, double mu, double sigma)
Definition: Miscellaneous.C:306
QUESO::ScalarSequence::subUniformlySampledMdf
void subUniformlySampledMdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &mdfValues) const
Definition: ScalarSequence.C:3364
QUESO::ScalarSequence::subCdfPercentageRange
void subCdfPercentageRange(unsigned int initialPos, unsigned int numPos, double range, T &lowerValue, T &upperValue) const
Definition: ScalarSequence.C:3730
QUESO::ScalarSequence::rawData
std::vector< T > & rawData()
The sequence of scalars. Access to private attribute m_seq.
Definition: ScalarSequence.C:3215
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::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::ScalarSequence::clear
void clear()
Clears the sequence of scalars.
Definition: ScalarSequence.C:123
QUESO::ScalarSequence::append
void append(const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
Appends the scalar sequence src to this sequence.
Definition: ScalarSequence.C:2456
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::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::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::ComputeCovCorrBetweenScalarSequences
void ComputeCovCorrBetweenScalarSequences(const ScalarSequence< T > &subPSeq, const ScalarSequence< T > &subQSeq, unsigned int subNumSamples, T &covValue, T &corrValue)
Definition: ScalarSequence.C:4040
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::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::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::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::subMeanCltStd
T subMeanCltStd(unsigned int initialPos, unsigned int numPos, const T &meanValue) const
Definition: ScalarSequence.C:3405
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:77
QUESO::ScalarSequence::geweke
T geweke(unsigned int initialPos, double ratioNa, double ratioNb) const
Definition: ScalarSequence.C:3653
QUESO::ScalarSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
Definition: ScalarSequence.C:143
QUESO::ScalarSequence::setName
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
Definition: ScalarSequence.C:115
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
copy
and that you are informed that you can do these things To protect your we need to make restrictions that forbid distributors to deny you these rights or to ask you to surrender these rights These restrictions translate to certain responsibilities for you if you distribute copies of the library or if you modify it For if you distribute copies of the whether gratis or for a you must give the recipients all the rights that we gave you You must make sure that receive or can get the source code If you link other code with the you must provide complete object files to the so that they can relink them with the library after making changes to the library and recompiling it And you must show them these terms so they know their rights We protect your rights with a two step which gives you legal permission to copy
Definition: License.txt:60
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::meanStacc
T meanStacc(unsigned int initialPos) const
Definition: ScalarSequence.C:3718
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::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::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::RawType_MPI_Status
MPI_Status RawType_MPI_Status
Definition: MpiComm.h:46
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::operator=
ScalarSequence< T > & operator=(const ScalarSequence< T > &rhs)
Assignment operator; it copies rhs to this.
Definition: ScalarSequence.C:61
QUESO::ScalarSequence::psd
void psd(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, std::vector< double > &psdSequence) const
Definition: ScalarSequence.C:3530
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::copy
void copy(const ScalarSequence< T > &src)
Copies the scalar sequence src to this.
Definition: ScalarSequence.C:3136
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::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::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::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::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::name
const std::string & name() const
Access to the name of the sequence of scalars.
Definition: ScalarSequence.C:108
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::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::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::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::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:282
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2388
QUESO::ScalarSequence::subMinPlain
const T & subMinPlain() const
Finds the minimum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:343
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::ScalarSequence
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
Definition: ScalarSequence.C:32
QUESO::BaseEnvironment::numSubEnvironments
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:335
Generated on Sat Apr 22 2017 14:04:35 for queso-0.57.0 by  doxygen 1.8.5