queso-0.56.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  double tmp = log((double) numPos)/log(2.);
1370  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1371  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1372  unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
1373 
1374  std::vector<double> rawDataVec(numPos,0.);
1375  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1376  Fft<T> fftObj(m_env);
1377 
1378  // Forward FFT
1379  this->extractRawData(initialPos,
1380  1, // spacing
1381  numPos,
1382  rawDataVec);
1383  T meanValue = this->subMeanExtra(initialPos,
1384  numPos);
1385  for (unsigned int j = 0; j < numPos; ++j) {
1386  rawDataVec[j] -= meanValue; // IMPORTANT
1387  }
1388 
1389  rawDataVec.resize(fftSize,0.);
1390 
1391  //if (m_env.subDisplayFile()) {
1392  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1393  // << ": about to call fftObj.forward()"
1394  // << " with rawDataVec.size() = " << rawDataVec.size()
1395  // << ", fftSize = " << fftSize
1396  // << ", resultData.size() = " << resultData.size()
1397  // << std::endl;
1398  //}
1399  fftObj.forward(rawDataVec,fftSize,resultData);
1400 
1401  // Inverse FFT
1402  for (unsigned int j = 0; j < fftSize; ++j) {
1403  rawDataVec[j] = std::norm(resultData[j]);
1404  }
1405  //if (m_env.subDisplayFile()) {
1406  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1407  // << ": about to call fftObj.inverse()"
1408  // << " with rawDataVec.size() = " << rawDataVec.size()
1409  // << ", fftSize = " << fftSize
1410  // << ", resultData.size() = " << resultData.size()
1411  // << std::endl;
1412  //}
1413  fftObj.inverse(rawDataVec,fftSize,resultData);
1414  //if (m_env.subDisplayFile()) {
1415  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft()"
1416  // << ": returned succesfully from fftObj.inverse()"
1417  // << std::endl;
1418  //}
1419 
1420  // Prepare return data
1421  autoCorrs.resize(maxLag+1,0.); // Yes, +1
1422  for (unsigned int j = 0; j < autoCorrs.size(); ++j) {
1423  double ratio = ((double) j)/((double) (numPos-1));
1424  autoCorrs[j] = ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1425  }
1426 
1427  return;
1428 }
1429 // --------------------------------------------------
1430 template <class T>
1431 void
1432 ScalarSequence<T>::autoCorrViaFft(
1433  unsigned int initialPos,
1434  unsigned int numPos,
1435  unsigned int numSum,
1436  T& autoCorrsSum) const
1437 {
1438  //if (m_env.subDisplayFile()) {
1439  // *m_env.subDisplayFile() << "Entering ScalarSequence<T>::autoCorrViaFft(), for sum"
1440  // << ": initialPos = " << initialPos
1441  // << ", numPos = " << numPos
1442  // << std::endl;
1443  //}
1444 
1445  double tmp = log((double) numPos)/log(2.);
1446  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
1447  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
1448  unsigned int fftSize = (unsigned int) std::pow(2.,tmp+1);
1449 
1450  std::vector<double> rawDataVec(numPos,0.);
1451  std::vector<std::complex<double> > resultData(0,std::complex<double>(0.,0.));
1452  Fft<T> fftObj(m_env);
1453 
1454  // Forward FFT
1455  this->extractRawData(initialPos,
1456  1, // spacing
1457  numPos,
1458  rawDataVec);
1459  T meanValue = this->subMeanExtra(initialPos,
1460  numPos);
1461  for (unsigned int j = 0; j < numPos; ++j) {
1462  rawDataVec[j] -= meanValue; // IMPORTANT
1463  }
1464  rawDataVec.resize(fftSize,0.);
1465 
1466  //if (m_env.subDisplayFile()) {
1467  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
1468  // << ": about to call fftObj.forward()"
1469  // << " with rawDataVec.size() = " << rawDataVec.size()
1470  // << ", fftSize = " << fftSize
1471  // << ", resultData.size() = " << resultData.size()
1472  // << std::endl;
1473  //}
1474  fftObj.forward(rawDataVec,fftSize,resultData);
1475 
1476  // Inverse FFT
1477  for (unsigned int j = 0; j < fftSize; ++j) {
1478  rawDataVec[j] = std::norm(resultData[j]);
1479  }
1480  fftObj.inverse(rawDataVec,fftSize,resultData);
1481 
1482  //if (m_env.subDisplayFile()) {
1483  // *m_env.subDisplayFile() << "In ScalarSequence<T>::autoCorrViaFft(), for sum"
1484  // << ": computed auto covariance for lag 0 = " << resultData[0].real()/((double) (numPos))
1485  // << ", computed resultData[0].imag() = " << resultData[0].imag()
1486  // << std::endl;
1487  //}
1488 
1489  // Prepare return data
1490  autoCorrsSum = 0.;
1491  for (unsigned int j = 0; j < numSum; ++j) { // Yes, begin at lag '0'
1492  double ratio = ((double) j)/((double) (numPos-1));
1493  autoCorrsSum += ( resultData[j].real()/resultData[0].real() )*(1.-ratio);
1494  }
1495 
1496  return;
1497 }
1498 // --------------------------------------------------
1499 template <class T>
1500 void
1501 ScalarSequence<T>::subMinMaxExtra(
1502  unsigned int initialPos,
1503  unsigned int numPos,
1504  T& minValue,
1505  T& maxValue) const
1506 {
1507  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
1508 
1509  seqScalarPositionConstIteratorTypedef pos1 = m_seq.begin();
1510  std::advance(pos1,initialPos);
1511 
1512  seqScalarPositionConstIteratorTypedef pos2 = m_seq.begin();
1513  std::advance(pos2,initialPos+numPos);
1514 
1515  if ((initialPos+numPos) == this->subSequenceSize()) {
1516  queso_require_msg(!(pos2 != m_seq.end()), "invalid state");
1517  }
1518 
1519  seqScalarPositionConstIteratorTypedef pos;
1520  pos = std::min_element(pos1, pos2);
1521  minValue = *pos;
1522  pos = std::max_element(pos1, pos2);
1523  maxValue = *pos;
1524 
1525  return;
1526 }
1527 // --------------------------------------------------
1528 template <class T>
1529 void
1530 ScalarSequence<T>::unifiedMinMaxExtra(
1531  bool useOnlyInter0Comm,
1532  unsigned int initialPos,
1533  unsigned int numPos,
1534  T& unifiedMinValue,
1535  T& unifiedMaxValue) const
1536 {
1537  if (m_env.numSubEnvironments() == 1) {
1538  return this->subMinMaxExtra(initialPos,
1539  numPos,
1540  unifiedMinValue,
1541  unifiedMaxValue);
1542  }
1543 
1544  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
1545 
1546  if (useOnlyInter0Comm) {
1547  if (m_env.inter0Rank() >= 0) {
1548  // Find local min and max
1549  T minValue;
1550  T maxValue;
1551  this->subMinMaxExtra(initialPos,
1552  numPos,
1553  minValue,
1554  maxValue);
1555 
1556  // Get overall min
1557  std::vector<double> sendBuf(1,0.);
1558  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1559  sendBuf[i] = minValue;
1560  }
1561  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMinValue, (int) sendBuf.size(), RawValue_MPI_MIN,
1562  "ScalarSequence<T>::unifiedMinMaxExtra()",
1563  "failed MPI.Allreduce() for min");
1564 
1565  // Get overall max
1566  for (unsigned int i = 0; i < sendBuf.size(); ++i) {
1567  sendBuf[i] = maxValue;
1568  }
1569  m_env.inter0Comm().template Allreduce<double>(&sendBuf[0], &unifiedMaxValue, (int) sendBuf.size(), RawValue_MPI_MAX,
1570  "ScalarSequence<T>::unifiedMinMaxExtra()",
1571  "failed MPI.Allreduce() for max");
1572 
1573  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1574  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedMinMaxExtra()"
1575  << ": localMinValue = " << minValue
1576  << ", localMaxValue = " << maxValue
1577  << ", unifiedMinValue = " << unifiedMinValue
1578  << ", unifiedMaxValue = " << unifiedMaxValue
1579  << std::endl;
1580  }
1581  }
1582  else {
1583  // Node not in the 'inter0' communicator
1584  this->subMinMaxExtra(initialPos,
1585  numPos,
1586  unifiedMinValue,
1587  unifiedMaxValue);
1588  }
1589  }
1590  else {
1591  queso_error_msg("parallel vectors not supported yet");
1592  }
1593 
1594  //m_env.fullComm().Barrier();
1595 
1596  return;
1597 }
1598 // --------------------------------------------------
1599 template <class T>
1600 void
1601 ScalarSequence<T>::subHistogram(
1602  unsigned int initialPos,
1603  const T& minHorizontalValue,
1604  const T& maxHorizontalValue,
1605  std::vector<T>& centers,
1606  std::vector<unsigned int>& bins) const
1607 {
1608  queso_require_equal_to_msg(centers.size(), bins.size(), "vectors 'centers' and 'bins' have different sizes");
1609 
1610  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1611 
1612  if (initialPos) {}; // just to remove compiler warning
1613 
1614  for (unsigned int j = 0; j < bins.size(); ++j) {
1615  centers[j] = 0.;
1616  bins[j] = 0;
1617  }
1618 
1619  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1620 
1621  double minCenter = minHorizontalValue - horizontalDelta/2.;
1622  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1623  for (unsigned int j = 0; j < centers.size(); ++j) {
1624  double factor = ((double) j)/(((double) centers.size()) - 1.);
1625  centers[j] = (1. - factor) * minCenter + factor * maxCenter;
1626  }
1627 
1628  unsigned int dataSize = this->subSequenceSize();
1629  for (unsigned int j = 0; j < dataSize; ++j) {
1630  double value = m_seq[j];
1631  if (value < minHorizontalValue) {
1632  bins[0]++;
1633  }
1634  else if (value >= maxHorizontalValue) {
1635  bins[bins.size()-1]++;
1636  }
1637  else {
1638  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1639  bins[index]++;
1640  }
1641  }
1642 
1643  return;
1644 }
1645 // --------------------------------------------------
1646 template <class T>
1647 void
1648 ScalarSequence<T>::unifiedHistogram(
1649  bool useOnlyInter0Comm,
1650  unsigned int initialPos,
1651  const T& unifiedMinHorizontalValue,
1652  const T& unifiedMaxHorizontalValue,
1653  std::vector<T>& unifiedCenters,
1654  std::vector<unsigned int>& unifiedBins) const
1655 {
1656  if (m_env.numSubEnvironments() == 1) {
1657  return this->subHistogram(initialPos,
1658  unifiedMinHorizontalValue,
1659  unifiedMaxHorizontalValue,
1660  unifiedCenters,
1661  unifiedBins);
1662  }
1663 
1664  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
1665 
1666  if (useOnlyInter0Comm) {
1667  if (m_env.inter0Rank() >= 0) {
1668  queso_require_equal_to_msg(unifiedCenters.size(), unifiedBins.size(), "vectors 'unifiedCenters' and 'unifiedBins' have different sizes");
1669 
1670  queso_require_greater_equal_msg(unifiedBins.size(), 3, "number of 'unifiedBins' is too small: should be at least 3");
1671 
1672  for (unsigned int j = 0; j < unifiedBins.size(); ++j) {
1673  unifiedCenters[j] = 0.;
1674  unifiedBins[j] = 0;
1675  }
1676 
1677  double unifiedHorizontalDelta = (unifiedMaxHorizontalValue - unifiedMinHorizontalValue)/(((double) unifiedBins.size()) - 2.); // IMPORTANT: -2
1678 
1679  double unifiedMinCenter = unifiedMinHorizontalValue - unifiedHorizontalDelta/2.;
1680  double unifiedMaxCenter = unifiedMaxHorizontalValue + unifiedHorizontalDelta/2.;
1681  for (unsigned int j = 0; j < unifiedCenters.size(); ++j) {
1682  double factor = ((double) j)/(((double) unifiedCenters.size()) - 1.);
1683  unifiedCenters[j] = (1. - factor) * unifiedMinCenter + factor * unifiedMaxCenter;
1684  }
1685 
1686  std::vector<unsigned int> localBins(unifiedBins.size(),0);
1687  unsigned int dataSize = this->subSequenceSize();
1688  for (unsigned int j = 0; j < dataSize; ++j) {
1689  double value = m_seq[j];
1690  if (value < unifiedMinHorizontalValue) {
1691  localBins[0]++;
1692  }
1693  else if (value >= unifiedMaxHorizontalValue) {
1694  localBins[localBins.size()-1]++;
1695  }
1696  else {
1697  unsigned int index = 1 + (unsigned int) ((value - unifiedMinHorizontalValue)/unifiedHorizontalDelta);
1698  localBins[index]++;
1699  }
1700  }
1701 
1702  m_env.inter0Comm().template Allreduce<unsigned int>(&localBins[0], &unifiedBins[0], (int) localBins.size(), RawValue_MPI_SUM,
1703  "ScalarSequence<T>::unifiedHistogram()",
1704  "failed MPI.Allreduce() for bins");
1705 
1706  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1707  for (unsigned int i = 0; i < unifiedCenters.size(); ++i) {
1708  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedHistogram()"
1709  << ": i = " << i
1710  << ", unifiedMinHorizontalValue = " << unifiedMinHorizontalValue
1711  << ", unifiedMaxHorizontalValue = " << unifiedMaxHorizontalValue
1712  << ", unifiedCenters = " << unifiedCenters[i]
1713  << ", unifiedBins = " << unifiedBins[i]
1714  << std::endl;
1715  }
1716  }
1717  }
1718  else {
1719  // Node not in the 'inter0' communicator
1720  this->subHistogram(initialPos,
1721  unifiedMinHorizontalValue,
1722  unifiedMaxHorizontalValue,
1723  unifiedCenters,
1724  unifiedBins);
1725  }
1726  }
1727  else {
1728  queso_error_msg("parallel vectors not supported yet");
1729  }
1730 
1731  //m_env.fullComm().Barrier();
1732 
1733  return;
1734 }
1735 // --------------------------------------------------
1736 template <class T>
1737 void
1738 ScalarSequence<T>::subBasicHistogram(
1739  unsigned int initialPos,
1740  const T& minHorizontalValue,
1741  const T& maxHorizontalValue,
1742  UniformOneDGrid<T>*& gridValues,
1743  std::vector<unsigned int>& bins) const
1744 {
1745  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1746 
1747  for (unsigned int j = 0; j < bins.size(); ++j) {
1748  bins[j] = 0;
1749  }
1750 
1751  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1752  double minCenter = minHorizontalValue - horizontalDelta/2.;
1753  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1754  gridValues = new UniformOneDGrid<T>(m_env,
1755  "",
1756  bins.size(),
1757  minCenter,
1758  maxCenter);
1759 
1760  unsigned int dataSize = this->subSequenceSize();
1761  for (unsigned int j = 0; j < dataSize; ++j) {
1762  double value = m_seq[j];
1763  if (value < minHorizontalValue) {
1764  bins[0] += value;
1765  }
1766  else if (value >= maxHorizontalValue) {
1767  bins[bins.size()-1] += value;
1768  }
1769  else {
1770  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1771  bins[index] += value;
1772  }
1773  }
1774 
1775  return;
1776 }
1777 // --------------------------------------------------
1778 template <class T>
1779 void
1780 ScalarSequence<T>::subWeightHistogram(
1781  unsigned int initialPos,
1782  const T& minHorizontalValue,
1783  const T& maxHorizontalValue,
1784  UniformOneDGrid<T>*& gridValues,
1785  std::vector<unsigned int>& bins) const
1786 {
1787  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1788 
1789  for (unsigned int j = 0; j < bins.size(); ++j) {
1790  bins[j] = 0;
1791  }
1792 
1793  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1794  double minCenter = minHorizontalValue - horizontalDelta/2.;
1795  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1796  gridValues = new UniformOneDGrid<T>(m_env,
1797  "",
1798  bins.size(),
1799  minCenter,
1800  maxCenter);
1801 
1802  unsigned int dataSize = this->subSequenceSize();
1803  for (unsigned int j = 0; j < dataSize; ++j) {
1804  double value = m_seq[j];
1805  if (value < minHorizontalValue) {
1806  bins[0] += value;
1807  }
1808  else if (value >= maxHorizontalValue) {
1809  bins[bins.size()-1] += value;
1810  }
1811  else {
1812  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1813  bins[index] += value;
1814  }
1815  }
1816 
1817  return;
1818 }
1819 // --------------------------------------------------
1820 template <class T>
1821 void
1822 ScalarSequence<T>::subWeightHistogram(
1823  unsigned int initialPos,
1824  const T& minHorizontalValue,
1825  const T& maxHorizontalValue,
1826  std::vector<T>& gridValues,
1827  std::vector<unsigned int>& bins) const
1828 {
1829  queso_require_greater_equal_msg(bins.size(), 3, "number of 'bins' is too small: should be at least 3");
1830 
1831  for (unsigned int j = 0; j < bins.size(); ++j) {
1832  bins[j] = 0;
1833  }
1834 
1835  double horizontalDelta = (maxHorizontalValue - minHorizontalValue)/(((double) bins.size()) - 2.); // IMPORTANT: -2
1836  double minCenter = minHorizontalValue - horizontalDelta/2.;
1837  double maxCenter = maxHorizontalValue + horizontalDelta/2.;
1838  UniformOneDGrid<T> tmpGrid(m_env,
1839  "",
1840  bins.size(),
1841  minCenter,
1842  maxCenter);
1843  gridValues.clear();
1844  gridValues.resize(tmpGrid.size(),0.);
1845  for (unsigned int i = 0; i < tmpGrid.size(); ++i) {
1846  gridValues[i] = tmpGrid[i];
1847  }
1848 
1849  unsigned int dataSize = this->subSequenceSize();
1850  for (unsigned int j = 0; j < dataSize; ++j) {
1851  double value = m_seq[j];
1852  if (value < minHorizontalValue) {
1853  bins[0] += value;
1854  }
1855  else if (value >= maxHorizontalValue) {
1856  bins[bins.size()-1] += value;
1857  }
1858  else {
1859  unsigned int index = 1 + (unsigned int) ((value - minHorizontalValue)/horizontalDelta);
1860  bins[index] += value;
1861  }
1862  }
1863 
1864  //std::cout << "In ScalarSequence<T>::subWeightHistogram():" << std::endl;
1865  //for (unsigned int j = 0; j < bins.size(); ++j) {
1866  // std::cout << "bins[" << j << "] = " << bins[j] << std::endl;
1867  //}
1868 
1869  return;
1870 }
1871 // --------------------------------------------------
1872 template <class T>
1873 void
1874 ScalarSequence<T>::subSort(
1875  unsigned int initialPos,
1876  ScalarSequence<T>& sortedSequence) const
1877 {
1878  unsigned int numPos = this->subSequenceSize() - initialPos;
1879  sortedSequence.resizeSequence(numPos);
1880  this->extractScalarSeq(initialPos,
1881  1,
1882  numPos,
1883  sortedSequence);
1884  sortedSequence.subSort();
1885 
1886  return;
1887 }
1888 // --------------------------------------------------
1889 template <class T>
1890 void
1891 ScalarSequence<T>::unifiedSort(
1892  bool useOnlyInter0Comm,
1893  unsigned int initialPos,
1894  ScalarSequence<T>& unifiedSortedSequence) const
1895 {
1896  if (m_env.numSubEnvironments() == 1) {
1897  return this->subSort(initialPos,unifiedSortedSequence);
1898  }
1899 
1900  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
1901 
1902  if (useOnlyInter0Comm) {
1903  if (m_env.inter0Rank() >= 0) {
1904  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1905 
1906  unsigned int localNumPos = this->subSequenceSize() - initialPos;
1907 
1908  std::vector<T> leafData(localNumPos,0.);
1909  this->extractRawData(0,
1910  1,
1911  localNumPos,
1912  leafData);
1913 
1914  if (m_env.inter0Rank() == 0) {
1915  int minus1NumTreeLevels = 0;
1916  int power2NumTreeNodes = 1;
1917 
1918  while (power2NumTreeNodes < m_env.inter0Comm().NumProc()) {
1919  power2NumTreeNodes += power2NumTreeNodes;
1920  minus1NumTreeLevels++;
1921  }
1922 
1923  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1924  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedSort()"
1925  << ": sorting tree has " << m_env.inter0Comm().NumProc()
1926  << " nodes and " << minus1NumTreeLevels+1
1927  << " levels"
1928  << std::endl;
1929  }
1930 
1931  this->parallelMerge(unifiedSortedSequence.rawData(),
1932  leafData,
1933  minus1NumTreeLevels);
1934  }
1935  else if (m_env.inter0Rank() > 0) { // KAUST
1936  unsigned int uintBuffer[1];
1937  RawType_MPI_Status status;
1938  m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, RawValue_MPI_ANY_SOURCE, SCALAR_SEQUENCE_INIT_MPI_MSG, &status,
1939  "ScalarSequence<T>::unifiedSort()",
1940  "failed MPI.Recv() for init");
1941  //if (status) {}; // just to remove compiler warning
1942 
1943  unsigned int treeLevel = uintBuffer[0];
1944  this->parallelMerge(unifiedSortedSequence.rawData(),
1945  leafData,
1946  treeLevel);
1947  }
1948 
1949  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), returned from parallelMerge()",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1950 
1951  // Broadcast
1952  unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
1953  m_env.inter0Comm().Bcast((void *) &unifiedDataSize, (int) 1, RawValue_MPI_UNSIGNED, 0,
1954  "ScalarSequence<T>::unifiedSort()",
1955  "failed MPI.Bcast() for unified data size");
1956 
1957  unsigned int sumOfNumPos = 0;
1958  m_env.inter0Comm().template Allreduce<unsigned int>(&localNumPos, &sumOfNumPos, (int) 1, RawValue_MPI_SUM,
1959  "ScalarSequence<T>::unifiedSort()",
1960  "failed MPI.Allreduce() for data size");
1961 
1962  queso_require_equal_to_msg(sumOfNumPos, unifiedDataSize, "incompatible unified sizes");
1963 
1964  unifiedSortedSequence.resizeSequence(unifiedDataSize);
1965  m_env.inter0Comm().Bcast((void *) &unifiedSortedSequence.rawData()[0], (int) unifiedDataSize, RawValue_MPI_DOUBLE, 0,
1966  "ScalarSequence<T>::unifiedSort()",
1967  "failed MPI.Bcast() for unified data");
1968 
1969  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
1970  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
1971  << ": tree node " << m_env.inter0Rank()
1972  << ", unifiedSortedSequence[0] = " << unifiedSortedSequence[0]
1973  << ", unifiedSortedSequence[" << unifiedSortedSequence.subSequenceSize()-1 << "] = " << unifiedSortedSequence[unifiedSortedSequence.subSequenceSize()-1]
1974  << std::endl;
1975  }
1976 
1977  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedSort(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
1978  }
1979  else {
1980  // Node not in the 'inter0' communicator
1981  this->subSort(initialPos,unifiedSortedSequence);
1982  }
1983  }
1984  else {
1985  queso_error_msg("parallel vectors not supported yet");
1986  }
1987 
1988  //m_env.fullComm().Barrier();
1989 
1990  return;
1991 }
1992 // --------------------------------------------------
1993 template <class T>
1994 T
1995 ScalarSequence<T>::subInterQuantileRange(unsigned int initialPos) const
1996 {
1997  queso_require_less_msg(initialPos, this->subSequenceSize(), "'initialPos' is too big");
1998 
1999  ScalarSequence sortedSequence(m_env,0,"");
2000  this->subSort(initialPos,
2001  sortedSequence);
2002 
2003  // The test above guarantees that 'dataSize >= 1'
2004  unsigned int dataSize = this->subSequenceSize() - initialPos;
2005 
2006  queso_require_equal_to_msg(dataSize, sortedSequence.subSequenceSize(), "inconsistent size variables");
2007 
2008  bool everythingOk = true;
2009 
2010  // pos1 = (dataSize+1)/4 - 1
2011  // pos1 >= 0 <==> dataSize >= 3
2012  // pos1 < (dataSize-1) <==> 3*dataSize > 1
2013  unsigned int pos1 = (unsigned int) ( (((double) dataSize) + 1.)*1./4. - 1. );
2014  if (pos1 > (dataSize-1)) {
2015  pos1 = 0;
2016  everythingOk = false;
2017  }
2018  unsigned int pos1inc = pos1+1;
2019  if (pos1inc > (dataSize-1)) {
2020  pos1inc = dataSize-1;
2021  everythingOk = false;
2022  }
2023 
2024  // pos3 = (dataSize+1)*3/4 - 1
2025  // pos3 >= 0 <==> dataSize >= 1/3
2026  // pos3 < (dataSize-1) <==> dataSize > 3
2027  unsigned int pos3 = (unsigned int) ( (((double) dataSize) + 1.)*3./4. - 1. );
2028  if (pos3 > (dataSize-1)) {
2029  pos3 = 0;
2030  everythingOk = false;
2031  }
2032  unsigned int pos3inc = pos3+1;
2033  if (pos3inc > (dataSize-1)) {
2034  pos3inc = dataSize-1;
2035  everythingOk = false;
2036  }
2037 
2038  double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
2039  if (fraction1 < 0.) {
2040  fraction1 = 0.;
2041  everythingOk = false;
2042  }
2043  double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
2044  if (fraction3 < 0.) {
2045  fraction3 = 0.;
2046  everythingOk = false;
2047  }
2048 
2049  if (everythingOk == false) {
2050  std::cerr << "In ScalarSequence<T>::subInterQuantileRange()"
2051  << ", worldRank = " << m_env.worldRank()
2052  << ": at least one adjustment was necessary"
2053  << std::endl;
2054  }
2055 
2056  //if (m_env.subDisplayFile()) {
2057  // *m_env.subDisplayFile() << "In ScalarSequence::subInterQuantileRange()"
2058  // << ", initialPos = " << initialPos
2059  // << ", this->subSequenceSize() = " << this->subSequenceSize()
2060  // << ", dataSize = " << dataSize
2061  // << ", sortedSequence.size() = " << sortedSequence.size()
2062  // << ", pos1 = " << pos1
2063  // << ", pos3 = " << pos3
2064  // << std::endl;
2065  //}
2066 
2067  T value1 = (1.-fraction1) * sortedSequence[pos1] + fraction1 * sortedSequence[pos1inc];
2068  T value3 = (1.-fraction3) * sortedSequence[pos3] + fraction3 * sortedSequence[pos3inc];
2069  T iqrValue = value3 - value1;
2070 
2071  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2072  *m_env.subDisplayFile() << "In ScalarSequence<T>::subInterQuantileRange()"
2073  << ": iqrValue = " << iqrValue
2074  << ", dataSize = " << dataSize
2075  << ", pos1 = " << pos1
2076  << ", pos3 = " << pos3
2077  << ", value1 = " << value1
2078  << ", value3 = " << value3
2079  << std::endl;
2080 
2081  // Save data only once into a separate file
2082  //std::ofstream* ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
2083  //if ((ofsvar == NULL ) ||
2084  // (ofsvar->is_open() == false)) {
2085  // delete ofsvar;
2086  // ofsvar = new std::ofstream(("sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
2087 
2088  // *ofsvar << "var_sort_sub" << m_env.subIdString() << " = zeros(" << 1
2089  // << "," << dataSize
2090  // << ");"
2091  // << std::endl;
2092  // for (unsigned int j = 0; j < dataSize; ++j) {
2093  // *ofsvar << "var_sort_sub" << m_env.subIdString() << "(" << 1
2094  // << "," << j+1
2095  // << ") = " << sortedSequence[j]
2096  // << ";"
2097  // << std::endl;
2098  // }
2099  //}
2100  //delete ofsvar;
2101  }
2102 
2103  return iqrValue;
2104 }
2105 // --------------------------------------------------
2106 template <class T>
2107 T
2108 ScalarSequence<T>::unifiedInterQuantileRange(
2109  bool useOnlyInter0Comm,
2110  unsigned int initialPos) const
2111 {
2112  T unifiedIqrValue = 0.;
2113 
2114  if (m_env.numSubEnvironments() == 1) {
2115  return this->subInterQuantileRange(initialPos);
2116  }
2117 
2118  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2119 
2120  if (useOnlyInter0Comm) {
2121  if (m_env.inter0Rank() >= 0) {
2122  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), beginning logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
2123 
2124  ScalarSequence unifiedSortedSequence(m_env,0,"");
2125  this->unifiedSort(useOnlyInter0Comm,
2126  initialPos,
2127  unifiedSortedSequence);
2128  unsigned int unifiedDataSize = unifiedSortedSequence.subSequenceSize();
2129 
2130  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2131  unsigned int sumOfLocalSizes = 0;
2132  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &sumOfLocalSizes, (int) 1, RawValue_MPI_SUM,
2133  "ScalarSequence<T>::unifiedInterQuantileRange()",
2134  "failed MPI.Allreduce() for data size");
2135 
2136  queso_require_equal_to_msg(sumOfLocalSizes, unifiedDataSize, "incompatible unified sizes");
2137 
2138  unsigned int pos1 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*1./4. - 1. );
2139  unsigned int pos3 = (unsigned int) ( (((double) unifiedDataSize) + 1.)*3./4. - 1. );
2140 
2141  double fraction1 = (((double) unifiedDataSize) + 1.)*1./4. - 1. - ((double) pos1);
2142  double fraction3 = (((double) unifiedDataSize) + 1.)*3./4. - 1. - ((double) pos3);
2143 
2144  T value1 = (1.-fraction1) * unifiedSortedSequence[pos1] + fraction1 * unifiedSortedSequence[pos1+1];
2145  T value3 = (1.-fraction3) * unifiedSortedSequence[pos3] + fraction3 * unifiedSortedSequence[pos3+1];
2146  unifiedIqrValue = value3 - value1;
2147 
2148  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2149  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedInterQuantileRange()"
2150  << ": unifiedIqrValue = " << unifiedIqrValue
2151  << ", localDataSize = " << localDataSize
2152  << ", unifiedDataSize = " << unifiedDataSize
2153  << ", pos1 = " << pos1
2154  << ", pos3 = " << pos3
2155  << ", value1 = " << value1
2156  << ", value3 = " << value3
2157  << std::endl;
2158 
2159  // Save data only once into a separate file
2160  //std::ofstream* ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::in | std::ofstream::ate);
2161  //if ((ofsvar == NULL ) ||
2162  // (ofsvar->is_open() == false)) {
2163  // delete ofsvar;
2164  // ofsvar = new std::ofstream(("unif_sort_sub"+m_env.subIdString()+".m").c_str(), std::ofstream::out | std::ofstream::trunc);
2165 
2166  // *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << " = zeros(" << 1
2167  // << "," << unifiedDataSize
2168  // << ");"
2169  // << std::endl;
2170  // for (unsigned int j = 0; j < unifiedDataSize; ++j) {
2171  // *ofsvar << "var_unif_sort_sub" << m_env.subIdString() << "(" << 1
2172  // << "," << j+1
2173  // << ") = " << unifiedSortedSequence[j]
2174  // << ";"
2175  // << std::endl;
2176  // }
2177  //}
2178  //delete ofsvar;
2179  }
2180 
2181  //m_env.syncPrintDebugMsg("In ScalarSequence<T>::unifiedInterQuantileRange(), ending logic",3,3000000,m_env.inter0Comm()); // Dangerous to barrier on inter0Comm ... // KAUST
2182  }
2183  else {
2184  // Node not in the 'inter0' communicator
2185  unifiedIqrValue = this->subInterQuantileRange(initialPos);
2186  }
2187  }
2188  else {
2189  queso_error_msg("parallel vectors not supported yet");
2190  }
2191 
2192  //m_env.fullComm().Barrier();
2193 
2194  return unifiedIqrValue;
2195 }
2196 // --------------------------------------------------
2197 template <class T>
2198 T
2199 ScalarSequence<T>::subScaleForKde(
2200  unsigned int initialPos,
2201  const T& iqrValue,
2202  unsigned int kdeDimension) const
2203 {
2204  bool bRC = (initialPos < this->subSequenceSize());
2205  queso_require_msg(bRC, "invalid input data");
2206 
2207  unsigned int dataSize = this->subSequenceSize() - initialPos;
2208 
2209  T meanValue = this->subMeanExtra(initialPos,
2210  dataSize);
2211 
2212  T samValue = this->subSampleVarianceExtra(initialPos,
2213  dataSize,
2214  meanValue);
2215 
2216  T scaleValue;
2217  if (iqrValue <= 0.) {
2218  scaleValue = 1.06*std::sqrt(samValue)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2219  }
2220  else {
2221  scaleValue = 1.06*std::min(std::sqrt(samValue),iqrValue/1.34)/std::pow(dataSize,1./(4. + ((double) kdeDimension)));
2222  }
2223 
2224  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2225  *m_env.subDisplayFile() << "In ScalarSequence<T>::subScaleForKde()"
2226  << ": iqrValue = " << iqrValue
2227  << ", meanValue = " << meanValue
2228  << ", samValue = " << samValue
2229  << ", dataSize = " << dataSize
2230  << ", scaleValue = " << scaleValue
2231  << std::endl;
2232  }
2233 
2234  return scaleValue;
2235 }
2236 // --------------------------------------------------
2237 template <class T>
2238 T
2239 ScalarSequence<T>::unifiedScaleForKde(
2240  bool useOnlyInter0Comm,
2241  unsigned int initialPos,
2242  const T& unifiedIqrValue,
2243  unsigned int kdeDimension) const
2244 {
2245  if (m_env.numSubEnvironments() == 1) {
2246  return this->subScaleForKde(initialPos,
2247  unifiedIqrValue,
2248  kdeDimension);
2249  }
2250 
2251  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2252 
2253  T unifiedScaleValue = 0.;
2254  if (useOnlyInter0Comm) {
2255  if (m_env.inter0Rank() >= 0) {
2256  bool bRC = (initialPos < this->subSequenceSize());
2257  queso_require_msg(bRC, "invalid input data");
2258 
2259  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2260 
2261  T unifiedMeanValue = this->unifiedMeanExtra(useOnlyInter0Comm,
2262  initialPos,
2263  localDataSize);
2264 
2265  T unifiedSamValue = this->unifiedSampleVarianceExtra(useOnlyInter0Comm,
2266  initialPos,
2267  localDataSize,
2268  unifiedMeanValue);
2269 
2270  unsigned int unifiedDataSize = 0;
2271  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &unifiedDataSize, (int) 1, RawValue_MPI_SUM,
2272  "ScalarSequence<T>::unifiedScaleForKde()",
2273  "failed MPI.Allreduce() for data size");
2274 
2275  if (unifiedIqrValue <= 0.) {
2276  unifiedScaleValue = 1.06*std::sqrt(unifiedSamValue)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
2277  }
2278  else {
2279  unifiedScaleValue = 1.06*std::min(std::sqrt(unifiedSamValue),unifiedIqrValue/1.34)/std::pow(unifiedDataSize,1./(4. + ((double) kdeDimension)));
2280  }
2281 
2282  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2283  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedScaleForKde()"
2284  << ": unifiedIqrValue = " << unifiedIqrValue
2285  << ", unifiedMeanValue = " << unifiedMeanValue
2286  << ", unifiedSamValue = " << unifiedSamValue
2287  << ", unifiedDataSize = " << unifiedDataSize
2288  << ", unifiedScaleValue = " << unifiedScaleValue
2289  << std::endl;
2290  }
2291  }
2292  else {
2293  // Node not in the 'inter0' communicator
2294  unifiedScaleValue = this->subScaleForKde(initialPos,
2295  unifiedIqrValue,
2296  kdeDimension);
2297  }
2298  }
2299  else {
2300  queso_error_msg("parallel vectors not supported yet");
2301  }
2302 
2303  //m_env.fullComm().Barrier();
2304 
2305  return unifiedScaleValue;
2306 }
2307 // --------------------------------------------------
2308 template <class T>
2309 void
2310 ScalarSequence<T>::subGaussian1dKde(
2311  unsigned int initialPos,
2312  double scaleValue,
2313  const std::vector<T>& evaluationPositions,
2314  std::vector<double>& densityValues) const
2315 {
2316  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2317  (0 < evaluationPositions.size()) &&
2318  (evaluationPositions.size() == densityValues.size() ));
2319  queso_require_msg(bRC, "invalid input data");
2320 
2321  unsigned int dataSize = this->subSequenceSize() - initialPos;
2322  unsigned int numEvals = evaluationPositions.size();
2323 
2324  double scaleInv = 1./scaleValue;
2325  for (unsigned int j = 0; j < numEvals; ++j) {
2326  double x = evaluationPositions[j];
2327  double value = 0.;
2328  for (unsigned int k = 0; k < dataSize; ++k) {
2329  double xk = m_seq[initialPos+k];
2330  value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
2331  }
2332  densityValues[j] = scaleInv * (value/(double) dataSize);
2333  }
2334 
2335  return;
2336 }
2337 // --------------------------------------------------
2338 template <class T>
2339 void
2340 ScalarSequence<T>::unifiedGaussian1dKde(
2341  bool useOnlyInter0Comm,
2342  unsigned int initialPos,
2343  double unifiedScaleValue,
2344  const std::vector<T>& unifiedEvaluationPositions,
2345  std::vector<double>& unifiedDensityValues) const
2346 {
2347  if (m_env.numSubEnvironments() == 1) {
2348  return this->subGaussian1dKde(initialPos,
2349  unifiedScaleValue,
2350  unifiedEvaluationPositions,
2351  unifiedDensityValues);
2352  }
2353 
2354  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
2355 
2356  if (useOnlyInter0Comm) {
2357  if (m_env.inter0Rank() >= 0) {
2358  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2359  (0 < unifiedEvaluationPositions.size()) &&
2360  (unifiedEvaluationPositions.size() == unifiedDensityValues.size() ));
2361  queso_require_msg(bRC, "invalid input data");
2362 
2363  unsigned int localDataSize = this->subSequenceSize() - initialPos;
2364  unsigned int unifiedDataSize = 0;
2365  m_env.inter0Comm().template Allreduce<unsigned int>(&localDataSize, &unifiedDataSize, (int) 1, RawValue_MPI_SUM,
2366  "ScalarSequence<T>::unifiedGaussian1dKde()",
2367  "failed MPI.Allreduce() for data size");
2368 
2369  unsigned int numEvals = unifiedEvaluationPositions.size();
2370 
2371  std::vector<double> densityValues(numEvals,0.);
2372  double unifiedScaleInv = 1./unifiedScaleValue;
2373  for (unsigned int j = 0; j < numEvals; ++j) {
2374  double x = unifiedEvaluationPositions[j];
2375  double value = 0.;
2376  for (unsigned int k = 0; k < localDataSize; ++k) {
2377  double xk = m_seq[initialPos+k];
2378  value += MiscGaussianDensity((x-xk)*unifiedScaleInv,0.,1.);
2379  }
2380  densityValues[j] = value;
2381  }
2382 
2383  for (unsigned int j = 0; j < numEvals; ++j) {
2384  unifiedDensityValues[j] = 0.;
2385  }
2386  m_env.inter0Comm().template Allreduce<double>(&densityValues[0], &unifiedDensityValues[0], (int) numEvals, RawValue_MPI_SUM,
2387  "ScalarSequence<T>::unifiedGaussian1dKde()",
2388  "failed MPI.Allreduce() for density values");
2389 
2390  for (unsigned int j = 0; j < numEvals; ++j) {
2391  unifiedDensityValues[j] *= unifiedScaleInv/((double) unifiedDataSize);
2392  }
2393 
2394  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
2395  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedGaussian1dKde()"
2396  << ": unifiedDensityValues[0] = " << unifiedDensityValues[0]
2397  << ", unifiedDensityValues[" << unifiedDensityValues.size()-1 << "] = " << unifiedDensityValues[unifiedDensityValues.size()-1]
2398  << std::endl;
2399  }
2400  }
2401  else {
2402  // Node not in the 'inter0' communicator
2403  this->subGaussian1dKde(initialPos,
2404  unifiedScaleValue,
2405  unifiedEvaluationPositions,
2406  unifiedDensityValues);
2407  }
2408  }
2409  else {
2410  queso_error_msg("parallel vectors not supported yet");
2411  }
2412 
2413  //m_env.fullComm().Barrier();
2414 
2415  return;
2416 }
2417 // --------------------------------------------------
2418 template <class T>
2419 void
2420 ScalarSequence<T>::filter(
2421  unsigned int initialPos,
2422  unsigned int spacing)
2423 {
2424  if (m_env.subDisplayFile()) {
2425  *m_env.subDisplayFile() << "Entering ScalarSequence<V,M>::filter()"
2426  << ": initialPos = " << initialPos
2427  << ", spacing = " << spacing
2428  << ", subSequenceSize = " << this->subSequenceSize()
2429  << std::endl;
2430  }
2431 
2432  unsigned int i = 0;
2433  unsigned int j = initialPos;
2434  unsigned int originalSubSequenceSize = this->subSequenceSize();
2435  while (j < originalSubSequenceSize) {
2436  if (i != j) {
2437  //*m_env.subDisplayFile() << i << "--" << j << " ";
2438  m_seq[i] = m_seq[j];
2439  }
2440  i++;
2441  j += spacing;
2442  }
2443 
2444  this->resizeSequence(i);
2445 
2446  if (m_env.subDisplayFile()) {
2447  *m_env.subDisplayFile() << "Leaving ScalarSequence<V,M>::filter()"
2448  << ": initialPos = " << initialPos
2449  << ", spacing = " << spacing
2450  << ", subSequenceSize = " << this->subSequenceSize()
2451  << std::endl;
2452  }
2453 
2454  return;
2455 }
2456 // --------------------------------------------------
2457 template <class T>
2458 T
2459 ScalarSequence<T>::brooksGelmanConvMeasure(
2460  bool useOnlyInter0Comm,
2461  unsigned int initialPos,
2462  unsigned int spacing) const
2463 {
2464  double resultValue = 0.;
2465 
2466  // FIX ME: put logic if (numSubEnvs == 1) ...
2467 
2468  if (useOnlyInter0Comm) {
2469  if (m_env.inter0Rank() >= 0) {
2470  queso_not_implemented();
2471  }
2472  else {
2473  // Node not in the 'inter0' communicator
2474  // Do nothing
2475  }
2476  }
2477  else {
2478  queso_error_msg("parallel vectors not supported yet");
2479  }
2480 
2481  //m_env.fullComm().Barrier();
2482 
2483  return resultValue;
2484 }
2485 // --------------------------------------------------
2486 template <class T>
2487 void
2488 ScalarSequence<T>::append(
2489  const ScalarSequence<T>& src,
2490  unsigned int srcInitialPos,
2491  unsigned int srcNumPos)
2492 {
2493  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+1), "srcInitialPos is too big");
2494 
2495  queso_require_greater_equal_msg(src.subSequenceSize(), (srcInitialPos+srcNumPos), "srcNumPos is too big");
2496 
2497  deleteStoredScalars();
2498  unsigned int currentSize = this->subSequenceSize();
2499  m_seq.resize(currentSize+srcNumPos,0.);
2500  for (unsigned int i = 0; i < srcNumPos; ++i) {
2501  m_seq[currentSize+i] = src.m_seq[srcInitialPos+i];
2502  }
2503 
2504  return;
2505 }
2506 // --------------------------------------------------
2507 template <class T>
2508 T
2509 ScalarSequence<T>::subPositionsOfMaximum(
2510  const ScalarSequence<T>& subCorrespondingScalarValues,
2511  ScalarSequence<T>& subPositionsOfMaximum)
2512 {
2513  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2514 
2515  T subMaxValue = subCorrespondingScalarValues.subMaxPlain();
2516  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2517 
2518  unsigned int subNumPos = 0;
2519  for (unsigned int i = 0; i < iMax; ++i) {
2520  if (subCorrespondingScalarValues[i] == subMaxValue) {
2521  subNumPos++;
2522  }
2523  }
2524 
2525  subPositionsOfMaximum.resizeSequence(subNumPos);
2526  unsigned int j = 0;
2527  for (unsigned int i = 0; i < iMax; ++i) {
2528  if (subCorrespondingScalarValues[i] == subMaxValue) {
2529  subPositionsOfMaximum[j] = (*this)[i];
2530  j++;
2531  }
2532  }
2533 
2534  return subMaxValue;
2535 }
2536 // --------------------------------------------------
2537 template <class T>
2538 T
2539 ScalarSequence<T>::unifiedPositionsOfMaximum( // rr0
2540  const ScalarSequence<T>& subCorrespondingScalarValues,
2541  ScalarSequence<T>& unifiedPositionsOfMaximum)
2542 {
2543  queso_require_equal_to_msg(subCorrespondingScalarValues.subSequenceSize(), this->subSequenceSize(), "invalid input");
2544 
2545  T maxValue = subCorrespondingScalarValues.subMaxPlain();
2546  unsigned int iMax = subCorrespondingScalarValues.subSequenceSize();
2547 
2548  unsigned int numPos = 0;
2549  for (unsigned int i = 0; i < iMax; ++i) {
2550  if (subCorrespondingScalarValues[i] == maxValue) {
2551  numPos++;
2552  }
2553  }
2554 
2555  unifiedPositionsOfMaximum.resizeSequence(numPos);
2556  unsigned int j = 0;
2557  for (unsigned int i = 0; i < iMax; ++i) {
2558  if (subCorrespondingScalarValues[i] == maxValue) {
2559  unifiedPositionsOfMaximum[j] = (*this)[i];
2560  j++;
2561  }
2562  }
2563 
2564  return maxValue;
2565 }
2566 // --------------------------------------------------
2567 template <class T>
2568 void
2569 ScalarSequence<T>::subWriteContents(
2570  unsigned int initialPos,
2571  unsigned int numPos,
2572  const std::string& fileName,
2573  const std::string& fileType,
2574  const std::set<unsigned int>& allowedSubEnvIds) const
2575 {
2576  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
2577 
2578  FilePtrSetStruct filePtrSet;
2579  if (m_env.openOutputFile(fileName,
2580  fileType,
2581  allowedSubEnvIds,
2582  false, // A 'true' causes problems when the user chooses (via options
2583  // in the input file) to use just one file for all outputs.
2584  filePtrSet)) {
2585 
2586  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT ||
2587  fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2588  this->subWriteContents(initialPos,
2589  numPos,
2590  *filePtrSet.ofsVar,
2591  fileType);
2592  }
2593 #ifdef QUESO_HAS_HDF5
2594  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2595 
2596  // Create dataspace
2597  hsize_t dims[1] = { this->subSequenceSize() };
2598  hid_t dataspace_id = H5Screate_simple(1, dims, dims);
2599  queso_require_greater_equal_msg(
2600  dataspace_id,
2601  0,
2602  "error create dataspace with id: " << dataspace_id);
2603 
2604  // Create dataset
2605  hid_t dataset_id = H5Dcreate(filePtrSet.h5Var,
2606  "data",
2607  H5T_IEEE_F64LE,
2608  dataspace_id,
2609  H5P_DEFAULT,
2610  H5P_DEFAULT,
2611  H5P_DEFAULT);
2612 
2613  queso_require_greater_equal_msg(
2614  dataset_id,
2615  0,
2616  "error creating dataset with id: " << dataset_id);
2617 
2618  // Write the dataset
2619  herr_t status = H5Dwrite(
2620  dataset_id,
2621  H5T_NATIVE_DOUBLE, // The type in memory
2622  H5S_ALL, // The dataspace in memory
2623  dataspace_id, // The file dataspace
2624  H5P_DEFAULT, // Xfer property list
2625  &m_seq[0]);
2626 
2627  queso_require_greater_equal_msg(
2628  status,
2629  0,
2630  "error writing dataset to file with id: " << filePtrSet.h5Var);
2631 
2632  // Clean up
2633  H5Dclose(dataset_id);
2634  H5Sclose(dataspace_id);
2635  }
2636 #endif
2637 
2638  m_env.closeFile(filePtrSet,fileType);
2639  }
2640 }
2641 // --------------------------------------------------
2642 template <class T>
2643 void
2644 ScalarSequence<T>::subWriteContents( // rr0
2645  unsigned int /* initialPos */,
2646  unsigned int /* numPos */,
2647  std::ofstream& ofs,
2648  const std::string& fileType) const
2649 {
2650  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2651  this->writeSubMatlabHeader(ofs, this->subSequenceSize());
2652  }
2653  else if (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
2654  this->writeTxtHeader(ofs, this->subSequenceSize());
2655  }
2656 
2657  unsigned int chainSize = this->subSequenceSize();
2658  for (unsigned int j = 0; j < chainSize; ++j) {
2659  ofs << m_seq[j]
2660  << std::endl;
2661  }
2662 
2663  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2664  ofs << "];\n";
2665  }
2666 
2667  return;
2668 }
2669 // --------------------------------------------------
2670 template <class T>
2671 void
2672 ScalarSequence<T>::unifiedWriteContents(
2673  const std::string& fileName,
2674  const std::string& inputFileType) const
2675 {
2676  std::string fileType(inputFileType);
2677 #ifdef QUESO_HAS_HDF5
2678  // Do nothing
2679 #else
2680  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2681  if (m_env.subDisplayFile()) {
2682  *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
2683  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2684  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2685  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2686  << "' instead..."
2687  << std::endl;
2688  }
2689  if (m_env.subRank() == 0) {
2690  std::cerr << "WARNING in ScalarSequence<T>::unifiedWriteContents()"
2691  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2692  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2693  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2694  << "' instead..."
2695  << std::endl;
2696  }
2697  fileType = UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT;
2698  }
2699 #endif
2700 
2701  // All processors in 'fullComm' should call this routine...
2702 
2703  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2704  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
2705  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedWriteContents()"
2706  << ": worldRank " << m_env.worldRank()
2707  << ", subEnvironment " << m_env.subId()
2708  << ", subRank " << m_env.subRank()
2709  << ", inter0Rank " << m_env.inter0Rank()
2710  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
2711  << ", fileName = " << fileName
2712  << ", fileType = " << fileType
2713  << std::endl;
2714  }
2715 
2716  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
2717 
2718  if (m_env.inter0Rank() >= 0) {
2719  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
2720  if (m_env.inter0Rank() == (int) r) {
2721  // My turn
2722  FilePtrSetStruct unifiedFilePtrSet;
2723  // bool writeOver = (r == 0);
2724  bool writeOver = false; // A 'true' causes problems when the user chooses (via options
2725  // in the input file) to use just one file for all outputs.
2726  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 000 \n" << std::endl;
2727  if (m_env.openUnifiedOutputFile(fileName,
2728  fileType, // "m or hdf"
2729  writeOver,
2730  unifiedFilePtrSet)) {
2731  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 001 \n" << std::endl;
2732 
2733  unsigned int chainSize = this->subSequenceSize();
2734  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2735  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2736 
2737  // Write header info
2738  if (r == 0) {
2739  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2740  this->writeUnifiedMatlabHeader(*unifiedFilePtrSet.ofsVar,
2741  this->subSequenceSize()*m_env.inter0Comm().NumProc());
2742  }
2743  else { // It's definitely txt if we get in here
2744  this->writeTxtHeader(*unifiedFilePtrSet.ofsVar,
2745  this->subSequenceSize()*m_env.inter0Comm().NumProc());
2746  }
2747  }
2748 
2749  for (unsigned int j = 0; j < chainSize; ++j) {
2750  *unifiedFilePtrSet.ofsVar << m_seq[j]
2751  << std::endl;
2752  }
2753 
2754  m_env.closeFile(unifiedFilePtrSet,fileType);
2755  }
2756 #ifdef QUESO_HAS_HDF5
2757  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2758  unsigned int numParams = 1; // m_vectorSpace.dimLocal();
2759  if (r == 0) {
2760  hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
2761  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type created" << std::endl;
2762  hsize_t dimsf[1];
2763  dimsf[0] = chainSize;
2764  hid_t dataspace = H5Screate_simple(1, dimsf, NULL); // HDF5_rank = 2
2765  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space created" << std::endl;
2766  hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var,
2767  "data",
2768  datatype,
2769  dataspace,
2770  H5P_DEFAULT, // Link creation property list
2771  H5P_DEFAULT, // Dataset creation property list
2772  H5P_DEFAULT); // Dataset access property list
2773  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set created" << std::endl;
2774 
2775  struct timeval timevalBegin;
2776  int iRC = UQ_OK_RC;
2777  iRC = gettimeofday(&timevalBegin,NULL);
2778  if (iRC) {}; // just to remove compiler warning
2779 
2780  herr_t status;
2781  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 002 \n" << std::endl;
2782  status = H5Dwrite(dataset,
2783  H5T_NATIVE_DOUBLE,
2784  H5S_ALL,
2785  H5S_ALL,
2786  H5P_DEFAULT,
2787  &m_seq[0]);
2788  if (status) {}; // just to remove compiler warning
2789 
2790  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 003 \n" << std::endl;
2791  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data written" << std::endl;
2792 
2793  double writeTime = MiscGetEllapsedSeconds(&timevalBegin);
2794  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2795  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedWriteContents()"
2796  << ": worldRank " << m_env.worldRank()
2797  << ", fullRank " << m_env.fullRank()
2798  << ", subEnvironment " << m_env.subId()
2799  << ", subRank " << m_env.subRank()
2800  << ", inter0Rank " << m_env.inter0Rank()
2801  << ", fileName = " << fileName
2802  << ", numParams = " << numParams
2803  << ", chainSize = " << chainSize
2804  << ", writeTime = " << writeTime << " seconds"
2805  << std::endl;
2806  }
2807 
2808  H5Dclose(dataset);
2809  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data set closed" << std::endl;
2810  H5Sclose(dataspace);
2811  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data space closed" << std::endl;
2812  H5Tclose(datatype);
2813  //std::cout << "In ScalarSequence<T>::unifiedWriteContents(): h5 case, data type closed" << std::endl;
2814  }
2815  else {
2816  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
2817  }
2818  }
2819 #endif
2820  } // if (m_env.openUnifiedOutputFile())
2821  //std::cout << "\n In ScalarSequence<T>::unifiedWriteContents(), pos 004 \n" << std::endl;
2822  } // if (m_env.inter0Rank() == (int) r)
2823  m_env.inter0Comm().Barrier();
2824  } // for r
2825 
2826  if (m_env.inter0Rank() == 0) {
2827  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2828  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2829  FilePtrSetStruct unifiedFilePtrSet;
2830  if (m_env.openUnifiedOutputFile(fileName,
2831  fileType,
2832  false, // Yes, 'writeOver = false' in order to close the array for matlab
2833  unifiedFilePtrSet)) {
2834 
2835  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
2836  *unifiedFilePtrSet.ofsVar << "];\n";
2837  }
2838 
2839  m_env.closeFile(unifiedFilePtrSet,fileType);
2840  }
2841  }
2842  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2843  // Do nothing
2844  }
2845  else {
2846  queso_error_msg("invalid file type");
2847  }
2848  }
2849  } // if (m_env.inter0Rank() >= 0)
2850 
2851  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
2852  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedWriteContents()"
2853  << ", fileName = " << fileName
2854  << ", fileType = " << fileType
2855  << std::endl;
2856  }
2857  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2858 
2859  return;
2860 }
2861 
2862 template <class T>
2863 void
2864 ScalarSequence<T>::writeUnifiedMatlabHeader(std::ofstream & ofs,
2865  double sequenceSize) const
2866 {
2867  ofs << m_name << "_unified" << " = zeros(" << sequenceSize
2868  << "," << 1
2869  << ");"
2870  << std::endl;
2871  ofs << m_name << "_unified" << " = [";
2872 }
2873 
2874 template <class T>
2875 void
2876 ScalarSequence<T>::writeSubMatlabHeader(std::ofstream & ofs,
2877  double sequenceSize) const
2878 {
2879  ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << sequenceSize
2880  << "," << 1
2881  << ");"
2882  << std::endl;
2883  ofs << m_name << "_sub" << m_env.subIdString() << " = [";
2884 }
2885 
2886 template <class T>
2887 void
2888 ScalarSequence<T>::writeTxtHeader(std::ofstream & ofs,
2889  double sequenceSize) const
2890 {
2891  ofs << sequenceSize << " " << 1 << std::endl;
2892 }
2893 
2894 // --------------------------------------------------
2895 template <class T>
2896 void
2897 ScalarSequence<T>::unifiedReadContents(
2898  const std::string& fileName,
2899  const std::string& inputFileType,
2900  const unsigned int subReadSize)
2901 {
2902  queso_require_not_equal_to_msg(inputFileType, std::string(UQ_FILE_EXTENSION_FOR_TXT_FORMAT), std::string("reading txt files is not yet supported"));
2903  std::string fileType(inputFileType);
2904 #ifdef QUESO_HAS_HDF5
2905  // Do nothing
2906 #else
2907  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
2908  if (m_env.subDisplayFile()) {
2909  *m_env.subDisplayFile() << "WARNING in ScalarSequence<T>::unifiedReadContents()"
2910  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2911  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2912  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2913  << "' instead..."
2914  << std::endl;
2915  }
2916  if (m_env.subRank() == 0) {
2917  std::cerr << "WARNING in ScalarSequence<T>::unifiedReadContents()"
2918  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2919  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
2920  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
2921  << "' instead..."
2922  << std::endl;
2923  }
2924  fileType = UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT;
2925  }
2926 #endif
2927 
2928  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2929  if (m_env.subDisplayFile()) {
2930  *m_env.subDisplayFile() << "Entering ScalarSequence<T>::unifiedReadContents()"
2931  << ": worldRank " << m_env.worldRank()
2932  << ", fullRank " << m_env.fullRank()
2933  << ", subEnvironment " << m_env.subId()
2934  << ", subRank " << m_env.subRank()
2935  << ", inter0Rank " << m_env.inter0Rank()
2936  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
2937  << ", fileName = " << fileName
2938  << ", subReadSize = " << subReadSize
2939  //<< ", unifiedReadSize = " << unifiedReadSize
2940  << std::endl;
2941  }
2942 
2943  this->resizeSequence(subReadSize);
2944 
2945  if (m_env.inter0Rank() >= 0) {
2946  double unifiedReadSize = subReadSize*m_env.inter0Comm().NumProc();
2947 
2948  // In the logic below, the id of a line' begins with value 0 (zero)
2949  unsigned int idOfMyFirstLine = 1 + m_env.inter0Rank()*subReadSize;
2950  unsigned int idOfMyLastLine = (1 + m_env.inter0Rank())*subReadSize;
2951  unsigned int numParams = 1; // this->vectorSizeLocal();
2952 
2953  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // "m or hdf"
2954  if (m_env.inter0Rank() == (int) r) {
2955  // My turn
2956  FilePtrSetStruct unifiedFilePtrSet;
2957  if (m_env.openUnifiedInputFile(fileName,
2958  fileType,
2959  unifiedFilePtrSet)) {
2960  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
2961  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
2962  if (r == 0) {
2963  // Read number of chain positions in the file by taking care of the first line,
2964  // which resembles something like 'variable_name = zeros(n_positions,m_params);'
2965  std::string tmpString;
2966 
2967  // Read 'variable name' string
2968  *unifiedFilePtrSet.ifsVar >> tmpString;
2969  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2970 
2971  // Read '=' sign
2972  *unifiedFilePtrSet.ifsVar >> tmpString;
2973  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2974  queso_require_equal_to_msg(tmpString, std::string("="), std::string("string should be the '=' sign"));
2975 
2976  // Read 'zeros(n_positions,n_params)' string
2977  *unifiedFilePtrSet.ifsVar >> tmpString;
2978  //std::cout << "Just read '" << tmpString << "'" << std::endl;
2979  unsigned int posInTmpString = 6;
2980 
2981  // Isolate 'n_positions' in a string
2982  //char nPositionsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
2983  std::string nPositionsString((size_t) (tmpString.size()-posInTmpString+1),' ');
2984  unsigned int posInPositionsString = 0;
2985  do {
2986  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
2987  nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
2988  } while (tmpString[posInTmpString] != ',');
2989  nPositionsString[posInPositionsString] = '\0';
2990 
2991  // Isolate 'n_params' in a string
2992  posInTmpString++; // Avoid reading ',' char
2993  //char nParamsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
2994  std::string nParamsString((size_t) (tmpString.size()-posInTmpString+1),' ');
2995  unsigned int posInParamsString = 0;
2996  do {
2997  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
2998  nParamsString[posInParamsString++] = tmpString[posInTmpString++];
2999  } while (tmpString[posInTmpString] != ')');
3000  nParamsString[posInParamsString] = '\0';
3001 
3002  // Convert 'n_positions' and 'n_params' strings to numbers
3003  unsigned int sizeOfChainInFile = (unsigned int) strtod(nPositionsString.c_str(),NULL);
3004  unsigned int numParamsInFile = (unsigned int) strtod(nParamsString.c_str(), NULL);
3005  if (m_env.subDisplayFile()) {
3006  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
3007  << ": worldRank " << m_env.worldRank()
3008  << ", fullRank " << m_env.fullRank()
3009  << ", sizeOfChainInFile = " << sizeOfChainInFile
3010  << ", numParamsInFile = " << numParamsInFile
3011  << std::endl;
3012  }
3013 
3014  // Check if [size of chain in file] >= [requested unified sequence size]
3015  queso_require_greater_equal_msg(sizeOfChainInFile, unifiedReadSize, "size of chain in file is not big enough");
3016 
3017  // Check if [num params in file] == [num params in current chain]
3018  queso_require_equal_to_msg(numParamsInFile, numParams, "number of parameters of chain in file is different than number of parameters in this chain object");
3019  } // if (r == 0)
3020 
3021  // Code common to any core in 'inter0Comm', including core of rank 0
3022  unsigned int maxCharsPerLine = 64*numParams; // Up to about 60 characters to represent each parameter value
3023 
3024  unsigned int lineId = 0;
3025  while (lineId < idOfMyFirstLine) {
3026  unifiedFilePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
3027  lineId++;
3028  };
3029 
3030  if (r == 0) {
3031  // Take care of initial part of the first data line,
3032  // which resembles something like 'variable_name = [value1 value2 ...'
3033  std::string tmpString;
3034 
3035  // Read 'variable name' string
3036  *unifiedFilePtrSet.ifsVar >> tmpString;
3037  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
3038 
3039  // Read '=' sign
3040  *unifiedFilePtrSet.ifsVar >> tmpString;
3041  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
3042  queso_require_equal_to_msg(tmpString, std::string("="), std::string("in core 0, string should be the '=' sign"));
3043 
3044  // Take into account the ' [' portion
3045  std::streampos tmpPos = unifiedFilePtrSet.ifsVar->tellg();
3046  unifiedFilePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
3047  }
3048 
3049  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3050  while (lineId <= idOfMyLastLine) {
3051  for (unsigned int i = 0; i < numParams; ++i) {
3052  *unifiedFilePtrSet.ifsVar >> tmpScalar;
3053  }
3054  m_seq[lineId - idOfMyFirstLine] = tmpScalar;
3055  lineId++;
3056  };
3057  }
3058 #ifdef QUESO_HAS_HDF5
3059  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
3060  if (r == 0) {
3061  hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var,
3062  "data",
3063  H5P_DEFAULT); // Dataset access property list
3064  hid_t datatype = H5Dget_type(dataset);
3065  H5T_class_t t_class = H5Tget_class(datatype);
3066  queso_require_equal_to_msg(t_class, H5T_FLOAT, "t_class is not H5T_DOUBLE");
3067  hid_t dataspace = H5Dget_space(dataset);
3068  int rank = H5Sget_simple_extent_ndims(dataspace);
3069  queso_require_equal_to_msg(rank, 2, "hdf rank is not 2");
3070  hsize_t dims_in[2];
3071  int status_n;
3072  status_n = H5Sget_simple_extent_dims(dataspace, dims_in, NULL);
3073  if (status_n) {}; // just to remove compiler warning
3074  //std::cout << "In ScalarSequence<T>::unifiedReadContents()"
3075  // << ": dims_in[0] = " << dims_in[0]
3076  // << ", dims_in[1] = " << dims_in[1]
3077  // << std::endl;
3078  queso_require_equal_to_msg(dims_in[0], numParams, "dims_in[0] is not equal to 'numParams'");
3079  queso_require_greater_equal_msg(dims_in[1], subReadSize, "dims_in[1] is smaller that requested 'subReadSize'");
3080 
3081  struct timeval timevalBegin;
3082  int iRC = UQ_OK_RC;
3083  iRC = gettimeofday(&timevalBegin,NULL);
3084  if (iRC) {}; // just to remove compiler warning
3085 
3086  unsigned int chainSizeIn = (unsigned int) dims_in[1];
3087  //double* dataIn[numParams]; // avoid compiler warning
3088  std::vector<double*> dataIn((size_t) numParams,NULL);
3089  dataIn[0] = (double*) malloc(numParams*chainSizeIn*sizeof(double));
3090  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
3091  dataIn[i] = dataIn[i-1] + chainSizeIn; // Yes, just 'chainSizeIn', not 'chainSizeIn*sizeof(double)'
3092  }
3093  //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, memory allocated" << std::endl;
3094  herr_t status;
3095  status = H5Dread(dataset,
3096  H5T_NATIVE_DOUBLE,
3097  H5S_ALL,
3098  dataspace,
3099  H5P_DEFAULT,
3100  dataIn[0]);
3101  if (status) {}; // just to remove compiler warning
3102  //std::cout << "In ScalarSequence<T>::unifiedReadContents(): h5 case, data read" << std::endl;
3103  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3104  for (unsigned int j = 0; j < subReadSize; ++j) { // Yes, 'subReadSize', not 'chainSizeIn'
3105  for (unsigned int i = 0; i < numParams; ++i) {
3106  tmpScalar = dataIn[i][j];
3107  }
3108  m_seq[j] = tmpScalar;
3109  }
3110 
3111  double readTime = MiscGetEllapsedSeconds(&timevalBegin);
3112  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
3113  *m_env.subDisplayFile() << "In ScalarSequence<T>::unifiedReadContents()"
3114  << ": worldRank " << m_env.worldRank()
3115  << ", fullRank " << m_env.fullRank()
3116  << ", subEnvironment " << m_env.subId()
3117  << ", subRank " << m_env.subRank()
3118  << ", inter0Rank " << m_env.inter0Rank()
3119  << ", fileName = " << fileName
3120  << ", numParams = " << numParams
3121  << ", chainSizeIn = " << chainSizeIn
3122  << ", subReadSize = " << subReadSize
3123  << ", readTime = " << readTime << " seconds"
3124  << std::endl;
3125  }
3126 
3127  H5Sclose(dataspace);
3128  H5Tclose(datatype);
3129  H5Dclose(dataset);
3130  //free(dataIn[0]); // related to the changes above for compiler warning
3131  for (unsigned int tmpIndex = 0; tmpIndex < dataIn.size(); tmpIndex++) {
3132  free (dataIn[tmpIndex]);
3133  }
3134  }
3135  else {
3136  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
3137  }
3138  }
3139 #endif
3140  else {
3141  queso_error_msg("invalid file type");
3142  }
3143  m_env.closeFile(unifiedFilePtrSet,fileType);
3144  } // if (m_env.openUnifiedInputFile())
3145  } // if (m_env.inter0Rank() == (int) r)
3146  m_env.inter0Comm().Barrier();
3147  } // for r
3148  } // if (m_env.inter0Rank() >= 0)
3149  else {
3150  T tmpScalar(0.); // V tmpVec(m_vectorSpace.zeroVector());
3151  for (unsigned int i = 1; i < subReadSize; ++i) {
3152  m_seq[i] = tmpScalar;
3153  }
3154  }
3155 
3156  if (m_env.subDisplayFile()) {
3157  *m_env.subDisplayFile() << "Leaving ScalarSequence<T>::unifiedReadContents()"
3158  << ", fileName = " << fileName
3159  << std::endl;
3160  }
3161  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
3162 
3163  return;
3164 }
3165 // Private methods-----------------------------------
3166 template <class T>
3167 void
3168 ScalarSequence<T>::copy(const ScalarSequence<T>& src)
3169 {
3170  m_name = src.m_name;
3171  m_seq.clear();
3172  m_seq.resize(src.subSequenceSize(),0.);
3173  for (unsigned int i = 0; i < m_seq.size(); ++i) {
3174  m_seq[i] = src.m_seq[i];
3175  }
3176  deleteStoredScalars();
3177 
3178  return;
3179 }
3180 // --------------------------------------------------
3181 template <class T>
3182 void
3183 ScalarSequence<T>::extractScalarSeq(
3184  unsigned int initialPos,
3185  unsigned int spacing,
3186  unsigned int numPos,
3187  ScalarSequence<T>& scalarSeq) const
3188 {
3189  scalarSeq.resizeSequence(numPos);
3190  if (spacing == 1) {
3191  for (unsigned int j = 0; j < numPos; ++j) {
3192  //if ((initialPos+j*spacing) > m_seq.size()) {
3193  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3194  // << ": initialPos = " << initialPos
3195  // << ", spacing = " << spacing
3196  // << ", numPos = " << numPos
3197  // << ", j = " << j
3198  // << ", position got too large"
3199  // << std::endl;
3200  //}
3201  scalarSeq[j] = m_seq[initialPos+j ];
3202  }
3203  }
3204  else {
3205  for (unsigned int j = 0; j < numPos; ++j) {
3206  //if ((initialPos+j*spacing) > m_seq.size()) {
3207  // std::cerr << "In ScalarSequence<T>::extraScalarSeq()"
3208  // << ": initialPos = " << initialPos
3209  // << ", spacing = " << spacing
3210  // << ", numPos = " << numPos
3211  // << ", j = " << j
3212  // << ", position got too large"
3213  // << std::endl;
3214  //}
3215  scalarSeq[j] = m_seq[initialPos+j*spacing];
3216  }
3217  }
3218 
3219  return;
3220 }
3221 // --------------------------------------------------
3222 template <class T>
3223 void
3224 ScalarSequence<T>::extractRawData(
3225  unsigned int initialPos,
3226  unsigned int spacing,
3227  unsigned int numPos,
3228  std::vector<double>& rawDataVec) const
3229 {
3230  rawDataVec.resize(numPos);
3231  if (spacing == 1) {
3232  for (unsigned int j = 0; j < numPos; ++j) {
3233  rawDataVec[j] = m_seq[initialPos+j ];
3234  }
3235  }
3236  else {
3237  for (unsigned int j = 0; j < numPos; ++j) {
3238  rawDataVec[j] = m_seq[initialPos+j*spacing];
3239  }
3240  }
3241 
3242  return;
3243 }
3244 // --------------------------------------------------
3245 template <class T>
3246 std::vector<T>&
3247 ScalarSequence<T>::rawData()
3248 {
3249  return m_seq;
3250 }
3251 
3252 // --------------------------------------------------
3253 template <class T>
3254 void
3255 ScalarSequence<T>::subSort()
3256 {
3257  std::sort(m_seq.begin(), m_seq.end());
3258  return;
3259 }
3260 // --------------------------------------------------
3261 // Acknowledgement: the tree idea in this routine was taken from
3262 // 'http://penguin.ewu.edu/~trolfe/ParallelMerge/ParallelMerge.html', as of March 08, 2009
3263 template <class T>
3264 void
3265 ScalarSequence<T>::parallelMerge(
3266  std::vector<T>& sortedBuffer,
3267  const std::vector<T>& leafData,
3268  unsigned int currentTreeLevel) const
3269 {
3270  int parentNode = m_env.inter0Rank() & ~(1 << currentTreeLevel);
3271 
3272  if (m_env.inter0Rank() >= 0) { // KAUST
3273  if (currentTreeLevel == 0) {
3274  // Leaf node: sort own local data.
3275  unsigned int leafDataSize = leafData.size();
3276  sortedBuffer.resize(leafDataSize,0.);
3277  for (unsigned int i = 0; i < leafDataSize; ++i) {
3278  sortedBuffer[i] = leafData[i];
3279  }
3280  std::sort(sortedBuffer.begin(), sortedBuffer.end());
3281  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3282  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3283  << ": tree node " << m_env.inter0Rank()
3284  << ", leaf sortedBuffer[0] = " << sortedBuffer[0]
3285  << ", leaf sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3286  << std::endl;
3287  }
3288  }
3289  else {
3290  int nextTreeLevel = currentTreeLevel - 1;
3291  int rightChildNode = m_env.inter0Rank() | (1 << nextTreeLevel);
3292 
3293  if (rightChildNode >= m_env.inter0Comm().NumProc()) { // No right child. Move down one level.
3294  this->parallelMerge(sortedBuffer,
3295  leafData,
3296  nextTreeLevel);
3297  }
3298  else {
3299  unsigned int uintBuffer[1];
3300  uintBuffer[0] = nextTreeLevel;
3301  m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_INIT_MPI_MSG,
3302  "ScalarSequence<T>::parallelMerge()",
3303  "failed MPI.Send() for init");
3304 
3305  this->parallelMerge(sortedBuffer,
3306  leafData,
3307  nextTreeLevel);
3308 
3309  // Prepare variable 'leftSortedBuffer': just copy own current sorted data.
3310  unsigned int leftSize = sortedBuffer.size();
3311  std::vector<T> leftSortedBuffer(leftSize,0.);
3312  for (unsigned int i = 0; i < leftSize; ++i) {
3313  leftSortedBuffer[i] = sortedBuffer[i];
3314  }
3315 
3316  // Prepare variable 'rightSortedBuffer': receive data from right child node.
3317  RawType_MPI_Status status;
3318  m_env.inter0Comm().Recv((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, rightChildNode, SCALAR_SEQUENCE_SIZE_MPI_MSG, &status,
3319  "ScalarSequence<T>::parallelMerge()",
3320  "failed MPI.Recv() for size");
3321  //if (status) {}; // just to remove compiler warning
3322 
3323  unsigned int rightSize = uintBuffer[0];
3324  std::vector<T> rightSortedBuffer(rightSize,0.);
3325  m_env.inter0Comm().Recv((void *) &rightSortedBuffer[0], (int) rightSize, RawValue_MPI_DOUBLE, rightChildNode, SCALAR_SEQUENCE_DATA_MPI_MSG, &status,
3326  "ScalarSequence<T>::parallelMerge()",
3327  "failed MPI.Recv() for data");
3328 
3329  // Merge the two results into 'sortedBuffer'.
3330  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3331  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3332  << ": tree node " << m_env.inter0Rank()
3333  << " is combining " << leftSortedBuffer.size()
3334  << " left doubles with " << rightSortedBuffer.size()
3335  << " right doubles"
3336  << std::endl;
3337  }
3338 
3339  sortedBuffer.clear();
3340  sortedBuffer.resize(leftSortedBuffer.size()+rightSortedBuffer.size(),0.);
3341  unsigned int i = 0;
3342  unsigned int j = 0;
3343  unsigned int k = 0;
3344  while ((i < leftSize ) &&
3345  (j < rightSize)) {
3346  if (leftSortedBuffer[i] > rightSortedBuffer[j]) sortedBuffer[k++] = rightSortedBuffer[j++];
3347  else sortedBuffer[k++] = leftSortedBuffer [i++];
3348  }
3349  while (i < leftSize ) sortedBuffer[k++] = leftSortedBuffer [i++];
3350  while (j < rightSize) sortedBuffer[k++] = rightSortedBuffer[j++];
3351  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 0)) {
3352  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3353  << ": tree node " << m_env.inter0Rank()
3354  << ", merged sortedBuffer[0] = " << sortedBuffer[0]
3355  << ", merged sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3356  << std::endl;
3357  }
3358  }
3359  }
3360 
3361  if (parentNode != m_env.inter0Rank()) {
3362  // Transmit data to parent node.
3363  unsigned int uintBuffer[1];
3364  uintBuffer[0] = sortedBuffer.size();
3365  m_env.inter0Comm().Send((void *) uintBuffer, 1, RawValue_MPI_UNSIGNED, parentNode, SCALAR_SEQUENCE_SIZE_MPI_MSG,
3366  "ScalarSequence<T>::parallelMerge()",
3367  "failed MPI.Send() for size");
3368 
3369  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
3370  *m_env.subDisplayFile() << "In ScalarSequence<T>::parallelMerge()"
3371  << ": tree node " << m_env.inter0Rank()
3372  << " is sending " << sortedBuffer.size()
3373  << " doubles to tree node " << parentNode
3374  << ", with sortedBuffer[0] = " << sortedBuffer[0]
3375  << " and sortedBuffer[" << sortedBuffer.size()-1 << "] = " << sortedBuffer[sortedBuffer.size()-1]
3376  << std::endl;
3377  }
3378 
3379  m_env.inter0Comm().Send((void *) &sortedBuffer[0], (int) sortedBuffer.size(), RawValue_MPI_DOUBLE, parentNode, SCALAR_SEQUENCE_DATA_MPI_MSG,
3380  "ScalarSequence<T>::parallelMerge()",
3381  "failed MPI.Send() for data");
3382  }
3383  } // KAUST
3384 
3385  return;
3386 }
3387 
3388 // --------------------------------------------------
3389 // Methods conditionally available ------------------
3390 // --------------------------------------------------
3391 // --------------------------------------------------
3392 
3393 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
3394 template <class T>
3395 void
3396 ScalarSequence<T>::subUniformlySampledMdf(
3397  unsigned int numEvaluationPoints,
3398  T& minDomainValue,
3399  T& maxDomainValue,
3400  std::vector<T>& mdfValues) const
3401 {
3402  T tmpMinValue;
3403  T tmpMaxValue;
3404  std::vector<T> centers(numEvaluationPoints,0.);
3405  std::vector<unsigned int> bins (numEvaluationPoints,0);
3406 
3407  subMinMaxExtra(0, // initialPos
3408  this->subSequenceSize(),
3409  tmpMinValue,
3410  tmpMaxValue);
3411  subHistogram(0, // initialPos,
3412  tmpMinValue,
3413  tmpMaxValue,
3414  centers,
3415  bins);
3416 
3417  minDomainValue = centers[0];
3418  maxDomainValue = centers[centers.size()-1];
3419  T binSize = (maxDomainValue - minDomainValue)/((double)(centers.size() - 1));
3420 
3421  unsigned int totalCount = 0;
3422  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
3423  totalCount += bins[i];
3424  }
3425 
3426  mdfValues.clear();
3427  mdfValues.resize(numEvaluationPoints);
3428  for (unsigned int i = 0; i < numEvaluationPoints; ++i) {
3429  mdfValues[i] = ((T) bins[i])/((T) totalCount)/binSize;
3430  }
3431 
3432  return;
3433 }
3434 // --------------------------------------------------
3435 template <class T>
3436 T
3437 ScalarSequence<T>::subMeanCltStd(
3438  unsigned int initialPos,
3439  unsigned int numPos,
3440  const T& meanValue) const
3441 {
3442  if (this->subSequenceSize() == 0) return 0.;
3443 
3444  bool bRC = ((initialPos < this->subSequenceSize()) &&
3445  (0 < numPos ) &&
3446  ((initialPos+numPos) <= this->subSequenceSize()));
3447  queso_require_msg(bRC, "invalid input data");
3448 
3449  unsigned int finalPosPlus1 = initialPos + numPos;
3450  T diff;
3451  T stdValue = 0.;
3452  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
3453  diff = m_seq[j] - meanValue;
3454  stdValue += diff*diff;
3455  }
3456 
3457  stdValue /= (((T) numPos) - 1.);
3458  stdValue /= (((T) numPos) - 1.);
3459  stdValue = sqrt(stdValue);
3460 
3461  return stdValue;
3462 }
3463 // --------------------------------------------------
3464 template <class T>
3465 T
3466 ScalarSequence<T>::unifiedMeanCltStd(
3467  bool useOnlyInter0Comm,
3468  unsigned int initialPos,
3469  unsigned int numPos,
3470  const T& unifiedMeanValue) const
3471 {
3472  if (m_env.numSubEnvironments() == 1) {
3473  return this->subMeanCltStd(initialPos,
3474  numPos,
3475  unifiedMeanValue);
3476  }
3477 
3478  // As of 14/Nov/2009, this routine does *not* require sub sequences to have equal size. Good.
3479 
3480  T unifiedStdValue = 0.;
3481  if (useOnlyInter0Comm) {
3482  if (m_env.inter0Rank() >= 0) {
3483  bool bRC = ((initialPos < this->subSequenceSize()) &&
3484  (0 < numPos ) &&
3485  ((initialPos+numPos) <= this->subSequenceSize()));
3486  queso_require_msg(bRC, "invalid input data");
3487 
3488  unsigned int finalPosPlus1 = initialPos + numPos;
3489  T diff;
3490  T localStdValue = 0.;
3491  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
3492  diff = m_seq[j] - unifiedMeanValue;
3493  localStdValue += diff*diff;
3494  }
3495 
3496  unsigned int unifiedNumPos = 0;
3497  m_env.inter0Comm().template Allreduce<unsigned int>(&numPos, &unifiedNumPos, (int) 1, RawValue_MPI_SUM,
3498  "ScalarSequence<T>::unifiedMeanCltStd()",
3499  "failed MPI.Allreduce() for numPos");
3500 
3501  m_env.inter0Comm().template Allreduce<double>(&localStdValue, &unifiedStdValue, (int) 1, RawValue_MPI_SUM,
3502  "ScalarSequence<T>::unifiedMeanCltStd()",
3503  "failed MPI.Allreduce() for stdValue");
3504 
3505  unifiedStdValue /= (((T) unifiedNumPos) - 1.);
3506  unifiedStdValue /= (((T) unifiedNumPos) - 1.);
3507  unifiedStdValue /= sqrt(unifiedStdValue);
3508  }
3509  else {
3510  // Node not in the 'inter0' communicator
3511  this->subMeanCltStd(initialPos,
3512  numPos,
3513  unifiedMeanValue);
3514  }
3515  }
3516  else {
3517  queso_error_msg("parallel vectors not supported yet");
3518  }
3519  //m_env.fullComm().Barrier();
3520 
3521  return unifiedStdValue;
3522 }
3523 //---------------------------------------------------
3524 template <class T>
3525 T
3526 ScalarSequence<T>::bmm(
3527  unsigned int initialPos,
3528  unsigned int batchLength) const
3529 {
3530  bool bRC = ((initialPos < this->subSequenceSize() ) &&
3531  (batchLength < (this->subSequenceSize()-initialPos)));
3532  queso_require_msg(bRC, "invalid input data");
3533 
3534  unsigned int numberOfBatches = (this->subSequenceSize() - initialPos)/batchLength;
3535  ScalarSequence<T> batchMeans(m_env,numberOfBatches,"");
3536 
3537  for (unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
3538  batchMeans[batchId] = this->subMeanExtra(initialPos + batchId*batchLength,
3539  batchLength);
3540  }
3541 
3542  T meanOfBatchMeans = batchMeans.subMeanExtra(0,
3543  batchMeans.subSequenceSize());
3544 
3545  //T covLag0OfBatchMeans = batchMeans.autoCovariance(0,
3546  // batchMeans.subSequenceSize(),
3547  // meanOfBatchMeans,
3548  // 0); // lag
3549 
3550  T bmmValue = batchMeans.subSampleVarianceExtra(0,
3551  batchMeans.subSequenceSize(),
3552  meanOfBatchMeans);
3553 
3554  bmmValue /= (T) batchMeans.subSequenceSize(); // CHECK
3555 //bmmValue *= (T) (this->subSequenceSize() - initialPos); // CHECK
3556 
3557  return bmmValue;
3558 }
3559 //---------------------------------------------------
3560 template <class T>
3561 void
3562 ScalarSequence<T>::psd(
3563  unsigned int initialPos,
3564  unsigned int numBlocks,
3565  double hopSizeRatio,
3566  std::vector<double>& psdResult) const
3567 {
3568  bool bRC = ((initialPos < this->subSequenceSize() ) &&
3569  (hopSizeRatio != 0. ) &&
3570  (numBlocks < (this->subSequenceSize() - initialPos)) &&
3571  (fabs(hopSizeRatio) < (double) (this->subSequenceSize() - initialPos)));
3572  queso_require_msg(bRC, "invalid input data");
3573 
3574  unsigned int dataSize = this->subSequenceSize() - initialPos;
3575 
3576  T meanValue = this->subMeanExtra(initialPos,
3577  dataSize);
3578 
3579  // Determine hopSize and blockSize
3580  unsigned int hopSize = 0;
3581  unsigned int blockSize = 0;
3582  if (hopSizeRatio <= -1.) {
3583  double overlapSize = -hopSizeRatio;
3584  double tmp = ((double) (dataSize + (numBlocks - 1)*overlapSize))/((double) numBlocks);
3585  blockSize = (unsigned int) tmp;
3586  }
3587  else if (hopSizeRatio < 0.) {
3588  double tmp = ((double) dataSize)/(( ((double) numBlocks) - 1. )*(-hopSizeRatio) + 1.);
3589  blockSize = (unsigned int) tmp;
3590  hopSize = (unsigned int) ( ((double) blockSize) * (-hopSizeRatio) );
3591  }
3592  else if (hopSizeRatio == 0.) {
3593  queso_error_msg("hopSizeRatio == 0");
3594  }
3595  else if (hopSizeRatio < 1.) {
3596  double tmp = ((double) dataSize)/(( ((double) numBlocks) - 1. )*hopSizeRatio + 1.);
3597  blockSize = (unsigned int) tmp;
3598  hopSize = (unsigned int) ( ((double) blockSize) * hopSizeRatio );
3599  }
3600  else {
3601  hopSize = (unsigned int) hopSizeRatio;
3602  blockSize = dataSize - (numBlocks - 1)*hopSize;
3603  }
3604 
3605  int numberOfDiscardedDataElements = dataSize - (numBlocks-1)*hopSize - blockSize;
3606  //if (m_env.subDisplayFile()) {
3607  // *m_env.subDisplayFile() << "initialPos = " << initialPos
3608  // << ", N = " << dataSize
3609  // << ", #Blocks = " << numBlocks
3610  // << ", R (hop size) = " << hopSize
3611  // << ", B (block size) = " << blockSize
3612  // << ", overlap = " << blockSize - hopSize
3613  // << ", [(#Blocks - 1) * R + B] = " << (numBlocks-1)*hopSize + blockSize
3614  // << ", numberOfDiscardedDataElements = " << numberOfDiscardedDataElements
3615  // << std::endl;
3616  //}
3617  queso_require_greater_equal_msg(numberOfDiscardedDataElements, 0., "eventual extra space for last block should not be negative");
3618 
3619  double tmp = log((double) blockSize)/log(2.);
3620  double fractionalPart = tmp - ((double) ((unsigned int) tmp));
3621  if (fractionalPart > 0.) tmp += (1. - fractionalPart);
3622  unsigned int fftSize = (unsigned int) std::pow(2.,tmp);
3623  //if (m_env.subDisplayFile()) {
3624  // *m_env.subDisplayFile() << "fractionalPart = " << fractionalPart
3625  // << ", B = " << blockSize
3626  // << ", fftSize = " << fftSize
3627  // << std::endl;
3628  //}
3629 
3630  double modificationScale = 0.;
3631  for (unsigned int j = 0; j < blockSize; ++j) {
3632  double tmpValue = MiscHammingWindow(blockSize-1,j);
3633  modificationScale += tmpValue*tmpValue;
3634  }
3635  modificationScale = 1./modificationScale;
3636 
3637  std::vector<double> blockData(blockSize,0.);
3638  Fft<T> fftObj(m_env);
3639  std::vector<std::complex<double> > fftResult(0,std::complex<double>(0.,0.));
3640 
3641 #if 0
3642  unsigned int halfFFTSize = fftSize/2;
3643  psdResult.clear();
3644  psdResult.resize(1+halfFFTSize,0.);
3645  double factor = 1./M_PI/((double) numBlocks); // /((double) blockSize);
3646 #else
3647  psdResult.clear();
3648  psdResult.resize(fftSize,0.);
3649  double factor = 1./2./M_PI/((double) numBlocks); // /((double) blockSize);
3650 #endif
3651 
3652  for (unsigned int blockId = 0; blockId < numBlocks; blockId++) {
3653  // Fill block using Hamming window
3654  unsigned int initialDataPos = initialPos + blockId*hopSize;
3655  for (unsigned int j = 0; j < blockSize; ++j) {
3656  unsigned int dataPos = initialDataPos + j;
3657  queso_require_less_msg(dataPos, dataSize, "too large position to be accessed in data");
3658  blockData[j] = MiscHammingWindow(blockSize-1,j) * ( m_seq[dataPos] - meanValue ); // IMPORTANT
3659  }
3660 
3661  fftObj.forward(blockData,fftSize,fftResult);
3662 
3663  //if (m_env.subDisplayFile()) {
3664  // *m_env.subDisplayFile() << "blockData.size() = " << blockData.size()
3665  // << ", fftSize = " << fftSize
3666  // << ", fftResult.size() = " << fftResult.size()
3667  // << ", psdResult.size() = " << psdResult.size()
3668  // << std::endl;
3669  //}
3670  // Normalized spectral density: power per radians per sample
3671  for (unsigned int j = 0; j < psdResult.size(); ++j) {
3672  psdResult[j] += norm(fftResult[j])*modificationScale;
3673  }
3674  }
3675 
3676  for (unsigned int j = 0; j < psdResult.size(); ++j) {
3677  psdResult[j] *= factor;
3678  }
3679 
3680  return;
3681 }
3682 //---------------------------------------------------
3683 template <class T>
3684 T
3685 ScalarSequence<T>::geweke(
3686  unsigned int initialPos,
3687  double ratioNa,
3688  double ratioNb) const
3689 {
3690  double doubleFullDataSize = (double) (this->subSequenceSize() - initialPos);
3691  ScalarSequence<T> tmpSeq(m_env,0,"");
3692  std::vector<double> psdResult(0,0.);
3693 
3694  unsigned int dataSizeA = (unsigned int) (doubleFullDataSize * ratioNa);
3695  double doubleDataSizeA = (double) dataSizeA;
3696  unsigned int initialPosA = initialPos;
3697  this->extractScalarSeq(initialPosA,
3698  1,
3699  dataSizeA,
3700  tmpSeq);
3701  double meanA = tmpSeq.subMeanExtra(0,
3702  dataSizeA);
3703  tmpSeq.psd(0,
3704  (unsigned int) std::sqrt((double) dataSizeA), // numBlocks
3705  .5, // hopSizeRatio
3706  psdResult);
3707  double psdA = psdResult[0];
3708  double varOfMeanA = 2.*M_PI*psdA/doubleDataSizeA;
3709 
3710  unsigned int dataSizeB = (unsigned int) (doubleFullDataSize * ratioNb);
3711  double doubleDataSizeB = (double) dataSizeB;
3712  unsigned int initialPosB = this->subSequenceSize() - dataSizeB;
3713  this->extractScalarSeq(initialPosB,
3714  1,
3715  dataSizeB,
3716  tmpSeq);
3717  double meanB = tmpSeq.subMeanExtra(0,
3718  dataSizeB);
3719  tmpSeq.psd(0,
3720  (unsigned int) std::sqrt((double) dataSizeB), // numBlocks
3721  .5, // hopSizeRatio
3722  psdResult);
3723  double psdB = psdResult[0];
3724  double varOfMeanB = 2.*M_PI*psdB/doubleDataSizeB;
3725 
3726  if (m_env.subDisplayFile()) {
3727  *m_env.subDisplayFile() << "In ScalarSequence<T>::geweke()"
3728  << ", before computation of gewCoef"
3729  << ":\n"
3730  << ", dataSizeA = " << dataSizeA
3731  << ", numBlocksA = " << (unsigned int) std::sqrt((double) dataSizeA)
3732  << ", meanA = " << meanA
3733  << ", psdA = " << psdA
3734  << ", varOfMeanA = " << varOfMeanA
3735  << "\n"
3736  << ", dataSizeB = " << dataSizeB
3737  << ", numBlocksB = " << (unsigned int) std::sqrt((double) dataSizeB)
3738  << ", meanB = " << meanB
3739  << ", psdB = " << psdB
3740  << ", varOfMeanB = " << varOfMeanB
3741  << std::endl;
3742  }
3743  double gewCoef = (meanA - meanB)/std::sqrt(varOfMeanA + varOfMeanB);
3744 
3745  return gewCoef;
3746 }
3747 //---------------------------------------------------
3748 template <class T>
3749 T
3750 ScalarSequence<T>::meanStacc(
3751  unsigned int initialPos) const
3752 {
3753  queso_not_implemented(); // ERNESTO
3754 
3755  double value = 0.;
3756 
3757  return value;
3758 }
3759 //---------------------------------------------------
3760 template <class T>
3761 void
3762 ScalarSequence<T>::subCdfPercentageRange(
3763  unsigned int initialPos,
3764  unsigned int numPos,
3765  double range, // \in [0,1]
3766  T& lowerValue,
3767  T& upperValue) const
3768 {
3769  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
3770 
3771  queso_require_msg(!((range < 0) || (range > 1.)), "invalid 'range' value");
3772 
3773  ScalarSequence<T> sortedSequence(m_env,0,"");;
3774  sortedSequence.resizeSequence(numPos);
3775  this->extractScalarSeq(initialPos,
3776  1,
3777  numPos,
3778  sortedSequence);
3779  sortedSequence.subSort();
3780 
3781  unsigned int lowerId = (unsigned int) round( 0.5*(1.-range)*((double) numPos) );
3782  lowerValue = sortedSequence[lowerId];
3783 
3784  unsigned int upperId = (unsigned int) round( 0.5*(1.+range)*((double) numPos) );
3785  if (upperId == numPos) {
3786  upperId = upperId-1;
3787  }
3788  queso_require_less_msg(upperId, numPos, "'upperId' got too big");
3789  upperValue = sortedSequence[upperId];
3790 
3791  return;
3792 }
3793 //---------------------------------------------------
3794 template <class T>
3795 void
3796 ScalarSequence<T>::unifiedCdfPercentageRange(
3797  bool useOnlyInter0Comm,
3798  unsigned int initialPos,
3799  unsigned int numPos,
3800  double range, // \in [0,1]
3801  T& unifiedLowerValue,
3802  T& unifiedUpperValue) const
3803 {
3804  if (m_env.numSubEnvironments() == 1) {
3805  return this->subCdfPercentageRange(initialPos,
3806  numPos,
3807  range,
3808  unifiedLowerValue,
3809  unifiedUpperValue);
3810  }
3811 
3812  // As of 07/Feb/2011, this routine does *not* require sub sequences to have equal size. Good.
3813 
3814  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid input");
3815 
3816  queso_require_msg(!((range < 0) || (range > 1.)), "invalid 'range' value");
3817 
3818  if (useOnlyInter0Comm) {
3819  if (m_env.inter0Rank() >= 0) {
3820  queso_error_msg("code not implemented yet");
3821  }
3822  else {
3823  // Node not in the 'inter0' communicator
3824  this->subCdfPercentageRange(initialPos,
3825  numPos,
3826  unifiedLowerValue,
3827  unifiedUpperValue);
3828  }
3829  }
3830  else {
3831  queso_error_msg("parallel vectors not supported yet");
3832  }
3833  //m_env.fullComm().Barrier();
3834 
3835  return;
3836 }
3837 //---------------------------------------------------
3838 template <class T>
3839 void
3840 ScalarSequence<T>::subCdfStacc(
3841  unsigned int initialPos,
3842  std::vector<double>& cdfStaccValues,
3843  std::vector<double>& cdfStaccValuesUp,
3844  std::vector<double>& cdfStaccValuesLow,
3845  const ScalarSequence<T>& sortedDataValues) const
3846 {
3847  queso_require_msg(!(false), "not implemented yet"); // Joseph
3848  bool bRC = (initialPos < this->subSequenceSize() );
3849  queso_require_msg(bRC, "invalid input data");
3850 
3851  unsigned int numPoints = subSequenceSize()-initialPos;
3852  double auxNumPoints = numPoints;
3853  double maxLamb = 0.;
3854  std::vector<double> ro (numPoints,0.);
3855  std::vector<double> Isam_mat(numPoints,0.);
3856 
3857  for (unsigned int pointId = 0; pointId < numPoints; pointId++) {
3858  double p = ( ((double) pointId) + 1.0 )/auxNumPoints;
3859  double ro0 = p*(1.0-p);
3860  cdfStaccValues[pointId] = p;
3861 
3862  //std::cout << "x-data" << data[pointId]
3863  // << std::endl;
3864 
3865  for (unsigned int k = 0; k < numPoints; k++) {
3866  if (m_seq[k] <= sortedDataValues[pointId]) {
3867  Isam_mat[k] = 1;
3868  }
3869  else {
3870  Isam_mat[k] = 0;
3871  }
3872  }
3873 
3874  for (unsigned int tau = 0; tau < (numPoints-1); tau++) {
3875  ro[tau] = 0.;
3876  for (unsigned int kk = 0; kk < numPoints-(tau+1); kk++) {
3877  ro[tau] += (Isam_mat[kk+tau+1]-p)*(Isam_mat[kk]-p);
3878  }
3879  ro[tau] *= 1.0/auxNumPoints;
3880  }
3881  double lamb = 0.;
3882  for (unsigned int tau = 0; tau < (numPoints-1); tau++) {
3883  double auxTau = tau;
3884  lamb += (1.-(auxTau+1.)/auxNumPoints)*ro[tau]/ro0;
3885  if (lamb > maxLamb) maxLamb = lamb;
3886  }
3887  lamb = 2.*maxLamb;
3888 
3889  //double ll=gsl_cdf_gaussian_Pinv(-0.05);
3890  cdfStaccValuesUp [pointId] = cdfStaccValues[pointId]+1.96*sqrt(ro0/auxNumPoints*(1.+lamb));
3891  cdfStaccValuesLow[pointId] = cdfStaccValues[pointId]-1.96*sqrt(ro0/auxNumPoints*(1.+lamb));
3892  if (cdfStaccValuesLow[pointId] < 0.) cdfStaccValuesLow[pointId] = 0.;
3893  if (cdfStaccValuesUp [pointId] > 1.) cdfStaccValuesUp [pointId] = 1.;
3894  //if (pointId==1 | pointId==100 | pointId==900) {
3895  // std::cout<<lamb<<ll<<std::endl;
3896  // std::cout<<lamb<<std::endl;
3897  //}
3898  }
3899 
3900 #if 0
3901  unsigned int dataSize = this->subSequenceSize() - initialPos;
3902  unsigned int numEvals = evaluationPositions.size();
3903 
3904  for (unsigned int j = 0; j < numEvals; ++j) {
3905  double x = evaluationPositions[j];
3906  double value = 0.;
3907  for (unsigned int k = 0; k < dataSize; ++k) {
3908  double xk = m_seq[initialPos+k];
3909  value += 0.;//MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
3910  }
3911  cdfStaccValues[j] = value/(double) dataSize;
3912  }
3913 #endif
3914 
3915  return;
3916 }
3917 //---------------------------------------------------
3918 template <class T>
3919 void
3920 ScalarSequence<T>::subCdfStacc(
3921  unsigned int initialPos,
3922  const std::vector<T>& evaluationPositions,
3923  std::vector<double>& cdfStaccValues) const
3924 {
3925  queso_not_implemented();
3926 
3927  bool bRC = ((initialPos < this->subSequenceSize() ) &&
3928  (0 < evaluationPositions.size()) &&
3929  (evaluationPositions.size() == cdfStaccValues.size() ));
3930  queso_require_msg(bRC, "invalid input data");
3931 
3932  // For Joseph:
3933  // Maybe sort first
3934  // For each of the evaluation positions:
3935  // --> 1) form temporary scalar seq that contains only 0 and 1
3936  // --> 2) call "temporary scalar seq"->meanStacc()
3937 
3938 #if 0
3939  unsigned int dataSize = this->subSequenceSize() - initialPos;
3940  unsigned int numEvals = evaluationPositions.size();
3941 
3942  for (unsigned int j = 0; j < numEvals; ++j) {
3943  //double x = evaluationPositions[j];
3944  double value = 0.;
3945  for (unsigned int k = 0; k < dataSize; ++k) {
3946  //double xk = m_seq[initialPos+k];
3947  value += 0.;//MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
3948  }
3949  cdfStaccValues[j] = value/(double) dataSize;
3950  }
3951 #endif
3952 
3953  return;
3954 }
3955 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
3956 
3957 // --------------------------------------------------
3958 // Additional methods -------------------------------
3959 // (outside class declaration) ----------------------
3960 //---------------------------------------------------
3961 template <class T>
3962 void
3963 ComputeSubGaussian2dKde(const ScalarSequence<T>& scalarSeq1,
3964  const ScalarSequence<T>& scalarSeq2,
3965  unsigned int initialPos,
3966  double scaleValue1,
3967  double scaleValue2,
3968  const std::vector<T>& evaluationPositions1,
3969  const std::vector<T>& evaluationPositions2,
3970  std::vector<double>& densityValues)
3971 {
3972  queso_require_equal_to_msg(initialPos, 0, "not implemented yet for initialPos != 0");
3973 
3974  double covValue = 0.;
3975  double corrValue = 0.;
3976  ComputeCovCorrBetweenScalarSequences(scalarSeq1,
3977  scalarSeq2,
3978  scalarSeq1.subSequenceSize(),
3979  covValue,
3980  corrValue);
3981 
3982  bool bRC = ((initialPos < scalarSeq1.subSequenceSize()) &&
3983  (scalarSeq1.subSequenceSize() == scalarSeq2.subSequenceSize()) &&
3984  (0 < evaluationPositions1.size() ) &&
3985  (evaluationPositions1.size() == evaluationPositions2.size() ) &&
3986  (evaluationPositions1.size() == densityValues.size() ));
3987  queso_require_msg(bRC, "invalid input data");
3988 
3989  unsigned int dataSize = scalarSeq1.subSequenceSize() - initialPos;
3990  unsigned int numEvals = evaluationPositions1.size();
3991 
3992  double scale1Inv = 1./scaleValue1;
3993  double scale2Inv = 1./scaleValue2;
3994  //corrValue = 0.;
3995  double r = 1. - corrValue*corrValue;
3996  if (r <= 0.) { // prudencio 2010-07-23
3997  std::cerr << "In ComputeSubGaussian2dKde()"
3998  << ": WARNING, r = " << r
3999  << std::endl;
4000  }
4001 
4002  // scalarSeq1.env().worldRank(),
4003  // "ComputeSubGaussian2dKde()",
4004  // "negative r");
4005  for (unsigned int j = 0; j < numEvals; ++j) {
4006  double x1 = evaluationPositions1[j];
4007  double x2 = evaluationPositions2[j];
4008  double value = 0.;
4009  for (unsigned int k = 0; k < dataSize; ++k) {
4010  double d1k = scale1Inv*(x1 - scalarSeq1[initialPos+k]);
4011  double d2k = scale2Inv*(x2 - scalarSeq2[initialPos+k]);
4012  value += exp(-.5*( d1k*d1k + 2*corrValue*d1k*d2k + d2k*d2k )/r);
4013  }
4014  densityValues[j] = scale1Inv * scale2Inv * (value/(double) dataSize) / 2. / M_PI / sqrt(r);
4015  }
4016 
4017  return;
4018 }
4019 //---------------------------------------------------
4020 template <class T>
4021 void
4022 ComputeUnifiedGaussian2dKde(bool useOnlyInter0Comm, // INPUT
4023  const ScalarSequence<T>& scalarSeq1, // INPUT
4024  const ScalarSequence<T>& scalarSeq2, // INPUT
4025  unsigned int initialPos, // INPUT
4026  double unifiedScaleValue1, // INPUT
4027  double unifiedScaleValue2, // INPUT
4028  const std::vector<T>& unifiedEvaluationPositions1, // INPUT
4029  const std::vector<T>& unifiedEvaluationPositions2, // INPUT
4030  std::vector<double>& unifiedDensityValues) // OUTPUT
4031 {
4032  if (scalarSeq1.env().numSubEnvironments() == 1) {
4033  return ComputeSubGaussian2dKde(scalarSeq1,
4034  scalarSeq2,
4035  initialPos,
4036  unifiedScaleValue1,
4037  unifiedScaleValue2,
4038  unifiedEvaluationPositions1,
4039  unifiedEvaluationPositions2,
4040  unifiedDensityValues);
4041  }
4042 
4043  // As of 14/Nov/2009, this routine needs to be checked if it requires sub sequences to have equal size. Good.
4044 
4045  if (useOnlyInter0Comm) {
4046  if (scalarSeq1.env().inter0Rank() >= 0) {
4047  queso_error_msg("inter0 case not supported yet");
4048  }
4049  else {
4050  // Node not in the 'inter0' communicator
4051  ComputeSubGaussian2dKde(scalarSeq1,
4052  scalarSeq2,
4053  initialPos,
4054  unifiedScaleValue1,
4055  unifiedScaleValue2,
4056  unifiedEvaluationPositions1,
4057  unifiedEvaluationPositions2,
4058  unifiedDensityValues);
4059  }
4060  }
4061  else {
4062  queso_error_msg("parallel vectors not supported yet");
4063  }
4064 
4065  //scalarSeq1.env().fullComm().Barrier();
4066 
4067  return;
4068 }
4069 // --------------------------------------------------
4070 template <class T>
4071 void
4072 ComputeCovCorrBetweenScalarSequences(
4073  const ScalarSequence<T>& subPSeq,
4074  const ScalarSequence<T>& subQSeq,
4075  unsigned int subNumSamples,
4076  T& covValue,
4077  T& corrValue)
4078 {
4079  // Check input data consistency
4080  const BaseEnvironment& env = subPSeq.env();
4081 
4082  queso_require_msg(!((subNumSamples > subPSeq.subSequenceSize()) || (subNumSamples > subQSeq.subSequenceSize())), "subNumSamples is too large");
4083 
4084  // For both P and Q vector sequences: fill them
4085  T tmpP = 0.;
4086  T tmpQ = 0.;
4087 
4088  // For both P and Q vector sequences: compute the unified mean
4089  T unifiedMeanP = subPSeq.unifiedMeanExtra(true,0,subNumSamples);
4090  T unifiedMeanQ = subQSeq.unifiedMeanExtra(true,0,subNumSamples);
4091 
4092  // Compute "sub" covariance matrix
4093  covValue = 0.;
4094  for (unsigned k = 0; k < subNumSamples; ++k) {
4095  // For both P and Q vector sequences: get the difference (wrt the unified mean) in them
4096  tmpP = subPSeq[k] - unifiedMeanP;
4097  tmpQ = subQSeq[k] - unifiedMeanQ;
4098  covValue += tmpP*tmpQ;
4099  }
4100 
4101  // For both P and Q vector sequences: compute the unified variance
4102  T unifiedSampleVarianceP = subPSeq.unifiedSampleVarianceExtra(true,
4103  0,
4104  subNumSamples,
4105  unifiedMeanP);
4106 
4107  T unifiedSampleVarianceQ = subQSeq.unifiedSampleVarianceExtra(true,
4108  0,
4109  subNumSamples,
4110  unifiedMeanQ);
4111 
4112  // Compute unified covariance
4113  if (env.inter0Rank() >= 0) {
4114  unsigned int unifiedNumSamples = 0;
4115  env.inter0Comm().template Allreduce<unsigned int>(&subNumSamples, &unifiedNumSamples, (int) 1, RawValue_MPI_SUM,
4116  "ComputeCovCorrBetweenScalarSequences()",
4117  "failed MPI.Allreduce() for subNumSamples");
4118 
4119  double aux = 0.;
4120  env.inter0Comm().template Allreduce<double>(&covValue, &aux, (int) 1, RawValue_MPI_SUM,
4121  "ComputeCovCorrBetweenScalarSequences()",
4122  "failed MPI.Allreduce() for a matrix position");
4123  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)
4124 
4125  corrValue = covValue/std::sqrt(unifiedSampleVarianceP)/std::sqrt(unifiedSampleVarianceQ);
4126 
4127  if ((corrValue < -1.) || (corrValue > 1.)) { // prudencio 2010-07-23
4128  std::cerr << "In ComputeCovCorrBetweenScalarSequences()"
4129  << ": computed correlation is out of range, corrValue = " << corrValue
4130  << std::endl;
4131  }
4132 
4133  // env.worldRank(),
4134  // "ComputeCovCorrBetweenScalarSequences()",
4135  // "computed correlation is out of range");
4136  }
4137  else {
4138  // Node not in the 'inter0' communicator: do nothing extra
4139  }
4140 
4141  return;
4142 }
4143 
4144 } // End namespace QUESO
4145 
4146 template class QUESO::ScalarSequence<double>;
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:2876
QUESO::ScalarSequence::erasePositions
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values of the sequence, starting at position initialPos.
Definition: ScalarSequence.C:206
QUESO::ScalarSequence::subMeanExtra
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
Definition: ScalarSequence.C:838
QUESO::UQ_OK_RC
const int UQ_OK_RC
Definition: Defines.h:89
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:2340
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:2509
QUESO::MiscGetEllapsedSeconds
double MiscGetEllapsedSeconds(struct timeval *timeval0)
Definition: Miscellaneous.C:282
QUESO::RawType_MPI_Status
int RawType_MPI_Status
Definition: MpiComm.h:62
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::ScalarSequence::unifiedSort
void unifiedSort(bool useOnlyInter0Comm, unsigned int initialPos, ScalarSequence< T > &unifiedSortedSequence) const
Sorts the unified sequence of scalars.
Definition: ScalarSequence.C:1891
RawValue_MPI_MIN
#define RawValue_MPI_MIN
Definition: MpiComm.h:69
RawValue_MPI_SUM
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
queso_require_greater_equal_msg
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:78
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:1530
QUESO::ScalarSequence::unifiedMeanPlain
const T & unifiedMeanPlain(bool useOnlyInter0Comm) const
Finds the mean value of the unified sequence of scalars.
Definition: ScalarSequence.C:407
QUESO::FilePtrSetStruct::ifsVar
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:88
QUESO::ScalarSequence::unifiedMinPlain
const T & unifiedMinPlain(bool useOnlyInter0Comm) const
Finds the minimum value of the unified sequence of scalars.
Definition: ScalarSequence.C:356
QUESO::ScalarSequence::autoCorrViaFft
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
Definition: ScalarSequence.C:1363
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::operator[]
const T & operator[](unsigned int posId) const
Access position posId of the sequence of scalars (const).
Definition: ScalarSequence.C:69
QUESO::ScalarSequence::unifiedReadContents
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
Definition: ScalarSequence.C:2897
QUESO::ScalarSequence::rawData
std::vector< T > & rawData()
The sequence of scalars. Access to private attribute m_seq.
Definition: ScalarSequence.C:3247
QUESO::ScalarSequence::seqScalarPositionConstIteratorTypedef
std::vector< T >::const_iterator seqScalarPositionConstIteratorTypedef
Definition: ScalarSequence.h:60
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:2239
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:4022
QUESO::BaseEnvironment
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:197
QUESO::ScalarSequence::clear
void clear()
Clears the sequence of scalars.
Definition: ScalarSequence.C:123
QUESO::ScalarSequence::append
void append(const ScalarSequence< T > &src, unsigned int srcInitialPos, unsigned int srcNumPos)
Appends the scalar sequence src to this sequence.
Definition: ScalarSequence.C:2488
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::operator=
ScalarSequence< T > & operator=(const ScalarSequence< T > &rhs)
Assignment operator; it copies rhs to this.
Definition: ScalarSequence.C:61
QUESO::BaseEnvironment::inter0Rank
int inter0Rank() const
Returns the process inter0 rank.
Definition: Environment.C:307
QUESO::MiscGaussianDensity
double MiscGaussianDensity(double x, double mu, double sigma)
Definition: Miscellaneous.C:306
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::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
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:1738
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:2459
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:3224
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_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:73
UQ_FILE_EXTENSION_FOR_HDF_FORMAT
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:104
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:3265
QUESO::ScalarSequence::name
const std::string & name() const
Access to the name of the sequence of scalars.
Definition: ScalarSequence.C:108
queso_require_less_msg
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:75
QUESO::ScalarSequence::unifiedSequenceSize
unsigned int unifiedSequenceSize(bool useOnlyInter0Comm) const
Size of the unified sequence of scalars.
Definition: ScalarSequence.C:143
QUESO::ScalarSequence::setName
void setName(const std::string &newName)
Sets a new name to the sequence of scalars.
Definition: ScalarSequence.C:115
RawValue_MPI_MAX
#define RawValue_MPI_MAX
Definition: MpiComm.h:70
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:1648
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:2539
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
Class for handling scalar samples.
Definition: ScalarSequence.h:54
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:2108
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
k
int k
Definition: ann_sample.cpp:53
UQ_FILE_EXTENSION_FOR_TXT_FORMAT
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
QUESO::FilePtrSetStruct
Struct for handling data input and output from files.
Definition: Environment.h:76
QUESO::ScalarSequence::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
Definition: ScalarSequence.C:2420
QUESO::ScalarSequence::subMinPlain
const T & subMinPlain() const
Finds the minimum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:343
QUESO::ScalarSequence::ScalarSequence
ScalarSequence(const BaseEnvironment &env, unsigned int subSequenceSize, const std::string &name)
Default constructor.
Definition: ScalarSequence.C:32
QUESO::MiscHammingWindow
double MiscHammingWindow(unsigned int N, unsigned int j)
Definition: Miscellaneous.C:298
QUESO::UniformOneDGrid
Class for accommodating uniform one-dimensional grids.
Definition: UniformOneDGrid.h:41
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
RawValue_MPI_DOUBLE
#define RawValue_MPI_DOUBLE
Definition: MpiComm.h:67
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:2310
RawValue_MPI_UNSIGNED
#define RawValue_MPI_UNSIGNED
Definition: MpiComm.h:68
QUESO::ScalarSequence::copy
void copy(const ScalarSequence< T > &src)
Copies the scalar sequence src to this.
Definition: ScalarSequence.C:3168
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::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:3963
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::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:2569
QUESO::ScalarSequence::subSort
void subSort(unsigned int initialPos, ScalarSequence< T > &sortedSequence) const
Sorts the sub-sequence of scalars.
Definition: ScalarSequence.C:1874
QUESO::ScalarSequence::subInterQuantileRange
T subInterQuantileRange(unsigned int initialPos) const
Returns the interquartile range of the values in the sub-sequence.
Definition: ScalarSequence.C:1995
QUESO::ScalarSequence::subMeanPlain
const T & subMeanPlain() const
Finds the mean value of the sub-sequence of scalars.
Definition: ScalarSequence.C:395
RawValue_MPI_ANY_SOURCE
#define RawValue_MPI_ANY_SOURCE
Definition: MpiComm.h:64
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::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::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_not_implemented
#define queso_not_implemented()
Definition: asserts.h:56
QUESO::ScalarSequence::unifiedMaxPlain
const T & unifiedMaxPlain(bool useOnlyInter0Comm) const
Finds the maximum value of the unified sequence of scalars.
Definition: ScalarSequence.C:382
QUESO::ScalarSequence::~ScalarSequence
~ScalarSequence()
Destructor.
Definition: ScalarSequence.C:54
QUESO::ScalarSequence::subScaleForKde
T subScaleForKde(unsigned int initialPos, const T &iqrValue, unsigned int kdeDimension) const
Selects the scales (output value) for the kernel density estimation, considering only the sub-sequenc...
Definition: ScalarSequence.C:2199
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:3183
QUESO::BaseEnvironment::numSubEnvironments
unsigned int numSubEnvironments() const
Access function to the number of sub-environments.
Definition: Environment.C:334
QUESO::ScalarSequence::unifiedWriteContents
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
Writes the unified sequence to a file.
Definition: ScalarSequence.C:2672
QUESO::FilePtrSetStruct::ofsVar
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:85
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
QUESO::ScalarSequence::subSampleVariancePlain
const T & subSampleVariancePlain() const
Finds the variance of a sample of the sub-sequence of scalars.
Definition: ScalarSequence.C:443
QUESO::ScalarSequence::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_require_not_equal_to_msg
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:74
QUESO::ScalarSequence::deleteStoredScalars
void deleteStoredScalars()
Deletes all stored scalars.
Definition: ScalarSequence.C:467
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:1780
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
SCALAR_SEQUENCE_INIT_MPI_MSG
#define SCALAR_SEQUENCE_INIT_MPI_MSG
Definition: ScalarSequence.h:37
dataIn
istream * dataIn
Definition: ann_sample.cpp:58
queso_require_less_equal_msg
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:77
QUESO::ScalarSequence::m_seq
std::vector< T > m_seq
Definition: ScalarSequence.h:527
QUESO::ScalarSequence::setGaussian
void setGaussian(const T &mean, const T &stdDev)
Sets the values of the sequence as a Gaussian distribution of mean given by meanVec and standard devi...
Definition: ScalarSequence.C:515
QUESO::ScalarSequence::seqScalarPositionIteratorTypedef
std::vector< T >::iterator seqScalarPositionIteratorTypedef
Definition: ScalarSequence.h:59
UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
queso_require_msg
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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:2864
QUESO::ScalarSequence::subSort
void subSort()
Sorts the sequence of scalars in the private attribute m_seq.
Definition: ScalarSequence.C:3255
SCALAR_SEQUENCE_SIZE_MPI_MSG
#define SCALAR_SEQUENCE_SIZE_MPI_MSG
Definition: ScalarSequence.h:38
QUESO::ComputeCovCorrBetweenScalarSequences
void ComputeCovCorrBetweenScalarSequences(const ScalarSequence< T > &subPSeq, const ScalarSequence< T > &subQSeq, unsigned int subNumSamples, T &covValue, T &corrValue)
Definition: ScalarSequence.C:4072
QUESO::ScalarSequence::unifiedMedianPlain
const T & unifiedMedianPlain(bool useOnlyInter0Comm) const
Finds the median value of the unified sequence of scalars.
Definition: ScalarSequence.C:431
QUESO::ScalarSequence::resetValues
void resetValues(unsigned int initialPos, unsigned int)
Sets numPos values of the sequence to zero, starting at position initialPos.
Definition: ScalarSequence.C:186
QUESO::ScalarSequence::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:1501
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_error_msg
#define queso_error_msg(msg)
Definition: asserts.h:47
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
QUESO::ScalarSequence::writeTxtHeader
void writeTxtHeader(std::ofstream &ofs, double sequenceSize) const
Helper function to write txt info for matlab files.
Definition: ScalarSequence.C:2888
QUESO::ScalarSequence::env
const BaseEnvironment & env() const
Access to QUESO environment.
Definition: ScalarSequence.C:101
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::autoCorrViaDef
T autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag) const
Calculates the autocorrelation via definition.
Definition: ScalarSequence.C:1334
SCALAR_SEQUENCE_DATA_MPI_MSG
#define SCALAR_SEQUENCE_DATA_MPI_MSG
Definition: ScalarSequence.h:39
QUESO::ScalarSequence::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of scalars.
Definition: ScalarSequence.C:136
QUESO::ScalarSequence::subMedianPlain
const T & subMedianPlain() const
Finds the median value of the sub-sequence of scalars.
Definition: ScalarSequence.C:419
QUESO::BaseEnvironment::inter0Comm
const MpiComm & inter0Comm() const
Access function for MpiComm inter0-communicator.
Definition: Environment.C:313
QUESO::Fft
Class for a Fast Fourier Transform (FFT) algorithm.
Definition: Fft.h:66
QUESO::ScalarSequence::subMaxPlain
const T & subMaxPlain() const
Finds the maximum value of the sub-sequence of scalars.
Definition: ScalarSequence.C:369
QUESO::ScalarSequence::subHistogram
void subHistogram(unsigned int initialPos, const T &minHorizontalValue, const T &maxHorizontalValue, std::vector< T > &centers, std::vector< unsigned int > &bins) const
Calculates the histogram of the sub-sequence.
Definition: ScalarSequence.C:1601
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
Generated on Tue Nov 29 2016 10:53:09 for queso-0.56.0 by  doxygen 1.8.5