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