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