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