queso-0.52.0
SequenceOfVectors.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/SequenceOfVectors.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Default constructor -----------------------------
32 template <class V, class M>
34  const VectorSpace<V,M>& vectorSpace,
35  unsigned int subSequenceSize,
36  const std::string& name)
37  :
38  BaseVectorSequence<V,M>(vectorSpace,subSequenceSize,name),
39  m_seq (subSequenceSize,NULL)
40 #ifdef UQ_CODE_HAS_MONITORS
41  ,
42  m_subMeanMonitorPosSeq (NULL),
43  m_subMeanVecSeq (NULL),
44  m_subMeanCltStdSeq (NULL),
45  m_subMeanInter0MonitorPosSeq (NULL),
46  m_subMeanInter0Mean (NULL),
47  m_subMeanInter0Clt95 (NULL),
48  m_subMeanInter0Empirical90 (NULL),
49  m_subMeanInter0Min (NULL),
50  m_subMeanInter0Max (NULL),
51  m_unifiedMeanMonitorPosSeq (NULL),
52  m_unifiedMeanVecSeq (NULL),
53  m_unifiedMeanCltStdSeq (NULL)
54 #endif
55 {
56  //if (m_env.subDisplayFile()) {
57  // *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::constructor()"
58  // << std::endl;
59  //}
60 
61  //if (m_env.subDisplayFile()) {
62  // *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::constructor()"
63  // << std::endl;
64  //}
65 }
66 // Destructor ---------------------------------------
67 template <class V, class M>
69 {
70 #ifdef UQ_CODE_HAS_MONITORS
71  if (m_subMeanMonitorPosSeq) delete m_subMeanMonitorPosSeq;
72  if (m_subMeanVecSeq ) delete m_subMeanVecSeq;
73  if (m_subMeanCltStdSeq ) delete m_subMeanCltStdSeq;
74 
75  if (m_subMeanInter0MonitorPosSeq) delete m_subMeanInter0MonitorPosSeq;
76  if (m_subMeanInter0Mean ) delete m_subMeanInter0Mean;
77  if (m_subMeanInter0Clt95 ) delete m_subMeanInter0Clt95;
78  if (m_subMeanInter0Empirical90 ) delete m_subMeanInter0Empirical90;
79  if (m_subMeanInter0Min ) delete m_subMeanInter0Min;
80  if (m_subMeanInter0Max ) delete m_subMeanInter0Max;
81 
82  if (m_unifiedMeanMonitorPosSeq) delete m_unifiedMeanMonitorPosSeq;
83  if (m_unifiedMeanVecSeq ) delete m_unifiedMeanVecSeq;
84  if (m_unifiedMeanCltStdSeq ) delete m_unifiedMeanCltStdSeq;
85 #endif
86 
87  for (unsigned int i = 0; i < (unsigned int) m_seq.size(); ++i) {
88  if (m_seq[i]) delete m_seq[i];
89  }
90 }
91 // Set methods --------------------------------------
92 template <class V, class M>
95 {
96  this->copy(rhs);
97  return *this;
98 }
99 
100 // Sequence methods ---------------------------------
101 template <class V, class M>
102 unsigned int
104 {
105  return m_seq.size();
106 }
107 //---------------------------------------------------
108 template <class V, class M>
109 void
110 SequenceOfVectors<V,M>::resizeSequence(unsigned int newSubSequenceSize)
111 {
112  if (newSubSequenceSize != this->subSequenceSize()) {
113  if (newSubSequenceSize < this->subSequenceSize()) {
114  this->resetValues(newSubSequenceSize,this->subSequenceSize()-newSubSequenceSize);
115  }
116  m_seq.resize(newSubSequenceSize,NULL);
117  std::vector<const V*>(m_seq).swap(m_seq);
119  }
120 
121  return;
122 }
123 //---------------------------------------------------
124 template <class V, class M>
125 void
126 SequenceOfVectors<V,M>::resetValues(unsigned int initialPos, unsigned int numPos)
127 {
128  bool bRC = ((initialPos < this->subSequenceSize()) &&
129  (0 < numPos ) &&
130  ((initialPos+numPos) <= this->subSequenceSize()));
131  if ((bRC == false) && (m_env.subDisplayFile())) {
132  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::resetValues()"
133  << ", initialPos = " << initialPos
134  << ", this->subSequenceSize() = " << this->subSequenceSize()
135  << ", numPos = " << numPos
136  << std::endl;
137  }
138  UQ_FATAL_TEST_MACRO(bRC == false,
139  m_env.worldRank(),
140  "SequenceOfVectors<V,M>::resetValues()",
141  "invalid input data");
142 
143  for (unsigned int j = 0; j < numPos; ++j) {
144  if (m_seq[initialPos+j] != NULL) {
145  delete m_seq[initialPos+j];
146  m_seq[initialPos+j] = NULL;
147  }
148  }
149 
151 
152  return;
153 }
154 //---------------------------------------------------
155 template <class V, class M>
156 void
157 SequenceOfVectors<V,M>::erasePositions(unsigned int initialPos, unsigned int numPos)
158 {
159  bool bRC = ((initialPos < this->subSequenceSize()) &&
160  (0 < numPos ) &&
161  ((initialPos+numPos) <= this->subSequenceSize()));
162  UQ_FATAL_TEST_MACRO(bRC == false,
163  m_env.worldRank(),
164  "SequenceOfVectors<V,M>::erasePositions()",
165  "invalid input data");
166 
167  for (unsigned int j = 0; j < numPos; ++j) {
168  if (m_seq[initialPos+j] != NULL) {
169  delete m_seq[initialPos+j];
170  m_seq[initialPos+j] = NULL;
171  }
172  }
173 
174  seqVectorPositionIteratorTypedef posIteratorBegin = m_seq.begin();
175  if (initialPos < this->subSequenceSize()) std::advance(posIteratorBegin,initialPos);
176  else posIteratorBegin = m_seq.end();
177 
178  unsigned int posEnd = initialPos + numPos - 1;
179  seqVectorPositionIteratorTypedef posIteratorEnd = m_seq.begin();
180  if (posEnd < this->subSequenceSize()) std::advance(posIteratorEnd,posEnd);
181  else posIteratorEnd = m_seq.end();
182 
183  unsigned int oldSubSequenceSize = this->subSequenceSize();
184  m_seq.erase(posIteratorBegin,posIteratorEnd);
185  UQ_FATAL_TEST_MACRO((oldSubSequenceSize - numPos) != this->subSequenceSize(),
186  m_env.worldRank(),
187  "SequenceOfVectors::erasePositions()",
188  "(oldSubSequenceSize - numPos) != this->subSequenceSize()");
189 
191 
192  return;
193 }
194 //---------------------------------------------------
195 template <class V, class M>
196 void
197 SequenceOfVectors<V,M>::getPositionValues(unsigned int posId, V& vec) const
198 {
199  UQ_FATAL_TEST_MACRO((posId >= this->subSequenceSize()),
200  m_env.worldRank(),
201  "SequenceOfVectorss<V,M>::getPositionValues()",
202  "posId > subSequenceSize()");
203 
204  UQ_FATAL_TEST_MACRO(m_seq[posId] == NULL,
205  m_env.worldRank(),
206  "SequenceOfVectorss<V,M>::getPositionValues()",
207  "posId is NULL");
208 
209  //if (posId == 0) { // mox
210  // std::cout << "In SequenceOfVectors<V,M>::getPositionValues(): m_seq[0] = " << m_seq[0] << ", *(m_seq[0]) = " << *(m_seq[0])
211  // << std::endl;
212  //}
213 
214  vec = *(m_seq[posId]); // *(const_cast<V*>(m_seq[posId])); // prudenci 2010-06-17 mox
215 
216  return;
217 }
218 //---------------------------------------------------
219 template <class V, class M>
220 void
221 SequenceOfVectors<V,M>::setPositionValues(unsigned int posId, const V& vec)
222 {
223  UQ_FATAL_TEST_MACRO((posId >= this->subSequenceSize()),
224  m_env.worldRank(),
225  "SequenceOfVectorss<V,M>::setPositionValues()",
226  "posId > subSequenceSize()");
227 
228  UQ_FATAL_TEST_MACRO(vec.sizeLocal() != m_vectorSpace.zeroVector().sizeLocal(),
229  m_env.worldRank(),
230  "SequenceOfVectorss<V,M>::setPositionValues()",
231  "invalid vec");
232 
233  if (m_seq[posId] != NULL) delete m_seq[posId];
234  m_seq[posId] = new V(vec);
235 
236  //if (posId == 0) { // mox
237  // std::cout << "In SequenceOfVectors<V,M>::setPositionValues(): m_seq[0] = " << m_seq[0] << ", *(m_seq[0]) = " << *(m_seq[0])
238  // << std::endl;
239  //}
240 
242 
243  return;
244 }
245 //---------------------------------------------------
246 template <class V, class M>
247 void
249  const V& numEvaluationPointsVec,
250  ArrayOfOneDGrids <V,M>& cdfGrids,
251  ArrayOfOneDTables<V,M>& cdfValues) const
252 {
253  V minDomainValues(m_vectorSpace.zeroVector());
254  V maxDomainValues(m_vectorSpace.zeroVector());
255 
256  ScalarSequence<double> data(m_env,0,"");
257 
258  unsigned int numParams = this->vectorSizeLocal();
259  for (unsigned int i = 0; i < numParams; ++i) {
260  this->extractScalarSeq(0, // initialPos
261  1, // spacing
262  subSequenceSize(), // numPos
263  i,
264  data);
265 
266  std::vector<double> aCdf(0);
267  data.subUniformlySampledCdf((unsigned int) numEvaluationPointsVec[i],
268  minDomainValues[i],
269  maxDomainValues[i],
270  aCdf);
271  cdfValues.setOneDTable(i,aCdf);
272  }
273 
274  cdfGrids.setUniformGrids(numEvaluationPointsVec,
275  minDomainValues,
276  maxDomainValues);
277 
278  return;
279 }
280 //---------------------------------------------------
281 template <class V, class M>
282 void
284  const V& numEvaluationPointsVec,
285  ArrayOfOneDGrids <V,M>& unifiedCdfGrids,
286  ArrayOfOneDTables<V,M>& unifiedCdfValues) const
287 {
288  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
289  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedUniformlySampledCdf()"
290  << std::endl;
291  }
292 
293  V unifiedMinDomainValues(m_vectorSpace.zeroVector());
294  V unifiedMaxDomainValues(m_vectorSpace.zeroVector());
295 
296  ScalarSequence<double> data(m_env,0,"");
297 
298  unsigned int numParams = this->vectorSizeLocal();
299  for (unsigned int i = 0; i < numParams; ++i) {
300  this->extractScalarSeq(0, // initialPos
301  1, // spacing
302  subSequenceSize(), // numPos
303  i,
304  data);
305 
306  std::vector<double> aCdf(0);
307  data.unifiedUniformlySampledCdf(m_vectorSpace.numOfProcsForStorage() == 1,
308  (unsigned int) numEvaluationPointsVec[i],
309  unifiedMinDomainValues[i],
310  unifiedMaxDomainValues[i],
311  aCdf);
312  unifiedCdfValues.setOneDTable(i,aCdf);
313  }
314 
315  unifiedCdfGrids.setUniformGrids(numEvaluationPointsVec,
316  unifiedMinDomainValues,
317  unifiedMaxDomainValues);
318 
319  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
320  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedUniformlySampledCdf()"
321  << std::endl;
322  }
323 
324  return;
325 }
326 //---------------------------------------------------
327 template <class V, class M>
328 void
330  unsigned int initialPos,
331  unsigned int numPos,
332  V& meanVec) const
333 {
334  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
335  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::subMeanExtra()"
336  << ": initialPos = " << initialPos
337  << ", numPos = " << numPos
338  << ", sub sequence size = " << this->subSequenceSize()
339  << std::endl;
340  }
341 
342  bool bRC = ((initialPos < this->subSequenceSize()) &&
343  (0 < numPos ) &&
344  ((initialPos+numPos) <= this->subSequenceSize()) &&
345  (this->vectorSizeLocal() == meanVec.sizeLocal() ));
346  if ((bRC == false) && (m_env.subDisplayFile())) {
347  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::subMeanExtra()"
348  << ", initialPos = " << initialPos
349  << ", this->subSequenceSize() = " << this->subSequenceSize()
350  << ", numPos = " << numPos
351  << ", this->vectorSizeLocal() = " << this->vectorSizeLocal()
352  << ", meanVec.sizeLocal() = " << meanVec.sizeLocal()
353  << std::endl;
354  }
355  UQ_FATAL_TEST_MACRO(bRC == false,
356  m_env.worldRank(),
357  "SequenceOfVectors<V,M>::subMeanExtra()",
358  "invalid input data");
359 
360  ScalarSequence<double> data(m_env,0,"");
361 
362  unsigned int numParams = this->vectorSizeLocal();
363  for (unsigned int i = 0; i < numParams; ++i) {
364  this->extractScalarSeq(initialPos,
365  1, // spacing
366  numPos,
367  i,
368  data);
369  meanVec[i] = data.subMeanExtra(0,
370  numPos);
371  }
372 
373  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
374  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::subMeanExtra()"
375  << ": initialPos = " << initialPos
376  << ", numPos = " << numPos
377  << ", sub sequence size = " << this->subSequenceSize()
378  << ", meanVec = " << meanVec
379  << std::endl;
380  }
381 
382  return;
383 }
384 //---------------------------------------------------
385 template <class V, class M>
386 void
388  unsigned int initialPos,
389  unsigned int numPos,
390  V& unifiedMeanVec) const
391 {
392  unsigned int tmpUnif = this->unifiedSequenceSize();
393  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
394  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedMeanExtra()"
395  << ": initialPos = " << initialPos
396  << ", numPos = " << numPos
397  << ", sub sequence size = " << this->subSequenceSize()
398  << ", unified sequence size = " << tmpUnif
399  << std::endl;
400  }
401 
402  bool bRC = ((initialPos < this->subSequenceSize() ) &&
403  (0 < numPos ) &&
404  ((initialPos+numPos) <= this->subSequenceSize() ) &&
405  (this->vectorSizeLocal() == unifiedMeanVec.sizeLocal()));
406  if ((bRC == false) && (m_env.subDisplayFile())) {
407  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedMeanExtra()"
408  << ", initialPos = " << initialPos
409  << ", this->subSequenceSize() = " << this->subSequenceSize()
410  << ", numPos = " << numPos
411  << ", this->vectorSizeLocal() = " << this->vectorSizeLocal()
412  << ", unifiedMeanVec.sizeLocal() = " << unifiedMeanVec.sizeLocal()
413  << std::endl;
414  }
415  UQ_FATAL_TEST_MACRO(bRC == false,
416  m_env.worldRank(),
417  "SequenceOfVectors<V,M>::unifiedMeanExtra()",
418  "invalid input data");
419 
420  ScalarSequence<double> data(m_env,0,"");
421 
422  unsigned int numParams = this->vectorSizeLocal();
423  for (unsigned int i = 0; i < numParams; ++i) {
424  this->extractScalarSeq(initialPos,
425  1, // spacing
426  numPos,
427  i,
428  data);
429  unifiedMeanVec[i] = data.unifiedMeanExtra(m_vectorSpace.numOfProcsForStorage() == 1,
430  0,
431  numPos);
432  }
433 
434  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
435  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedMeanExtra()"
436  << ": initialPos = " << initialPos
437  << ", numPos = " << numPos
438  << ", sub sequence size = " << this->subSequenceSize()
439  << ", unified sequence size = " << tmpUnif
440  << ", unifiedMeanVec = " << unifiedMeanVec
441  << std::endl;
442  }
443 
444  return;
445 }
446 //---------------------------------------------------
447 template <class V, class M>
448 void
450  unsigned int initialPos,
451  unsigned int numPos,
452  V& medianVec) const
453 {
454  if (this->subSequenceSize() == 0) return;
455 
456  bool bRC = ((initialPos < this->subSequenceSize()) &&
457  (0 < numPos ) &&
458  ((initialPos+numPos) <= this->subSequenceSize()));
459  if (bRC == false) {
460  std::cerr << "In SequenceOfVectors<V,M>::subMedianExtra()"
461  << ": ERROR at fullRank " << m_env.fullRank()
462  << ", initialPos = " << initialPos
463  << ", numPos = " << numPos
464  << ", this->subSequenceSize() = " << this->subSequenceSize()
465  << std::endl;
466  if (m_env.subDisplayFile()) {
467  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::subMedianExtra()"
468  << ": ERROR at fullRank " << m_env.fullRank()
469  << ", initialPos = " << initialPos
470  << ", numPos = " << numPos
471  << ", this->subSequenceSize() = " << this->subSequenceSize()
472  << std::endl;
473  }
474  }
475  UQ_FATAL_TEST_MACRO(bRC == false,
476  m_env.worldRank(),
477  "SequenceOfVectors<V,M>::subMedianExtra()",
478  "invalid input data");
479 
480  ScalarSequence<double> data(m_env,0,"");
481 
482  unsigned int numParams = this->vectorSizeLocal();
483  for (unsigned int i = 0; i < numParams; ++i) {
484  this->extractScalarSeq(initialPos,
485  1, // spacing
486  numPos,
487  i,
488  data);
489  medianVec[i] = data.subMedianExtra(0,
490  numPos);
491  }
492 
493  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
494  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::subMedianExtra()"
495  << ": initialPos = " << initialPos
496  << ", numPos = " << numPos
497  << ", sub sequence size = " << this->subSequenceSize()
498  << ", medianVec = " << medianVec
499  << std::endl;
500  }
501 
502  return;
503 }
504 //---------------------------------------------------
505 template <class V, class M>
506 void
508  unsigned int initialPos,
509  unsigned int numPos,
510  V& unifiedMedianVec) const
511 {
512  unsigned int tmpUnif = this->unifiedSequenceSize();
513  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
514  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedMedianExtra()"
515  << ": initialPos = " << initialPos
516  << ", numPos = " << numPos
517  << ", sub sequence size = " << this->subSequenceSize()
518  << ", unified sequence size = " << tmpUnif
519  << std::endl;
520  }
521 
522  bool bRC = ((initialPos < this->subSequenceSize() ) &&
523  (0 < numPos ) &&
524  ((initialPos+numPos) <= this->subSequenceSize() ) &&
525  (this->vectorSizeLocal() == unifiedMedianVec.sizeLocal()));
526  if ((bRC == false) && (m_env.subDisplayFile())) {
527  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedMedianExtra()"
528  << ", initialPos = " << initialPos
529  << ", this->subSequenceSize() = " << this->subSequenceSize()
530  << ", numPos = " << numPos
531  << ", this->vectorSizeLocal() = " << this->vectorSizeLocal()
532  << ", unifiedMedianVec.sizeLocal() = " << unifiedMedianVec.sizeLocal()
533  << std::endl;
534  }
535  UQ_FATAL_TEST_MACRO(bRC == false,
536  m_env.worldRank(),
537  "SequenceOfVectors<V,M>::unifiedMedianExtra()",
538  "invalid input data");
539 
540  ScalarSequence<double> data(m_env,0,"");
541 
542  unsigned int numParams = this->vectorSizeLocal();
543  for (unsigned int i = 0; i < numParams; ++i) {
544  this->extractScalarSeq(initialPos,
545  1, // spacing
546  numPos,
547  i,
548  data);
549  unifiedMedianVec[i] = data.unifiedMedianExtra(m_vectorSpace.numOfProcsForStorage() == 1,
550  0,
551  numPos);
552  }
553 
554  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
555  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedMedianExtra()"
556  << ": initialPos = " << initialPos
557  << ", numPos = " << numPos
558  << ", sub sequence size = " << this->subSequenceSize()
559  << ", unified sequence size = " << tmpUnif
560  << ", unifiedMedianVec = " << unifiedMedianVec
561  << std::endl;
562  }
563 
564  return;
565 }
566 //---------------------------------------------------
567 template <class V, class M>
568 void
570  unsigned int initialPos,
571  unsigned int numPos,
572  const V& meanVec,
573  V& samVec) const
574 {
575  bool bRC = ((initialPos < this->subSequenceSize()) &&
576  (0 < numPos ) &&
577  ((initialPos+numPos) <= this->subSequenceSize()) &&
578  (this->vectorSizeLocal() == meanVec.sizeLocal() ) &&
579  (this->vectorSizeLocal() == samVec.sizeLocal() ));
580  UQ_FATAL_TEST_MACRO(bRC == false,
581  m_env.worldRank(),
582  "SequenceOfVectors<V,M>::subSampleVarianceExtra()",
583  "invalid input data");
584 
585  ScalarSequence<double> data(m_env,0,"");
586 
587  unsigned int numParams = this->vectorSizeLocal();
588  for (unsigned int i = 0; i < numParams; ++i) {
589  this->extractScalarSeq(initialPos,
590  1, // spacing
591  numPos,
592  i,
593  data);
594  samVec[i] = data.subSampleVarianceExtra(0,
595  numPos,
596  meanVec[i]);
597  }
598 
599  return;
600 }
601 //---------------------------------------------------
602 template <class V, class M>
603 void
605  unsigned int initialPos,
606  unsigned int numPos,
607  const V& unifiedMeanVec,
608  V& unifiedSamVec) const
609 {
610  bool bRC = ((initialPos < this->subSequenceSize() ) &&
611  (0 < numPos ) &&
612  ((initialPos+numPos) <= this->subSequenceSize() ) &&
613  (this->vectorSizeLocal() == unifiedMeanVec.sizeLocal()) &&
614  (this->vectorSizeLocal() == unifiedSamVec.sizeLocal() ));
615  UQ_FATAL_TEST_MACRO(bRC == false,
616  m_env.worldRank(),
617  "SequenceOfVectors<V,M>::unifiedSampleVarianceExtra()",
618  "invalid input data");
619 
620  ScalarSequence<double> data(m_env,0,"");
621 
622  unsigned int numParams = this->vectorSizeLocal();
623  for (unsigned int i = 0; i < numParams; ++i) {
624  this->extractScalarSeq(initialPos,
625  1, // spacing
626  numPos,
627  i,
628  data);
629  unifiedSamVec[i] = data.unifiedSampleVarianceExtra(m_vectorSpace.numOfProcsForStorage() == 1,
630  0,
631  numPos,
632  unifiedMeanVec[i]);
633  }
634 
635  return;
636 }
637 //---------------------------------------------------
638 template <class V, class M>
639 void
641  unsigned int initialPos,
642  unsigned int numPos,
643  const V& meanVec,
644  V& stdvec) const
645 {
646  bool bRC = ((initialPos < this->subSequenceSize()) &&
647  (0 < numPos ) &&
648  ((initialPos+numPos) <= this->subSequenceSize()) &&
649  (this->vectorSizeLocal() == meanVec.sizeLocal() ) &&
650  (this->vectorSizeLocal() == stdvec.sizeLocal() ));
651  UQ_FATAL_TEST_MACRO(bRC == false,
652  m_env.worldRank(),
653  "SequenceOfVectors<V,M>::subSampleStd()",
654  "invalid input data");
655 
656  ScalarSequence<double> data(m_env,0,"");
657 
658  unsigned int numParams = this->vectorSizeLocal();
659  for (unsigned int i = 0; i < numParams; ++i) {
660  this->extractScalarSeq(initialPos,
661  1, // spacing
662  numPos,
663  i,
664  data);
665  stdvec[i] = data.subSampleStd(0,
666  numPos,
667  meanVec[i]);
668  }
669 
670  return;
671 }
672 //---------------------------------------------------
673 template <class V, class M>
674 void
676  unsigned int initialPos,
677  unsigned int numPos,
678  const V& unifiedMeanVec,
679  V& unifiedStdVec) const
680 {
681  bool bRC = ((initialPos < this->subSequenceSize() ) &&
682  (0 < numPos ) &&
683  ((initialPos+numPos) <= this->subSequenceSize() ) &&
684  (this->vectorSizeLocal() == unifiedMeanVec.sizeLocal()) &&
685  (this->vectorSizeLocal() == unifiedStdVec.sizeLocal() ));
686  UQ_FATAL_TEST_MACRO(bRC == false,
687  m_env.worldRank(),
688  "SequenceOfVectors<V,M>::unifiedSampleStd()",
689  "invalid input data");
690 
691  ScalarSequence<double> data(m_env,0,"");
692 
693  unsigned int numParams = this->vectorSizeLocal();
694  for (unsigned int i = 0; i < numParams; ++i) {
695  this->extractScalarSeq(initialPos,
696  1, // spacing
697  numPos,
698  i,
699  data);
700  unifiedStdVec[i] = data.unifiedSampleStd(m_vectorSpace.numOfProcsForStorage() == 1,
701  0,
702  numPos,
703  unifiedMeanVec[i]);
704  }
705 
706  return;
707 }
708 //---------------------------------------------------
709 template <class V, class M>
710 void
712  unsigned int initialPos,
713  unsigned int numPos,
714  const V& meanVec,
715  V& popVec) const
716 {
717  bool bRC = ((initialPos < this->subSequenceSize()) &&
718  (0 < numPos ) &&
719  ((initialPos+numPos) <= this->subSequenceSize()) &&
720  (this->vectorSizeLocal() == meanVec.sizeLocal() ) &&
721  (this->vectorSizeLocal() == popVec.sizeLocal() ));
722  UQ_FATAL_TEST_MACRO(bRC == false,
723  m_env.worldRank(),
724  "SequenceOfVectors<V,M>::subPopulationVariance()",
725  "invalid input data");
726 
727  ScalarSequence<double> data(m_env,0,"");
728 
729  unsigned int numParams = this->vectorSizeLocal();
730  for (unsigned int i = 0; i < numParams; ++i) {
731  this->extractScalarSeq(initialPos,
732  1, // spacing
733  numPos,
734  i,
735  data);
736  popVec[i] = data.subPopulationVariance(0,
737  numPos,
738  meanVec[i]);
739  }
740 
741  return;
742 }
743 //---------------------------------------------------
744 template <class V, class M>
745 void
747  unsigned int initialPos,
748  unsigned int numPos,
749  const V& unifiedMeanVec,
750  V& unifiedPopVec) const
751 {
752  bool bRC = ((initialPos < this->subSequenceSize() ) &&
753  (0 < numPos ) &&
754  ((initialPos+numPos) <= this->subSequenceSize() ) &&
755  (this->vectorSizeLocal() == unifiedMeanVec.sizeLocal()) &&
756  (this->vectorSizeLocal() == unifiedPopVec.sizeLocal() ));
757  UQ_FATAL_TEST_MACRO(bRC == false,
758  m_env.worldRank(),
759  "SequenceOfVectors<V,M>::unifiedPopulationVariance()",
760  "invalid input data");
761 
762  ScalarSequence<double> data(m_env,0,"");
763 
764  unsigned int numParams = this->vectorSizeLocal();
765  for (unsigned int i = 0; i < numParams; ++i) {
766  this->extractScalarSeq(initialPos,
767  1, // spacing
768  numPos,
769  i,
770  data);
771  unifiedPopVec[i] = data.unifiedPopulationVariance(m_vectorSpace.numOfProcsForStorage() == 1,
772  0,
773  numPos,
774  unifiedMeanVec[i]);
775  }
776 
777  return;
778 }
779 //---------------------------------------------------
780 template <class V, class M>
781 void
783  unsigned int initialPos,
784  unsigned int numPos,
785  const V& meanVec,
786  unsigned int lag,
787  V& covVec) const
788 {
789  bool bRC = ((initialPos < this->subSequenceSize()) &&
790  (0 < numPos ) &&
791  ((initialPos+numPos) <= this->subSequenceSize()) &&
792  (this->vectorSizeLocal() == meanVec.sizeLocal() ) &&
793  (lag < numPos ) && // lag should not be too large
794  (this->vectorSizeLocal() == covVec.sizeLocal() ));
795  UQ_FATAL_TEST_MACRO(bRC == false,
796  m_env.worldRank(),
797  "SequenceOfVectors<V,M>::autoCovariance()",
798  "invalid input data");
799 
800  ScalarSequence<double> data(m_env,0,"");
801 
802  unsigned int numParams = this->vectorSizeLocal();
803  for (unsigned int i = 0; i < numParams; ++i) {
804  this->extractScalarSeq(initialPos,
805  1, // spacing
806  numPos,
807  i,
808  data);
809  covVec[i] = data.autoCovariance(0,
810  numPos,
811  meanVec[i],
812  lag);
813  }
814 
815  return;
816 }
817 //---------------------------------------------------
818 template <class V, class M>
819 void
821  unsigned int initialPos,
822  unsigned int numPos,
823  unsigned int lag,
824  V& corrVec) const
825 {
826  bool bRC = ((initialPos < this->subSequenceSize()) &&
827  (0 < numPos ) &&
828  ((initialPos+numPos) <= this->subSequenceSize()) &&
829  (lag < numPos ) && // lag should not be too large
830  (this->vectorSizeLocal() == corrVec.sizeLocal() ));
831  UQ_FATAL_TEST_MACRO(bRC == false,
832  m_env.worldRank(),
833  "SequenceOfVectors<V,M>::autoCorrViaDef()",
834  "invalid input data");
835 
836  ScalarSequence<double> data(m_env,0,"");
837 
838  unsigned int numParams = this->vectorSizeLocal();
839  for (unsigned int i = 0; i < numParams; ++i) {
840  this->extractScalarSeq(initialPos,
841  1, // spacing
842  numPos,
843  i,
844  data);
845  corrVec[i] = data.autoCorrViaDef(0,
846  numPos,
847  lag);
848  }
849 
850  return;
851 }
852 //---------------------------------------------------
853 template <class V, class M>
854 void
856  unsigned int initialPos,
857  unsigned int numPos,
858  const std::vector<unsigned int>& lags,
859  std::vector<V*>& corrVecs) const
860 {
861  bool bRC = ((initialPos < this->subSequenceSize()) &&
862  (0 < numPos ) &&
863  ((initialPos+numPos) <= this->subSequenceSize()) &&
864  (0 < lags.size() ) &&
865  (lags[lags.size()-1] < numPos )); // lag should not be too large
866  UQ_FATAL_TEST_MACRO(bRC == false,
867  m_env.worldRank(),
868  "SequenceOfVectors<V,M>::autoCorrViaFft()",
869  "invalid input data");
870 
871  for (unsigned int j = lags.size(); j < corrVecs.size(); ++j) {
872  if (corrVecs[j] != NULL) {
873  delete corrVecs[j];
874  corrVecs[j] = NULL;
875  }
876  }
877  corrVecs.resize(lags.size(),NULL);
878  for (unsigned int j = 0; j < corrVecs.size(); ++j) {
879  if (corrVecs[j] == NULL) corrVecs[j] = new V(m_vectorSpace.zeroVector());
880  }
881 
882  ScalarSequence<double> data(m_env,0,"");
883  unsigned int maxLag = lags[lags.size()-1];
884  std::vector<double> autoCorrs(maxLag+1,0.); // Yes, +1
885 
886  unsigned int numParams = this->vectorSizeLocal();
887  for (unsigned int i = 0; i < numParams; ++i) {
888  this->extractScalarSeq(initialPos,
889  1, // spacing
890  numPos,
891  i,
892  data);
893 
894  //if (m_env.subDisplayFile()) {
895  // *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::autoCorrViaFft()"
896  // << ": about to call data.autoCorrViaFft() for paramId = " << i
897  // << ", with numPos = " << numPos
898  // << ", maxLag = " << maxLag
899  // << ", autoCorrs.size() = " << autoCorrs.size()
900  // << std::endl;
901  //}
902  data.autoCorrViaFft(0,
903  numPos,
904  maxLag,
905  autoCorrs);
906 
907  for (unsigned int j = 0; j < lags.size(); ++j) {
908  (*(corrVecs[j]))[i] = autoCorrs[lags[j]];
909  }
910  }
911 
912  return;
913 }
914 //---------------------------------------------------
915 template <class V, class M>
916 void
918  unsigned int initialPos,
919  unsigned int numPos,
920  unsigned int numSum,
921  V& autoCorrsSumVec) const
922 {
923  bool bRC = ((initialPos < this->subSequenceSize()) &&
924  (0 < numPos ) &&
925  ((initialPos+numPos) <= this->subSequenceSize()) &&
926  (0 < numSum ) &&
927  (numSum <= numPos ) &&
928  (autoCorrsSumVec.sizeLocal() == this->vectorSizeLocal()));
929  UQ_FATAL_TEST_MACRO(bRC == false,
930  m_env.worldRank(),
931  "SequenceOfVectors<V,M>::autoCorrViaFft(), for sum",
932  "invalid input data");
933 
934  ScalarSequence<double> data(m_env,0,"");
935 
936  unsigned int numParams = this->vectorSizeLocal();
937  for (unsigned int i = 0; i < numParams; ++i) {
938  this->extractScalarSeq(initialPos,
939  1, // spacing
940  numPos,
941  i,
942  data);
943 
944  data.autoCorrViaFft(0,
945  numPos,
946  numSum,
947  autoCorrsSumVec[i]);
948  }
949 
950  return;
951 }
952 //---------------------------------------------------
953 template <class V, class M>
954 void
956  unsigned int initialPos,
957  unsigned int numPos,
958  V& minVec,
959  V& maxVec) const
960 {
961  bool bRC = ((0 < numPos ) &&
962  ((initialPos+numPos) <= this->subSequenceSize()) &&
963  (this->vectorSizeLocal() == minVec.sizeLocal() ) &&
964  (this->vectorSizeLocal() == maxVec.sizeLocal() ));
965  UQ_FATAL_TEST_MACRO(bRC == false,
966  m_env.worldRank(),
967  "SequenceOfVectors<V,M>::subMinMaxExtra()",
968  "invalid input data");
969 
970  //unsigned int numPos = this->subSequenceSize() - initialPos;
971  unsigned int numParams = this->vectorSizeLocal();
972  ScalarSequence<double> data(m_env,0,"");
973 
974  for (unsigned int i = 0; i < numParams; ++i) {
975  this->extractScalarSeq(initialPos,
976  1, // spacing
977  numPos,
978  i,
979  data);
980  data.subMinMaxExtra(0,numPos,minVec[i],maxVec[i]);
981  }
982 
983  return;
984 }
985 //---------------------------------------------------
986 template <class V, class M>
987 void
989  unsigned int initialPos,
990  unsigned int numPos,
991  V& unifiedMinVec,
992  V& unifiedMaxVec) const
993 {
994  bool bRC = ((0 < numPos ) &&
995  ((initialPos+numPos) <= this->subSequenceSize() ) &&
996  (this->vectorSizeLocal() == unifiedMinVec.sizeLocal()) &&
997  (this->vectorSizeLocal() == unifiedMaxVec.sizeLocal()));
998  UQ_FATAL_TEST_MACRO(bRC == false,
999  m_env.worldRank(),
1000  "SequenceOfVectors<V,M>::unifiedMinMaxExtra()",
1001  "invalid input data");
1002 
1003  //unsigned int numPos = this->subSequenceSize() - initialPos;
1004  unsigned int numParams = this->vectorSizeLocal();
1005  ScalarSequence<double> data(m_env,0,"");
1006 
1007  for (unsigned int i = 0; i < numParams; ++i) {
1008  this->extractScalarSeq(initialPos,
1009  1, // spacing
1010  numPos,
1011  i,
1012  data);
1013  data.unifiedMinMaxExtra(m_vectorSpace.numOfProcsForStorage() == 1,
1014  0,
1015  numPos,
1016  unifiedMinVec[i],
1017  unifiedMaxVec[i]);
1018  }
1019 
1020  return;
1021 }
1022 //---------------------------------------------------
1023 template <class V, class M>
1024 void
1026  unsigned int initialPos,
1027  const V& minVec,
1028  const V& maxVec,
1029  std::vector<V*>& centersForAllBins,
1030  std::vector<V*>& quanttsForAllBins) const
1031 {
1032  bool bRC = ((initialPos < this->subSequenceSize() ) &&
1033  (this->vectorSizeLocal() == minVec.sizeLocal() ) &&
1034  (this->vectorSizeLocal() == maxVec.sizeLocal() ) &&
1035  (0 < centersForAllBins.size()) &&
1036  (centersForAllBins.size() == quanttsForAllBins.size()));
1037  UQ_FATAL_TEST_MACRO(bRC == false,
1038  m_env.worldRank(),
1039  "SequenceOfVectors<V,M>::subHistogram()",
1040  "invalid input data");
1041 
1042  for (unsigned int j = 0; j < quanttsForAllBins.size(); ++j) {
1043  centersForAllBins[j] = new V(m_vectorSpace.zeroVector());
1044  quanttsForAllBins [j] = new V(m_vectorSpace.zeroVector());
1045  }
1046 
1047  unsigned int dataSize = this->subSequenceSize() - initialPos;
1048  unsigned int numParams = this->vectorSizeLocal();
1049  for (unsigned int i = 0; i < numParams; ++i) {
1050  ScalarSequence<double> data(m_env,dataSize,"");
1051  for (unsigned int j = 0; j < dataSize; ++j) {
1052  data[j] = (*(m_seq[initialPos+j]))[i];
1053  }
1054 
1055  std::vector<double > centers(centersForAllBins.size(),0.);
1056  std::vector<unsigned int> quantts(quanttsForAllBins.size(), 0 );
1057  data.subHistogram(0,
1058  minVec[i],
1059  maxVec[i],
1060  centers,
1061  quantts);
1062 
1063  for (unsigned int j = 0; j < quantts.size(); ++j) {
1064  (*(centersForAllBins[j]))[i] = centers[j];
1065  (*(quanttsForAllBins[j]))[i] = (double) quantts[j];
1066  }
1067  }
1068 
1069  return;
1070 }
1071 //---------------------------------------------------
1072 template <class V, class M>
1073 void
1075  unsigned int initialPos,
1076  const V& unifiedMinVec,
1077  const V& unifiedMaxVec,
1078  std::vector<V*>& unifiedCentersForAllBins,
1079  std::vector<V*>& unifiedQuanttsForAllBins) const
1080 {
1081  bool bRC = ((initialPos < this->subSequenceSize() ) &&
1082  (this->vectorSizeLocal() == unifiedMinVec.sizeLocal() ) &&
1083  (this->vectorSizeLocal() == unifiedMaxVec.sizeLocal() ) &&
1084  (0 < unifiedCentersForAllBins.size()) &&
1085  (unifiedCentersForAllBins.size() == unifiedQuanttsForAllBins.size()));
1086  UQ_FATAL_TEST_MACRO(bRC == false,
1087  m_env.worldRank(),
1088  "SequenceOfVectors<V,M>::unifiedHistogram()",
1089  "invalid input data");
1090 
1091  for (unsigned int j = 0; j < unifiedQuanttsForAllBins.size(); ++j) {
1092  unifiedCentersForAllBins[j] = new V(m_vectorSpace.zeroVector());
1093  unifiedQuanttsForAllBins [j] = new V(m_vectorSpace.zeroVector());
1094  }
1095 
1096  unsigned int dataSize = this->subSequenceSize() - initialPos;
1097  unsigned int numParams = this->vectorSizeLocal();
1098  for (unsigned int i = 0; i < numParams; ++i) {
1099  ScalarSequence<double> data(m_env,dataSize,"");
1100  for (unsigned int j = 0; j < dataSize; ++j) {
1101  data[j] = (*(m_seq[initialPos+j]))[i];
1102  }
1103 
1104  std::vector<double > unifiedCenters(unifiedCentersForAllBins.size(),0.);
1105  std::vector<unsigned int> unifiedQuantts(unifiedQuanttsForAllBins.size(), 0 );
1106  data.unifiedHistogram(m_vectorSpace.numOfProcsForStorage() == 1,
1107  0,
1108  unifiedMinVec[i],
1109  unifiedMaxVec[i],
1110  unifiedCenters,
1111  unifiedQuantts);
1112 
1113  for (unsigned int j = 0; j < unifiedQuantts.size(); ++j) {
1114  (*(unifiedCentersForAllBins[j]))[i] = unifiedCenters[j];
1115  (*(unifiedQuanttsForAllBins[j]))[i] = (double) unifiedQuantts[j];
1116  }
1117  }
1118 
1119  return;
1120 }
1121 //---------------------------------------------------
1122 template <class V, class M>
1123 void
1125  unsigned int initialPos,
1126  V& iqrVec) const
1127 {
1128  bool bRC = ((initialPos < this->subSequenceSize()) &&
1129  (this->vectorSizeLocal() == iqrVec.sizeLocal() ));
1130  UQ_FATAL_TEST_MACRO(bRC == false,
1131  m_env.worldRank(),
1132  "SequenceOfVectors<V,M>::subInterQuantileRange()",
1133  "invalid input data");
1134 
1135  unsigned int numPos = this->subSequenceSize() - initialPos;
1136  ScalarSequence<double> data(m_env,0,"");
1137 
1138  unsigned int numParams = this->vectorSizeLocal();
1139  for (unsigned int i = 0; i < numParams; ++i) {
1140  this->extractScalarSeq(initialPos,
1141  1, // spacing
1142  numPos,
1143  i,
1144  data);
1145  iqrVec[i] = data.subInterQuantileRange(0);
1146  }
1147 
1148  return;
1149 }
1150 //---------------------------------------------------
1151 template <class V, class M>
1152 void
1154  unsigned int initialPos,
1155  V& unifiedIqrVec) const
1156 {
1157  bool bRC = ((initialPos < this->subSequenceSize() ) &&
1158  (this->vectorSizeLocal() == unifiedIqrVec.sizeLocal()));
1159  UQ_FATAL_TEST_MACRO(bRC == false,
1160  m_env.worldRank(),
1161  "SequenceOfVectors<V,M>::unifiedInterQuantileRange()",
1162  "invalid input data");
1163 
1164  unsigned int numPos = this->subSequenceSize() - initialPos;
1165  ScalarSequence<double> data(m_env,0,"");
1166 
1167  unsigned int numParams = this->vectorSizeLocal();
1168  for (unsigned int i = 0; i < numParams; ++i) {
1169  this->extractScalarSeq(initialPos,
1170  1, // spacing
1171  numPos,
1172  i,
1173  data);
1174  unifiedIqrVec[i] = data.unifiedInterQuantileRange(m_vectorSpace.numOfProcsForStorage() == 1,
1175  0);
1176  }
1177 
1178  return;
1179 }
1180 //---------------------------------------------------
1181 template <class V, class M>
1182 void
1184  unsigned int initialPos,
1185  const V& iqrVec,
1186  unsigned int kdeDimension,
1187  V& scaleVec) const
1188 {
1189  bool bRC = ((initialPos < this->subSequenceSize()) &&
1190  (this->vectorSizeLocal() == iqrVec.sizeLocal() ) &&
1191  (this->vectorSizeLocal() == scaleVec.sizeLocal() ));
1192  UQ_FATAL_TEST_MACRO(bRC == false,
1193  m_env.worldRank(),
1194  "SequenceOfVectors<V,M>::subScalesForKde()",
1195  "invalid input data");
1196 
1197  unsigned int numPos = this->subSequenceSize() - initialPos;
1198  ScalarSequence<double> data(m_env,0,"");
1199 
1200  unsigned int numParams = this->vectorSizeLocal();
1201  for (unsigned int i = 0; i < numParams; ++i) {
1202  this->extractScalarSeq(initialPos,
1203  1, // spacing
1204  numPos,
1205  i,
1206  data);
1207  scaleVec[i] = data.subScaleForKde(0,
1208  iqrVec[i],
1209  kdeDimension);
1210  }
1211 
1212  return;
1213 }
1214 //---------------------------------------------------
1215 template <class V, class M>
1216 void
1218  unsigned int initialPos,
1219  const V& unifiedIqrVec,
1220  unsigned int kdeDimension,
1221  V& unifiedScaleVec) const
1222 {
1223  bool bRC = ((initialPos < this->subSequenceSize() ) &&
1224  (this->vectorSizeLocal() == unifiedIqrVec.sizeLocal() ) &&
1225  (this->vectorSizeLocal() == unifiedScaleVec.sizeLocal()));
1226  UQ_FATAL_TEST_MACRO(bRC == false,
1227  m_env.worldRank(),
1228  "SequenceOfVectors<V,M>::unifiedScalesForKde()",
1229  "invalid input data");
1230 
1231  unsigned int numPos = this->subSequenceSize() - initialPos;
1232  ScalarSequence<double> data(m_env,0,"");
1233 
1234  unsigned int numParams = this->vectorSizeLocal();
1235  for (unsigned int i = 0; i < numParams; ++i) {
1236  this->extractScalarSeq(initialPos,
1237  1, // spacing
1238  numPos,
1239  i,
1240  data);
1241  unifiedScaleVec[i] = data.unifiedScaleForKde(m_vectorSpace.numOfProcsForStorage() == 1,
1242  0,
1243  unifiedIqrVec[i],
1244  kdeDimension);
1245  }
1246 
1247  return;
1248 }
1249 //---------------------------------------------------
1250 template <class V, class M>
1251 void
1253  unsigned int initialPos,
1254  const V& scaleVec,
1255  const std::vector<V*>& evalParamVecs,
1256  std::vector<V*>& densityVecs) const
1257 {
1258  bool bRC = ((initialPos < this->subSequenceSize()) &&
1259  (this->vectorSizeLocal() == scaleVec.sizeLocal() ) &&
1260  (0 < evalParamVecs.size() ) &&
1261  (evalParamVecs.size() == densityVecs.size() ));
1262  UQ_FATAL_TEST_MACRO(bRC == false,
1263  m_env.worldRank(),
1264  "SequenceOfVectors<V,M>::subGaussian1dKde()",
1265  "invalid input data");
1266 
1267  unsigned int numPos = this->subSequenceSize() - initialPos;
1268  ScalarSequence<double> data(m_env,0,"");
1269 
1270  unsigned int numEvals = evalParamVecs.size();
1271  for (unsigned int j = 0; j < numEvals; ++j) {
1272  densityVecs[j] = new V(m_vectorSpace.zeroVector());
1273  }
1274  std::vector<double> evalParams(numEvals,0.);
1275  std::vector<double> densities (numEvals,0.);
1276 
1277  unsigned int numParams = this->vectorSizeLocal();
1278  for (unsigned int i = 0; i < numParams; ++i) {
1279  this->extractScalarSeq(initialPos,
1280  1, // spacing
1281  numPos,
1282  i,
1283  data);
1284 
1285  for (unsigned int j = 0; j < numEvals; ++j) {
1286  evalParams[j] = (*evalParamVecs[j])[i];
1287  }
1288 
1289  data.subGaussian1dKde(0,
1290  scaleVec[i],
1291  evalParams,
1292  densities);
1293 
1294  for (unsigned int j = 0; j < numEvals; ++j) {
1295  (*densityVecs[j])[i] = densities[j];
1296  }
1297  }
1298 
1299  return;
1300 }
1301 //---------------------------------------------------
1302 template <class V, class M>
1303 void
1305  unsigned int initialPos,
1306  const V& unifiedScaleVec,
1307  const std::vector<V*>& unifiedEvalParamVecs,
1308  std::vector<V*>& unifiedDensityVecs) const
1309 {
1310  bool bRC = ((initialPos < this->subSequenceSize() ) &&
1311  (this->vectorSizeLocal() == unifiedScaleVec.sizeLocal()) &&
1312  (0 < unifiedEvalParamVecs.size()) &&
1313  (unifiedEvalParamVecs.size() == unifiedDensityVecs.size() ));
1314  UQ_FATAL_TEST_MACRO(bRC == false,
1315  m_env.worldRank(),
1316  "SequenceOfVectors<V,M>::unifiedGaussian1dKde()",
1317  "invalid input data");
1318 
1319  unsigned int numPos = this->subSequenceSize() - initialPos;
1320  ScalarSequence<double> data(m_env,0,"");
1321 
1322  unsigned int numEvals = unifiedEvalParamVecs.size();
1323  for (unsigned int j = 0; j < numEvals; ++j) {
1324  unifiedDensityVecs[j] = new V(m_vectorSpace.zeroVector());
1325  }
1326  std::vector<double> unifiedEvalParams(numEvals,0.);
1327  std::vector<double> unifiedDensities (numEvals,0.);
1328 
1329  unsigned int numParams = this->vectorSizeLocal();
1330  for (unsigned int i = 0; i < numParams; ++i) {
1331  this->extractScalarSeq(initialPos,
1332  1, // spacing
1333  numPos,
1334  i,
1335  data);
1336 
1337  for (unsigned int j = 0; j < numEvals; ++j) {
1338  unifiedEvalParams[j] = (*unifiedEvalParamVecs[j])[i];
1339  }
1340 
1341  data.unifiedGaussian1dKde(m_vectorSpace.numOfProcsForStorage() == 1,
1342  0,
1343  unifiedScaleVec[i],
1344  unifiedEvalParams,
1345  unifiedDensities);
1346 
1347  for (unsigned int j = 0; j < numEvals; ++j) {
1348  (*unifiedDensityVecs[j])[i] = unifiedDensities[j];
1349  }
1350  }
1351 
1352  return;
1353 }
1354 //---------------------------------------------------
1355 template <class V, class M>
1356 void
1358  unsigned int initialPos,
1359  unsigned int numPos,
1360  const std::string& fileName,
1361  const std::string& fileType,
1362  const std::set<unsigned int>& allowedSubEnvIds) const
1363 {
1364  UQ_FATAL_TEST_MACRO(m_env.subRank() < 0,
1365  m_env.worldRank(),
1366  "SequenceOfVectors<V,M>::subWriteContents(1)",
1367  "unexpected subRank");
1368 
1369  FilePtrSetStruct filePtrSet;
1370  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1371  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::subWriteContents()"
1372  << ": about to try to open file '" << fileName << "." << fileType
1373  << "'"
1374  << ", initialPos = " << initialPos
1375  << ", numPos = " << numPos
1376  << std::endl;
1377  }
1378  if (m_env.openOutputFile(fileName,
1379  fileType,
1380  allowedSubEnvIds,
1381  false, // A 'true' causes problems when the user chooses (via options
1382  // in the input file) to use just one file for all outputs.
1383  filePtrSet)) {
1384  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1385  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::subWriteContents()"
1386  << ": successfully opened file '" << fileName << "." << fileType
1387  << "'"
1388  << std::endl;
1389  }
1390  this->subWriteContents(initialPos,
1391  numPos,
1392  filePtrSet,
1393  fileType);
1394  m_env.closeFile(filePtrSet,fileType);
1395  }
1396  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1397  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::subWriteContents()"
1398  << ": before Barrier()"
1399  << std::endl;
1400  }
1401  m_env.subComm().Barrier();
1402 
1403  return;
1404 }
1405 //---------------------------------------------------
1406 template <class V, class M>
1407 void
1409  unsigned int initialPos,
1410  unsigned int numPos,
1411  FilePtrSetStruct& filePtrSet,
1412  const std::string& fileType) const // "m or hdf"
1413 {
1414  UQ_FATAL_TEST_MACRO(filePtrSet.ofsVar == NULL,
1415  m_env.worldRank(),
1416  "SequenceOfVectors<V,M>::subWriteContents(2)",
1417  "filePtrSet.ofsVar should not be NULL");
1418 
1419  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1420  this->subWriteContents(initialPos,
1421  numPos,
1422  *filePtrSet.ofsVar,
1423  fileType);
1424  }
1425  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1426  UQ_FATAL_TEST_MACRO(true,
1427  m_env.worldRank(),
1428  "SequenceOfVectors<V,M>::subWriteContents(2)",
1429  "hdf file type not supported yet");
1430  }
1431  else {
1432  UQ_FATAL_TEST_MACRO(true,
1433  m_env.worldRank(),
1434  "SequenceOfVectors<V,M>::subWriteContents(2)",
1435  "invalid file type");
1436  }
1437 
1438  return;
1439 }
1440 //---------------------------------------------------
1441 template <class V, class M>
1442 void
1444  unsigned int initialPos,
1445  unsigned int numPos,
1446  std::ofstream& ofs,
1447  const std::string& fileType) const // "m or hdf"
1448 {
1449  UQ_FATAL_TEST_MACRO((initialPos+numPos) > this->subSequenceSize(),
1450  m_env.worldRank(),
1451  "SequenceOfVectors<V,M>::subWriteContents(3)",
1452  "invalid routine input parameters");
1453 
1454  if (fileType.c_str()) {}; // just to avoid compiler warning
1455 
1456  if (initialPos == 0) {
1457  ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << this->subSequenceSize()
1458  << "," << this->vectorSizeLocal()
1459  << ");"
1460  << std::endl;
1461  ofs << m_name << "_sub" << m_env.subIdString() << " = [";
1462  }
1463 
1464  for (unsigned int j = initialPos; j < initialPos+numPos; ++j) {
1465  bool savedVectorPrintScientific = m_seq[j]->getPrintScientific();
1466  bool savedVectorPrintState = m_seq[j]->getPrintHorizontally();
1467  m_seq[j]->setPrintScientific (true);
1468  m_seq[j]->setPrintHorizontally(true);
1469 
1470  ofs << *(m_seq[j])
1471  << std::endl;
1472 
1473  m_seq[j]->setPrintHorizontally(savedVectorPrintState);
1474  m_seq[j]->setPrintScientific (savedVectorPrintScientific);
1475  }
1476  if ((initialPos+numPos) == this->subSequenceSize()) {
1477  ofs << "];\n";
1478  }
1479 
1480  return;
1481 }
1482 //---------------------------------------------------
1483 template <class V, class M>
1484 void
1486  const std::string& fileName,
1487  const std::string& inputFileType) const
1488 {
1489  std::string fileType(inputFileType);
1490 #ifdef QUESO_HAS_HDF5
1491  // Do nothing
1492 #else
1493  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1494  if (m_env.subDisplayFile()) {
1495  *m_env.subDisplayFile() << "WARNING in SequenceOfVectors<V,M>::unifiedWriteContents()"
1496  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1497  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1498  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1499  << "' instead..."
1500  << std::endl;
1501  }
1502  if (m_env.subRank() == 0) {
1503  std::cerr << "WARNING in SequenceOfVectors<V,M>::unifiedWriteContents()"
1504  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1505  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1506  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1507  << "' instead..."
1508  << std::endl;
1509  }
1511  }
1512 #endif
1513 
1514  // All processors in 'fullComm' should call this routine...
1515 
1516  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ... // prudenci-2011-01-17
1517  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1518  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedWriteContents()"
1519  << ": worldRank " << m_env.worldRank()
1520  << ", fullRank " << m_env.fullRank()
1521  << ", subEnvironment " << m_env.subId()
1522  << ", subRank " << m_env.subRank()
1523  << ", inter0Rank " << m_env.inter0Rank()
1524  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1525  << ", fileName = " << fileName
1526  << std::endl;
1527  }
1528 
1529  if (m_env.inter0Rank() >= 0) {
1530  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
1531  if (m_env.inter0Rank() == (int) r) {
1532  // My turn
1533  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1534  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1535  << ": worldRank " << m_env.worldRank()
1536  << ", fullRank " << m_env.fullRank()
1537  << ", subEnvironment " << m_env.subId()
1538  << ", subRank " << m_env.subRank()
1539  << ", inter0Rank " << m_env.inter0Rank()
1540  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1541  << ", fileName = " << fileName
1542  << ", about to open file for r = " << r
1543  << std::endl;
1544  }
1545 
1546  // bool writeOver = (r == 0);
1547  bool writeOver = false; // A 'true' causes problems when the user chooses (via options
1548  // in the input file) to use just one file for all outputs.
1549  FilePtrSetStruct unifiedFilePtrSet;
1550  if (m_env.openUnifiedOutputFile(fileName,
1551  fileType, // "m or hdf"
1552  writeOver,
1553  unifiedFilePtrSet)) {
1554  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) { // 2013-02-23
1555  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1556  << ": worldRank " << m_env.worldRank()
1557  << ", fullRank " << m_env.fullRank()
1558  << ", subEnvironment " << m_env.subId()
1559  << ", subRank " << m_env.subRank()
1560  << ", inter0Rank " << m_env.inter0Rank()
1561  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1562  << ", fileName = " << fileName
1563  << ", just opened file for r = " << r
1564  << std::endl;
1565  }
1566 
1567  unsigned int chainSize = this->subSequenceSize();
1568  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1569  if (r == 0) {
1570  *unifiedFilePtrSet.ofsVar << m_name << "_unified" << " = zeros(" << this->subSequenceSize()*m_env.inter0Comm().NumProc()
1571  << "," << this->vectorSizeLocal()
1572  << ");"
1573  << std::endl;
1574  *unifiedFilePtrSet.ofsVar << m_name << "_unified" << " = [";
1575  }
1576 
1577  for (unsigned int j = 0; j < chainSize; ++j) { // 2013-02-23
1578  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): m_seq[" << j << "] = " << m_seq[j]
1579  // << std::endl;
1580  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): &(m_seq[" << j << "].map()) = " << &(m_seq[j]->map())
1581  // << std::endl;
1582  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): (m_seq[" << j << "].map().NumMyElements = " << m_seq[j]->map().NumMyElements()
1583  // << std::endl;
1584  //V tmpVec(*(m_seq[j]));
1585  //std::cout << "*(m_seq[" << j << "]) = " << tmpVec
1586  // << std::endl;
1587  //std::cout << "*(m_seq[" << j << "]) = " << *(m_seq[j])
1588  // << std::endl;
1589 
1590  bool savedVectorPrintScientific = m_seq[j]->getPrintScientific();
1591  bool savedVectorPrintState = m_seq[j]->getPrintHorizontally();
1592  m_seq[j]->setPrintScientific (true);
1593  m_seq[j]->setPrintHorizontally(true);
1594 
1595  *unifiedFilePtrSet.ofsVar << *(m_seq[j])
1596  << std::endl;
1597 
1598  m_seq[j]->setPrintHorizontally(savedVectorPrintState);
1599  m_seq[j]->setPrintScientific (savedVectorPrintScientific);
1600  }
1601  }
1602 #ifdef QUESO_HAS_HDF5
1603  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1604  unsigned int numParams = m_vectorSpace.dimLocal();
1605  if (r == 0) {
1606  hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
1607  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data type created" << std::endl;
1608  hsize_t dimsf[2];
1609  dimsf[0] = numParams;
1610  dimsf[1] = chainSize;
1611  hid_t dataspace = H5Screate_simple(2, dimsf, NULL); // HDF5_rank = 2
1612  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data space created" << std::endl;
1613  hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var,
1614  "seq_of_vectors",
1615  datatype,
1616  dataspace,
1617  H5P_DEFAULT, // Link creation property list
1618  H5P_DEFAULT, // Dataset creation property list
1619  H5P_DEFAULT); // Dataset access property list
1620  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data set created" << std::endl;
1621 
1622  struct timeval timevalBegin;
1623  int iRC = UQ_OK_RC;
1624  iRC = gettimeofday(&timevalBegin,NULL);
1625  if (iRC) {}; // just to remover compiler warning
1626 
1627  //double* dataOut[numParams]; // avoid compiler warning
1628  std::vector<double*> dataOut((size_t) numParams,NULL);
1629  dataOut[0] = (double*) malloc(numParams*chainSize*sizeof(double));
1630  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
1631  dataOut[i] = dataOut[i-1] + chainSize; // Yes, just 'chainSize', not 'chainSize*sizeof(double)'
1632  }
1633  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, memory allocated" << std::endl;
1634  for (unsigned int j = 0; j < chainSize; ++j) {
1635  V tmpVec(*(m_seq[j]));
1636  for (unsigned int i = 0; i < numParams; ++i) {
1637  dataOut[i][j] = tmpVec[i];
1638  }
1639  }
1640  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, memory filled" << std::endl;
1641 
1642  herr_t status;
1643  status = H5Dwrite(dataset,
1644  H5T_NATIVE_DOUBLE,
1645  H5S_ALL,
1646  H5S_ALL,
1647  H5P_DEFAULT,
1648  (void*) dataOut[0]);
1649  if (status) {}; // just to remover compiler warning
1650  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data written" << std::endl;
1651 
1652  double writeTime = MiscGetEllapsedSeconds(&timevalBegin);
1653  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
1654  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1655  << ": worldRank " << m_env.worldRank()
1656  << ", fullRank " << m_env.fullRank()
1657  << ", subEnvironment " << m_env.subId()
1658  << ", subRank " << m_env.subRank()
1659  << ", inter0Rank " << m_env.inter0Rank()
1660  << ", fileName = " << fileName
1661  << ", numParams = " << numParams
1662  << ", chainSize = " << chainSize
1663  << ", writeTime = " << writeTime << " seconds"
1664  << std::endl;
1665  }
1666 
1667  H5Dclose(dataset);
1668  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data set closed" << std::endl;
1669  H5Sclose(dataspace);
1670  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data space closed" << std::endl;
1671  H5Tclose(datatype);
1672  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data type closed" << std::endl;
1673  //free(dataOut[0]); // related to the changes above for compiler warning
1674  for (unsigned int tmpIndex = 0; tmpIndex < dataOut.size(); tmpIndex++) {
1675  free (dataOut[tmpIndex]);
1676  }
1677  }
1678  else {
1679  UQ_FATAL_TEST_MACRO(true,
1680  m_env.worldRank(),
1681  "SequenceOfVectors<V,M>::unifiedWriteContents()",
1682  "hdf file type not supported for multiple sub-environments yet");
1683  }
1684  }
1685 #endif
1686  else {
1687  UQ_FATAL_TEST_MACRO(true,
1688  m_env.worldRank(),
1689  "SequenceOfVectors<V,M>::unifiedWriteContents()",
1690  "invalid file type");
1691  }
1692 
1693  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1694  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1695  << ": worldRank " << m_env.worldRank()
1696  << ", fullRank " << m_env.fullRank()
1697  << ", subEnvironment " << m_env.subId()
1698  << ", subRank " << m_env.subRank()
1699  << ", inter0Rank " << m_env.inter0Rank()
1700  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1701  << ", fileName = " << fileName
1702  << ", about to close file for r = " << r
1703  << std::endl;
1704  }
1705 
1706  m_env.closeFile(unifiedFilePtrSet,fileType);
1707 
1708  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1709  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1710  << ": worldRank " << m_env.worldRank()
1711  << ", fullRank " << m_env.fullRank()
1712  << ", subEnvironment " << m_env.subId()
1713  << ", subRank " << m_env.subRank()
1714  << ", inter0Rank " << m_env.inter0Rank()
1715  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1716  << ", fileName = " << fileName
1717  << ", just closed file for r = " << r
1718  << std::endl;
1719  }
1720  } // if (m_env.openUnifiedOutputFile())
1721  } // if (m_env.inter0Rank() == (int) r)
1722  m_env.inter0Comm().Barrier();
1723  } // for r
1724 
1725  if (m_env.inter0Rank() == 0) {
1726  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1727  FilePtrSetStruct unifiedFilePtrSet;
1728  if (m_env.openUnifiedOutputFile(fileName,
1729  fileType,
1730  false, // Yes, 'writeOver = false' in order to close the array for matlab
1731  unifiedFilePtrSet)) {
1732  *unifiedFilePtrSet.ofsVar << "];\n";
1733  m_env.closeFile(unifiedFilePtrSet,fileType);
1734  }
1735  }
1736  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1737  // Do nothing
1738  }
1739  else {
1740  UQ_FATAL_TEST_MACRO(true,
1741  m_env.worldRank(),
1742  "SequenceOfVectors<V,M>::unifiedWriteContents(), final",
1743  "invalid file type");
1744  }
1745  }
1746  } // if (m_env.inter0Rank() >= 0)
1747 
1748  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1749  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedWriteContents()"
1750  << ", fileName = " << fileName
1751  << std::endl;
1752  }
1753  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ... // prudenci-2011-01-17
1754 
1755  return;
1756 }
1757 //---------------------------------------------------
1758 template <class V, class M>
1759 void
1761  const std::string& fileName,
1762  const std::string& inputFileType,
1763  const unsigned int subReadSize)
1764 {
1765  std::string fileType(inputFileType);
1766 #ifdef QUESO_HAS_HDF5
1767  // Do nothing
1768 #else
1769  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1770  if (m_env.subDisplayFile()) {
1771  *m_env.subDisplayFile() << "WARNING in SequenceOfVectors<V,M>::unifiedReadContents()"
1772  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1773  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1774  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1775  << "' instead..."
1776  << std::endl;
1777  }
1778  if (m_env.subRank() == 0) {
1779  std::cerr << "WARNING in SequenceOfVectors<V,M>::unifiedReadContents()"
1780  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1781  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1782  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1783  << "' instead..."
1784  << std::endl;
1785  }
1787  }
1788 #endif
1789 
1790  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1791  if (m_env.subDisplayFile()) {
1792  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedReadContents()"
1793  << ": worldRank " << m_env.worldRank()
1794  << ", fullRank " << m_env.fullRank()
1795  << ", subEnvironment " << m_env.subId()
1796  << ", subRank " << m_env.subRank()
1797  << ", inter0Rank " << m_env.inter0Rank()
1798  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1799  << ", fileName = " << fileName
1800  << ", subReadSize = " << subReadSize
1801  //<< ", unifiedReadSize = " << unifiedReadSize
1802  << std::endl;
1803  }
1804 
1805  this->resizeSequence(subReadSize);
1806 
1807  if (m_env.inter0Rank() >= 0) {
1808  double unifiedReadSize = subReadSize*m_env.inter0Comm().NumProc();
1809 
1810  // In the logic below, the id of a line' begins with value 0 (zero)
1811  unsigned int idOfMyFirstLine = 1 + m_env.inter0Rank()*subReadSize;
1812  unsigned int idOfMyLastLine = (1 + m_env.inter0Rank())*subReadSize;
1813  unsigned int numParams = this->vectorSizeLocal();
1814 
1815  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // "m or hdf"
1816  if (m_env.inter0Rank() == (int) r) {
1817  // My turn
1818  FilePtrSetStruct unifiedFilePtrSet;
1819  if (m_env.openUnifiedInputFile(fileName,
1820  fileType,
1821  unifiedFilePtrSet)) {
1822  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1823  if (r == 0) {
1824  // Read number of chain positions in the file by taking care of the first line,
1825  // which resembles something like 'variable_name = zeros(n_positions,m_params);'
1826  std::string tmpString;
1827 
1828  // Read 'variable name' string
1829  *unifiedFilePtrSet.ifsVar >> tmpString;
1830  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1831 
1832  // Read '=' sign
1833  *unifiedFilePtrSet.ifsVar >> tmpString;
1834  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1835  UQ_FATAL_TEST_MACRO(tmpString != "=",
1836  m_env.worldRank(),
1837  "SequenceOfVectors<V,M>::unifiedReadContents()",
1838  "string should be the '=' sign");
1839 
1840  // Read 'zeros(n_positions,n_params)' string
1841  *unifiedFilePtrSet.ifsVar >> tmpString;
1842  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1843  unsigned int posInTmpString = 6;
1844 
1845  // Isolate 'n_positions' in a string
1846  //char nPositionsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
1847  std::string nPositionsString((size_t) (tmpString.size()-posInTmpString+1),' ');
1848  unsigned int posInPositionsString = 0;
1849  do {
1850  UQ_FATAL_TEST_MACRO(posInTmpString >= tmpString.size(),
1851  m_env.worldRank(),
1852  "SequenceOfVectors<V,M>::unifiedReadContents()",
1853  "symbol ',' not found in first line of file");
1854  nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
1855  } while (tmpString[posInTmpString] != ',');
1856  nPositionsString[posInPositionsString] = '\0';
1857 
1858  // Isolate 'n_params' in a string
1859  posInTmpString++; // Avoid reading ',' char
1860  //char nParamsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
1861  std::string nParamsString((size_t) (tmpString.size()-posInTmpString+1),' ');
1862  unsigned int posInParamsString = 0;
1863  do {
1864  UQ_FATAL_TEST_MACRO(posInTmpString >= tmpString.size(),
1865  m_env.worldRank(),
1866  "SequenceOfVectors<V,M>::unifiedReadContents()",
1867  "symbol ')' not found in first line of file");
1868  nParamsString[posInParamsString++] = tmpString[posInTmpString++];
1869  } while (tmpString[posInTmpString] != ')');
1870  nParamsString[posInParamsString] = '\0';
1871 
1872  // Convert 'n_positions' and 'n_params' strings to numbers
1873  unsigned int sizeOfChainInFile = (unsigned int) strtod(nPositionsString.c_str(),NULL);
1874  unsigned int numParamsInFile = (unsigned int) strtod(nParamsString.c_str(), NULL);
1875  if (m_env.subDisplayFile()) {
1876  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedReadContents()"
1877  << ": worldRank " << m_env.worldRank()
1878  << ", fullRank " << m_env.fullRank()
1879  << ", sizeOfChainInFile = " << sizeOfChainInFile
1880  << ", numParamsInFile = " << numParamsInFile
1881  << std::endl;
1882  }
1883 
1884  // Check if [size of chain in file] >= [requested unified sequence size]
1885  UQ_FATAL_TEST_MACRO(sizeOfChainInFile < unifiedReadSize,
1886  m_env.worldRank(),
1887  "SequenceOfVectors<V,M>::unifiedReadContents()",
1888  "size of chain in file is not big enough");
1889 
1890  // Check if [num params in file] == [num params in current chain]
1891  UQ_FATAL_TEST_MACRO(numParamsInFile != numParams,
1892  m_env.worldRank(),
1893  "SequenceOfVectors<V,M>::unifiedReadContents()",
1894  "number of parameters of chain in file is different than number of parameters in this chain object");
1895  } // if (r == 0)
1896 
1897  // Code common to any core in 'inter0Comm', including core of rank 0
1898  unsigned int maxCharsPerLine = 64*numParams; // Up to about 60 characters to represent each parameter value
1899 
1900  unsigned int lineId = 0;
1901  while (lineId < idOfMyFirstLine) {
1902  unifiedFilePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
1903  lineId++;
1904  };
1905 
1906  if (r == 0) {
1907  // Take care of initial part of the first data line,
1908  // which resembles something like 'variable_name = [value1 value2 ...'
1909  std::string tmpString;
1910 
1911  // Read 'variable name' string
1912  *unifiedFilePtrSet.ifsVar >> tmpString;
1913  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1914 
1915  // Read '=' sign
1916  *unifiedFilePtrSet.ifsVar >> tmpString;
1917  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1918  UQ_FATAL_TEST_MACRO(tmpString != "=",
1919  m_env.worldRank(),
1920  "SequenceOfVectors<V,M>::unifiedReadContents()",
1921  "in core 0, string should be the '=' sign");
1922 
1923  // Take into account the ' [' portion
1924  std::streampos tmpPos = unifiedFilePtrSet.ifsVar->tellg();
1925  unifiedFilePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
1926  }
1927 
1928  V tmpVec(m_vectorSpace.zeroVector());
1929  while (lineId <= idOfMyLastLine) {
1930  for (unsigned int i = 0; i < numParams; ++i) {
1931  *unifiedFilePtrSet.ifsVar >> tmpVec[i];
1932  }
1933  this->setPositionValues(lineId - idOfMyFirstLine, tmpVec);
1934  lineId++;
1935  };
1936  }
1937 #ifdef QUESO_HAS_HDF5
1938  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1939  if (r == 0) {
1940  hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var,
1941  "seq_of_vectors",
1942  H5P_DEFAULT); // Dataset access property list
1943  hid_t datatype = H5Dget_type(dataset);
1944  H5T_class_t t_class = H5Tget_class(datatype);
1945  UQ_FATAL_TEST_MACRO(t_class != H5T_FLOAT,
1946  m_env.worldRank(),
1947  "SequenceOfVectors<V,M>::unifiedReadContents()",
1948  "t_class is not H5T_DOUBLE");
1949  hid_t dataspace = H5Dget_space(dataset);
1950  int rank = H5Sget_simple_extent_ndims(dataspace);
1951  UQ_FATAL_TEST_MACRO(rank != 2,
1952  m_env.worldRank(),
1953  "SequenceOfVectors<V,M>::unifiedReadContents()",
1954  "hdf rank is not 2");
1955  hsize_t dims_in[2];
1956  int status_n;
1957  status_n = H5Sget_simple_extent_dims(dataspace, dims_in, NULL);
1958  if (status_n) {}; // just to remover compiler warning
1959  //std::cout << "In SequenceOfVectors<V,M>::unifiedReadContents()"
1960  // << ": dims_in[0] = " << dims_in[0]
1961  // << ", dims_in[1] = " << dims_in[1]
1962  // << std::endl;
1963  UQ_FATAL_TEST_MACRO(dims_in[0] != numParams,
1964  m_env.worldRank(),
1965  "SequenceOfVectors<V,M>::unifiedReadContents()",
1966  "dims_in[0] is not equal to 'numParams'");
1967  UQ_FATAL_TEST_MACRO(dims_in[1] < subReadSize,
1968  m_env.worldRank(),
1969  "SequenceOfVectors<V,M>::unifiedReadContents()",
1970  "dims_in[1] is smaller that requested 'subReadSize'");
1971 
1972  struct timeval timevalBegin;
1973  int iRC = UQ_OK_RC;
1974  iRC = gettimeofday(&timevalBegin,NULL);
1975  if (iRC) {}; // just to remover compiler warning
1976 
1977  unsigned int chainSizeIn = (unsigned int) dims_in[1];
1978  //double* dataIn[numParams]; // avoid compiler warning
1979  std::vector<double*> dataIn((size_t) numParams,NULL);
1980  dataIn[0] = (double*) malloc(numParams*chainSizeIn*sizeof(double));
1981  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
1982  dataIn[i] = dataIn[i-1] + chainSizeIn; // Yes, just 'chainSizeIn', not 'chainSizeIn*sizeof(double)'
1983  }
1984  //std::cout << "In SequenceOfVectors<V,M>::unifiedReadContents(): h5 case, memory allocated" << std::endl;
1985  herr_t status;
1986  status = H5Dread(dataset,
1987  H5T_NATIVE_DOUBLE,
1988  H5S_ALL,
1989  dataspace,
1990  H5P_DEFAULT,
1991  dataIn[0]);
1992  if (status) {}; // just to remover compiler warning
1993  //std::cout << "In SequenceOfVectors<V,M>::unifiedReadContents(): h5 case, data read" << std::endl;
1994  V tmpVec(m_vectorSpace.zeroVector());
1995  for (unsigned int j = 0; j < subReadSize; ++j) { // Yes, 'subReadSize', not 'chainSizeIn'
1996  for (unsigned int i = 0; i < numParams; ++i) {
1997  tmpVec[i] = dataIn[i][j];
1998  }
1999  this->setPositionValues(j, tmpVec);
2000  }
2001 
2002  double readTime = MiscGetEllapsedSeconds(&timevalBegin);
2003  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
2004  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedReadContents()"
2005  << ": worldRank " << m_env.worldRank()
2006  << ", fullRank " << m_env.fullRank()
2007  << ", subEnvironment " << m_env.subId()
2008  << ", subRank " << m_env.subRank()
2009  << ", inter0Rank " << m_env.inter0Rank()
2010  << ", fileName = " << fileName
2011  << ", numParams = " << numParams
2012  << ", chainSizeIn = " << chainSizeIn
2013  << ", subReadSize = " << subReadSize
2014  << ", readTime = " << readTime << " seconds"
2015  << std::endl;
2016  }
2017 
2018  H5Sclose(dataspace);
2019  H5Tclose(datatype);
2020  H5Dclose(dataset);
2021  //free(dataIn[0]); // related to the changes above for compiler warning
2022  for (unsigned int tmpIndex = 0; tmpIndex < dataIn.size(); tmpIndex++) {
2023  free (dataIn[tmpIndex]);
2024  }
2025  }
2026  else {
2027  UQ_FATAL_TEST_MACRO(true,
2028  m_env.worldRank(),
2029  "SequenceOfVectors<V,M>::unifiedReadContents()",
2030  "hdf file type not supported for multiple sub-environments yet");
2031  }
2032  }
2033 #endif
2034  else {
2035  UQ_FATAL_TEST_MACRO(true,
2036  m_env.worldRank(),
2037  "SequenceOfVectors<V,M>::unifiedReadContents()",
2038  "invalid file type");
2039  }
2040  m_env.closeFile(unifiedFilePtrSet,fileType);
2041  } // if (m_env.openUnifiedInputFile())
2042  } // if (m_env.inter0Rank() == (int) r)
2043  m_env.inter0Comm().Barrier();
2044  } // for r
2045  } // if (m_env.inter0Rank() >= 0)
2046  else {
2047  V tmpVec(m_vectorSpace.zeroVector());
2048  for (unsigned int i = 1; i < subReadSize; ++i) {
2049  this->setPositionValues(i,tmpVec);
2050  }
2051  }
2052 
2053  if (m_env.subDisplayFile()) {
2054  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedReadContents()"
2055  << ", fileName = " << fileName
2056  << std::endl;
2057  }
2058  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2059 
2060  return;
2061 }
2062 //---------------------------------------------------
2063 template <class V, class M>
2064 void
2065 SequenceOfVectors<V,M>::select(const std::vector<unsigned int>& idsOfUniquePositions)
2066 {
2067  UQ_FATAL_RC_MACRO(true,
2068  m_env.worldRank(),
2069  "SequenceOfVectors<V,M>::select()",
2070  "Code is not complete yet");
2071 
2072  if (&idsOfUniquePositions) {}; // just to remove compiler warning
2073 
2074  return;
2075 }
2076 //---------------------------------------------------
2077 template <class V, class M>
2078 void
2080  unsigned int initialPos,
2081  unsigned int spacing)
2082 {
2083  if (m_env.subDisplayFile()) {
2084  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::filter()"
2085  << ": initialPos = " << initialPos
2086  << ", spacing = " << spacing
2087  << ", subSequenceSize = " << this->subSequenceSize()
2088  << std::endl;
2089  }
2090 
2091  unsigned int i = 0;
2092  unsigned int j = initialPos;
2093  unsigned int originalSubSequenceSize = this->subSequenceSize();
2094  while (j < originalSubSequenceSize) {
2095  if (i != j) {
2096  //*m_env.subDisplayFile() << i << "--" << j << " ";
2097  delete m_seq[i];
2098  m_seq[i] = new V(*(m_seq[j]));
2099  }
2100  i++;
2101  j += spacing;
2102  }
2103 
2104  this->resetValues(i,originalSubSequenceSize-i);
2105  this->resizeSequence(i);
2106 
2107  if (m_env.subDisplayFile()) {
2108  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::filter()"
2109  << ": initialPos = " << initialPos
2110  << ", spacing = " << spacing
2111  << ", subSequenceSize = " << this->subSequenceSize()
2112  << std::endl;
2113  }
2114 
2115  return;
2116 }
2117 //---------------------------------------------------
2118 template <class V, class M>
2119 double
2121  unsigned int initialPos,
2122  unsigned int numPos) const
2123 {
2124  // This method requires *at least* two sequences. Error if there is only one.
2125  UQ_FATAL_RC_MACRO((m_env.numSubEnvironments() < 2),
2126  m_env.worldRank(),
2127  "SequenceOfVectors<V,M>::estimateConvBrooksGelman()",
2128  "At least two sequences required for Brooks-Gelman convergence test.");
2129 
2130  // TODO: Need special case for 1-dimensional parameter space.
2131 
2132  // Initialize with garbage to give the user a clue something is funky.
2133  double convMeasure = -1.0;
2134 
2135  // We only do the work on the subenvironment where the sequence data
2136  // is stashed.
2137  if( m_env.inter0Rank() >= 0 )
2138  {
2139  // Sanity Checks
2140 
2141  // REMEMBER: \psi is a *vector* of parameters
2142  // Get quantities we will use several times
2143  V psi_j_dot = m_vectorSpace.zeroVector();
2144  V psi_dot_dot = m_vectorSpace.zeroVector();
2145  V work = m_vectorSpace.zeroVector();
2146 
2147  // m = number of chains > 1
2148  // n = number of steps for which we are computing the metric
2149  int m = m_env.numSubEnvironments();
2150  int n = numPos;
2151 
2152  this->subMeanExtra ( initialPos, numPos, psi_j_dot );
2153  this->unifiedMeanExtra( initialPos, numPos, psi_dot_dot );
2154 
2155 #if 0
2156  std::cout << "psi_j_dot = " << psi_j_dot << std::endl;
2157  std::cout << "psi_dot_dot = " << psi_dot_dot << std::endl;
2158 #endif
2159 
2160  /* W = \frac{1}{m*(n-1)}*\sum_{j=1}^m \sum{t=1}^n
2161  (\psi_{jt} - \overline{\psi_{j\cdot}})*(\psi_{jt} - \overline{\psi_{j\cdot}})^T
2162  This corresponds to the "within-sequence" covariance matrix. */
2163  M* W_local = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2164  M* W = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2165  V psi_j_t = m_vectorSpace.zeroVector();
2166 
2167  // Sum within the chain
2168  for( unsigned int t = initialPos; t < initialPos+numPos; ++t )
2169  {
2170  psi_j_t = *(m_seq[t]);
2171 
2172  work = psi_j_t - psi_j_dot;
2173 
2174  (*W_local) += matrixProduct( work, work );
2175  }
2176 
2177  // Now do the sum over the chains
2178  // W will be available on all inter0 processors
2179  W_local->mpiSum( m_env.inter0Comm(), (*W) );
2180 
2181  (*W) = 1.0/(double(m)*(double(n)-1.0)) * (*W);
2182 
2183 #if 0
2184  std::cout << "n, m = " << n << ", " << m << std::endl;
2185  std::cout << "W_local = " << *W_local << std::endl;
2186  std::cout << "W = " << *W << std::endl;
2187 #endif
2188 
2189  // Need to delete pointers to temporary covariance matrices
2190  delete W_local;
2191 
2192  /* B/n = \frac{1}{m-1}\sum_{j=1}^m
2193  (\overline{\psi_{j\cdot}} - \overline{\psi_{\cdot \cdot}})*
2194  (\overline{\psi_{j\cdot}} - \overline{\psi_{\cdot \cdot}})^T
2195  This corresponds to the "between-sequence" covariance matrix. */
2196  M* B_over_n_local = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2197  M* B_over_n = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2198 
2199  work = psi_j_dot - psi_dot_dot;
2200  (*B_over_n_local) = matrixProduct( work, work );
2201 
2202  B_over_n_local->mpiSum( m_env.inter0Comm(), (*B_over_n) );
2203 
2204  // Need to delete pointers to temporary covariance matrices
2205  delete B_over_n_local;
2206 
2207  (*B_over_n) = 1.0/(double(m)-1.0) * (*B_over_n);
2208 
2209 #if 0
2210  std::cout << "B_over_n = " << *B_over_n << std::endl;
2211 #endif
2212 
2213 
2214  /* R_p = (n-1)/n + (m+1)/m * \lambda
2215  \lambda = largest eigenvalue of W^{-1}*B/n */
2216  M* A = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2217 
2218  W->invertMultiply( *B_over_n, *A );
2219 
2220 #if 0
2221  std::cout << "A = " << *A << std::endl;
2222  std::cout.flush();
2223 #endif
2224  // Need to delete pointers to temporary covariance matrices
2225  delete W;
2226  delete B_over_n;
2227 
2228  double eigenValue;
2229  V eigenVector = m_vectorSpace.zeroVector();
2230 
2231  A->largestEigen( eigenValue, eigenVector );
2232 
2233  // Need to delete pointers to temporary covariance matrices
2234  delete A;
2235 
2236  // Now, finally compute the final convMeasure
2237  convMeasure = (double(n)-1.0)/double(n) + (double(m)+1.0)/double(m)*eigenValue;
2238 
2239  } // End of check on inter0Rank
2240 
2241  //TODO: Do we need a Barrier here?
2242  //TODO: Error checking on MPI.
2243  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2244 
2245  return convMeasure;
2246 }
2247 //---------------------------------------------------
2248 template <class V, class M>
2249 void
2251  unsigned int initialPos,
2252  unsigned int spacing,
2253  unsigned int numPos,
2254  unsigned int paramId,
2255  ScalarSequence<double>& scalarSeq) const
2256 {
2257  scalarSeq.resizeSequence(numPos);
2258  if (spacing == 1) {
2259  for (unsigned int j = 0; j < numPos; ++j) {
2260  scalarSeq[j] = (*(m_seq[initialPos+j ]))[paramId];
2261  }
2262  }
2263  else {
2264  for (unsigned int j = 0; j < numPos; ++j) {
2265  scalarSeq[j] = (*(m_seq[initialPos+j*spacing]))[paramId];
2266  }
2267  }
2268 
2269  return;
2270 }
2271 // Private methods ------------------------------------
2272 template <class V, class M>
2273 void
2275 {
2277  for (unsigned int i = 0; i < (unsigned int) m_seq.size(); ++i) {
2278  if (m_seq[i]) {
2279  delete m_seq[i];
2280  m_seq[i] = NULL;
2281  }
2282  }
2283  m_seq.resize(src.subSequenceSize(),NULL);
2284  for (unsigned int i = 0; i < m_seq.size(); ++i) {
2285  m_seq[i] = new V(*(src.m_seq[i]));
2286  }
2287 
2288  return;
2289 }
2290 //---------------------------------------------------
2291 template <class V, class M>
2292 void
2294  unsigned int initialPos,
2295  unsigned int spacing,
2296  unsigned int numPos,
2297  unsigned int paramId,
2298  std::vector<double>& rawData) const
2299 {
2300  rawData.resize(numPos);
2301  if (spacing == 1) {
2302  for (unsigned int j = 0; j < numPos; ++j) {
2303  rawData[j] = (*(m_seq[initialPos+j ]))[paramId];
2304  }
2305  }
2306  else {
2307  for (unsigned int j = 0; j < numPos; ++j) {
2308  rawData[j] = (*(m_seq[initialPos+j*spacing]))[paramId];
2309  }
2310  }
2311 
2312  return;
2313 }
2314 
2315 // --------------------------------------------------
2316 // Methods conditionally available ------------------
2317 // --------------------------------------------------
2318 // --------------------------------------------------
2319 #ifdef UQ_SEQ_VEC_USES_OPERATOR
2320 template <class V, class M>
2321 const V*
2322 SequenceOfVectors<V,M>::operator[](unsigned int posId) const
2323 {
2324  UQ_FATAL_TEST_MACRO((posId >= this->subSequenceSize()),
2325  m_env.worldRank(),
2326  "SequenceOfVectorss<V,M>::operator[] const",
2327  "posId > subSequenceSize()");
2328 
2329  return (const V*) (m_seq[posId]);
2330 }
2331 // --------------------------------------------------
2332 template <class V, class M>
2333 const V*&
2334 SequenceOfVectors<V,M>::operator[](unsigned int posId)
2335 {
2336  UQ_FATAL_TEST_MACRO((posId >= this->subSequenceSize()),
2337  m_env.worldRank(),
2338  "SequenceOfVectorss<V,M>::operator[] const",
2339  "posId > subSequenceSize()");
2340 
2341  return m_seq[posId];
2342 }
2343 #endif
2344 
2345 // --------------------------------------------------
2346 // --------------------------------------------------
2347 // --------------------------------------------------
2348 
2349 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
2350 template <class V, class M>
2351 void
2352 SequenceOfVectors<V,M>::subUniformlySampledMdf(
2353  const V& numEvaluationPointsVec,
2354  ArrayOfOneDGrids <V,M>& mdfGrids,
2355  ArrayOfOneDTables<V,M>& mdfValues) const
2356 {
2357  V minDomainValues(m_vectorSpace.zeroVector());
2358  V maxDomainValues(m_vectorSpace.zeroVector());
2359 
2360  ScalarSequence<double> data(m_env,0,"");
2361 
2362  unsigned int numParams = this->vectorSizeLocal();
2363  for (unsigned int i = 0; i < numParams; ++i) {
2364  this->extractScalarSeq(0, // initialPos
2365  1, // spacing
2366  subSequenceSize(), // numPos
2367  i,
2368  data);
2369 
2370  std::vector<double> aMdf(0);
2371  data.subUniformlySampledMdf((unsigned int) numEvaluationPointsVec[i],
2372  minDomainValues[i],
2373  maxDomainValues[i],
2374  aMdf);
2375  mdfValues.setOneDTable(i,aMdf);
2376  }
2377 
2378  mdfGrids.setUniformGrids(numEvaluationPointsVec,
2379  minDomainValues,
2380  maxDomainValues);
2381 
2382  return;
2383 }
2384 #endif
2385 
2386 // --------------------------------------------------
2387 // --------------------------------------------------
2388 // --------------------------------------------------
2389 
2390 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2391 template <class V, class M>
2392 void
2393 SequenceOfVectors<V,M>::subMeanCltStd(
2394  unsigned int initialPos,
2395  unsigned int numPos,
2396  const V& meanVec,
2397  V& stdVec) const
2398 {
2399  bool bRC = ((initialPos < this->subSequenceSize()) &&
2400  (0 < numPos ) &&
2401  ((initialPos+numPos) <= this->subSequenceSize()) &&
2402  (this->vectorSizeLocal() == meanVec.sizeLocal() ) &&
2403  (this->vectorSizeLocal() == stdVec.sizeLocal() ));
2404  UQ_FATAL_TEST_MACRO(bRC == false,
2405  m_env.worldRank(),
2406  "SequenceOfVectors<V,M>::subMeanCltStd()",
2407  "invalid input data");
2408 
2409  ScalarSequence<double> data(m_env,0,"");
2410 
2411  unsigned int numParams = this->vectorSizeLocal();
2412  for (unsigned int i = 0; i < numParams; ++i) {
2413  this->extractScalarSeq(initialPos,
2414  1, // spacing
2415  numPos,
2416  i,
2417  data);
2418  stdVec[i] = data.subMeanCltStd(0,
2419  numPos,
2420  meanVec[i]);
2421  }
2422 
2423  return;
2424 }
2425 
2426 template <class V, class M>
2427 void
2428 SequenceOfVectors<V,M>::unifiedMeanCltStd(
2429  unsigned int initialPos,
2430  unsigned int numPos,
2431  const V& unifiedMeanVec,
2432  V& unifiedSamVec) const
2433 {
2434  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2435  (0 < numPos ) &&
2436  ((initialPos+numPos) <= this->subSequenceSize() ) &&
2437  (this->vectorSizeLocal() == unifiedMeanVec.sizeLocal()) &&
2438  (this->vectorSizeLocal() == unifiedSamVec.sizeLocal() ));
2439  UQ_FATAL_TEST_MACRO(bRC == false,
2440  m_env.worldRank(),
2441  "SequenceOfVectors<V,M>::unifiedMeanCltStd()",
2442  "invalid input data");
2443 
2444  ScalarSequence<double> data(m_env,0,"");
2445 
2446  unsigned int numParams = this->vectorSizeLocal();
2447  for (unsigned int i = 0; i < numParams; ++i) {
2448  this->extractScalarSeq(initialPos,
2449  1, // spacing
2450  numPos,
2451  i,
2452  data);
2453  unifiedSamVec[i] = data.unifiedMeanCltStd(m_vectorSpace.numOfProcsForStorage() == 1,
2454  0,
2455  numPos,
2456  unifiedMeanVec[i]);
2457  }
2458 
2459  return;
2460 }
2461 //---------------------------------------------------
2462 template <class V, class M>
2463 void
2464 SequenceOfVectors<V,M>::bmm(
2465  unsigned int initialPos,
2466  unsigned int batchLength,
2467  V& bmmVec) const
2468 {
2469  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2470  (batchLength < (this->subSequenceSize()-initialPos)) &&
2471  (this->vectorSizeLocal() == bmmVec.sizeLocal() ));
2472  UQ_FATAL_TEST_MACRO(bRC == false,
2473  m_env.worldRank(),
2474  "SequenceOfVectors<V,M>::bmm()",
2475  "invalid input data");
2476 
2477  ScalarSequence<double> data(m_env,0,"");
2478 
2479  unsigned int numParams = this->vectorSizeLocal();
2480  for (unsigned int i = 0; i < numParams; ++i) {
2481  this->extractScalarSeq(initialPos,
2482  1, // spacing
2483  this->subSequenceSize()-initialPos,
2484  i,
2485  data);
2486  bmmVec[i] = data.bmm(0,
2487  batchLength);
2488  }
2489 
2490  return;
2491 }
2492 //---------------------------------------------------
2493 template <class V, class M>
2494 void
2495 SequenceOfVectors<V,M>::fftForward(
2496  unsigned int initialPos,
2497  unsigned int fftSize,
2498  unsigned int paramId,
2499  std::vector<std::complex<double> >& fftResult) const
2500 {
2501  bool bRC = ((initialPos < this->subSequenceSize()) &&
2502  (paramId < this->vectorSizeLocal()) &&
2503  (0 < fftSize ) &&
2504  ((initialPos+fftSize) <= this->subSequenceSize()) &&
2505  (fftSize < this->subSequenceSize()));
2506  UQ_FATAL_TEST_MACRO(bRC == false,
2507  m_env.worldRank(),
2508  "SequenceOfVectors<V,M>::fftForward()",
2509  "invalid input data");
2510 
2511  std::vector<double> rawData(fftSize,0.);
2512  this->extractRawData(initialPos,
2513  1, // spacing
2514  fftSize,
2515  paramId,
2516  rawData);
2517 
2518  m_fftObj->forward(rawData,fftSize,fftResult);
2519 
2520  return;
2521 }
2522 //---------------------------------------------------
2523 template <class V, class M>
2524 void
2525 SequenceOfVectors<V,M>::psd(
2526  unsigned int initialPos,
2527  unsigned int numBlocks,
2528  double hopSizeRatio,
2529  unsigned int paramId,
2530  std::vector<double>& psdResult) const
2531 {
2532  bool bRC = ((initialPos < this->subSequenceSize()) &&
2533  (paramId < this->vectorSizeLocal()));
2534  UQ_FATAL_TEST_MACRO(bRC == false,
2535  m_env.worldRank(),
2536  "SequenceOfVectors<V,M>::psd()",
2537  "invalid input data");
2538 
2539  ScalarSequence<double> data(m_env,0,"");
2540 
2541  this->extractScalarSeq(initialPos,
2542  1, // spacing
2543  this->subSequenceSize()-initialPos,
2544  paramId,
2545  data);
2546  data.psd(0,
2547  numBlocks,
2548  hopSizeRatio,
2549  psdResult);
2550 
2551  return;
2552 }
2553 //---------------------------------------------------
2554 template <class V, class M>
2555 void
2556 SequenceOfVectors<V,M>::psdAtZero(
2557  unsigned int initialPos,
2558  unsigned int numBlocks,
2559  double hopSizeRatio,
2560  V& psdVec) const
2561 {
2562  bool bRC = ((initialPos < this->subSequenceSize()) &&
2563  (this->vectorSizeLocal() == psdVec.sizeLocal()));
2564  UQ_FATAL_TEST_MACRO(bRC == false,
2565  m_env.worldRank(),
2566  "SequenceOfVectors<V,M>::psdAtZero()",
2567  "invalid input data");
2568 
2569  ScalarSequence<double> data(m_env,0,"");
2570  std::vector<double> psdResult(0,0.); // size will be determined by 'data.psd()'
2571 
2572  unsigned int numParams = this->vectorSizeLocal();
2573  for (unsigned int i = 0; i < numParams; ++i) {
2574  this->extractScalarSeq(initialPos,
2575  1, // spacing
2576  this->subSequenceSize()-initialPos,
2577  i,
2578  data);
2579  data.psd(0,
2580  numBlocks,
2581  hopSizeRatio,
2582  psdResult);
2583  psdVec[i] = psdResult[0];
2584  //*m_env.subDisplayFile() << "psdResult[0] = " << psdResult[0] << std::endl;
2585  }
2586 
2587  return;
2588 }
2589 //---------------------------------------------------
2590 template <class V, class M>
2591 void
2592 SequenceOfVectors<V,M>::geweke(
2593  unsigned int initialPos,
2594  double ratioNa,
2595  double ratioNb,
2596  V& gewVec) const
2597 {
2598  bool bRC = ((initialPos < this->subSequenceSize()) &&
2599  (this->vectorSizeLocal() == gewVec.sizeLocal() ));
2600  UQ_FATAL_TEST_MACRO(bRC == false,
2601  m_env.worldRank(),
2602  "SequenceOfVectors<V,M>::geweke()",
2603  "invalid input data");
2604 
2605  unsigned int numPos = this->subSequenceSize() - initialPos;
2606  ScalarSequence<double> data(m_env,0,"");
2607 
2608  unsigned int numParams = this->vectorSizeLocal();
2609  for (unsigned int i = 0; i < numParams; ++i) {
2610  this->extractScalarSeq(initialPos,
2611  1, // spacing
2612  numPos,
2613  i,
2614  data);
2615  gewVec[i] = data.geweke(0,
2616  ratioNa,
2617  ratioNb);
2618  }
2619 
2620  return;
2621 }
2622 //---------------------------------------------------
2623 template <class V, class M>
2624 void
2625 SequenceOfVectors<V,M>::meanStacc(
2626  unsigned int initialPos,
2627  V& meanStaccVec) const
2628 {
2629  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2630  (this->vectorSizeLocal() == meanStaccVec.sizeLocal()));
2631  UQ_FATAL_TEST_MACRO(bRC == false,
2632  m_env.worldRank(),
2633  "SequenceOfVectors<V,M>::meanStacc()",
2634  "invalid input data");
2635 
2636  unsigned int numPos = this->subSequenceSize() - initialPos;
2637  ScalarSequence<double> data(m_env,0,"");
2638 
2639  unsigned int numParams = this->vectorSizeLocal();
2640  for (unsigned int i = 0; i < numParams; ++i) {
2641  this->extractScalarSeq(initialPos,
2642  1, // spacing
2643  numPos,
2644  i,
2645  data);
2646  meanStaccVec[i] = data.meanStacc(0);
2647  }
2648 
2649  return;
2650 }
2651 //---------------------------------------------------
2652 template <class V, class M>
2653 void
2654 SequenceOfVectors<V,M>::subCdfPercentageRange(
2655  unsigned int initialPos,
2656  unsigned int numPos,
2657  double range, // \in [0,1]
2658  V& lowerVec,
2659  V& upperVec) const
2660 {
2661  bool bRC = ((0 < numPos ) &&
2662  ((initialPos+numPos) <= this->subSequenceSize()) &&
2663  (this->vectorSizeLocal() == lowerVec.sizeLocal() ) &&
2664  (this->vectorSizeLocal() == upperVec.sizeLocal() ));
2665  UQ_FATAL_TEST_MACRO(bRC == false,
2666  m_env.worldRank(),
2667  "SequenceOfVectors<V,M>::subCdfPercentageRange()",
2668  "invalid input data");
2669 
2670  unsigned int numParams = this->vectorSizeLocal();
2671  ScalarSequence<double> data(m_env,0,"");
2672 
2673  for (unsigned int i = 0; i < numParams; ++i) {
2674  this->extractScalarSeq(initialPos,
2675  1, // spacing
2676  numPos,
2677  i,
2678  data);
2679  data.subCdfPercentageRange(0,
2680  numPos,
2681  range,
2682  lowerVec[i],
2683  upperVec[i]);
2684  }
2685 
2686  return;
2687 }
2688 //---------------------------------------------------
2689 template <class V, class M>
2690 void
2691 SequenceOfVectors<V,M>::unifiedCdfPercentageRange(
2692  unsigned int initialPos,
2693  unsigned int numPos,
2694  double range, // \in [0,1]
2695  V& lowerVec,
2696  V& upperVec) const
2697 {
2698  bool bRC = ((0 < numPos ) &&
2699  ((initialPos+numPos) <= this->subSequenceSize()) &&
2700  (this->vectorSizeLocal() == lowerVec.sizeLocal() ) &&
2701  (this->vectorSizeLocal() == upperVec.sizeLocal() ));
2702  UQ_FATAL_TEST_MACRO(bRC == false,
2703  m_env.worldRank(),
2704  "SequenceOfVectors<V,M>::unifiedCdfPercentageRange()",
2705  "invalid input data");
2706 
2707  unsigned int numParams = this->vectorSizeLocal();
2708  ScalarSequence<double> data(m_env,0,"");
2709 
2710  for (unsigned int i = 0; i < numParams; ++i) {
2711  this->extractScalarSeq(initialPos,
2712  1, // spacing
2713  numPos,
2714  i,
2715  data);
2716  data.unifiedCdfPercentageRange(m_vectorSpace.numOfProcsForStorage() == 1,
2717  0,
2718  numPos,
2719  range,
2720  lowerVec[i],
2721  upperVec[i]);
2722  }
2723 
2724  return;
2725 }
2726 //---------------------------------------------------
2727 template <class V, class M>
2728 void
2729 SequenceOfVectors<V,M>::subCdfStacc(
2730  unsigned int initialPos,
2731  std::vector<V*>& cdfStaccVecs,
2732  std::vector<V*>& cdfStaccVecsUp,
2733  std::vector<V*>& cdfStaccVecsLow,
2734  std::vector<V*>& sortedDataVecs) const
2735 {
2736  bool bRC = (initialPos < this->subSequenceSize());
2737  UQ_FATAL_TEST_MACRO(bRC == false,
2738  m_env.worldRank(),
2739  "SequenceOfVectors<V,M>::subCdfStacc()",
2740  "invalid input data");
2741 
2742  unsigned int numPos = this->subSequenceSize() - initialPos;
2743  unsigned int numEvals = numPos;
2744  for (unsigned int j = 0; j < numEvals; ++j) {
2745  cdfStaccVecs [j] = new V(m_vectorSpace.zeroVector());
2746  cdfStaccVecsUp [j] = new V(m_vectorSpace.zeroVector());
2747  cdfStaccVecsLow[j] = new V(m_vectorSpace.zeroVector());
2748  sortedDataVecs [j] = new V(m_vectorSpace.zeroVector());
2749  }
2750  std::vector<double> cdfStaccs (numEvals,0.);
2751  std::vector<double> cdfStaccsup (numEvals,0.);
2752  std::vector<double> cdfStaccslow(numEvals,0.);
2753 
2754  ScalarSequence<double> data (m_env,0,"");
2755  ScalarSequence<double> sortedData(m_env,0,"");
2756  unsigned int numParams = this->vectorSizeLocal();
2757  for (unsigned int i = 0; i < numParams; ++i) {
2758  this->extractScalarSeq(initialPos,
2759  1, // spacing
2760  numPos,
2761  i,
2762  data);
2763  //std::cout << "x-data" << data<< std::endl;
2764  data.subSort(initialPos,sortedData);
2765  data.subCdfStacc(initialPos,
2766  cdfStaccs,
2767  cdfStaccsup,
2768  cdfStaccslow,
2769  sortedData);
2770 
2771  for (unsigned int j = 0; j < numEvals; ++j) {
2772  (*sortedDataVecs [j])[i] = sortedData [j];
2773  (*cdfStaccVecs [j])[i] = cdfStaccs [j];
2774  (*cdfStaccVecsUp [j])[i] = cdfStaccsup [j];
2775  (*cdfStaccVecsLow[j])[i] = cdfStaccslow[j];
2776  }
2777  }
2778 
2779  return;
2780 }
2781 //---------------------------------------------------
2782 template <class V, class M>
2783 void
2784 SequenceOfVectors<V,M>::subCdfStacc(
2785  unsigned int initialPos,
2786  const std::vector<V*>& evalPositionsVecs,
2787  std::vector<V*>& cdfStaccVecs) const
2788 {
2789  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2790  (0 < evalPositionsVecs.size()) &&
2791  (evalPositionsVecs.size() == cdfStaccVecs.size() ));
2792  UQ_FATAL_TEST_MACRO(bRC == false,
2793  m_env.worldRank(),
2794  "SequenceOfVectors<V,M>::subCdfStacc()",
2795  "invalid input data");
2796 
2797  unsigned int numPos = this->subSequenceSize() - initialPos;
2798  ScalarSequence<double> data(m_env,0,"");
2799 
2800  unsigned int numEvals = evalPositionsVecs.size();
2801  for (unsigned int j = 0; j < numEvals; ++j) {
2802  cdfStaccVecs[j] = new V(m_vectorSpace.zeroVector());
2803  }
2804  std::vector<double> evalPositions(numEvals,0.);
2805  std::vector<double> cdfStaccs (numEvals,0.);
2806 
2807  unsigned int numParams = this->vectorSizeLocal();
2808  for (unsigned int i = 0; i < numParams; ++i) {
2809  this->extractScalarSeq(initialPos,
2810  1, // spacing
2811  numPos,
2812  i,
2813  data);
2814 
2815  for (unsigned int j = 0; j < numEvals; ++j) {
2816  evalPositions[j] = (*evalPositionsVecs[j])[i];
2817  }
2818 
2819  data.subCdfStacc(0,
2820  evalPositions,
2821  cdfStaccs);
2822 
2823  for (unsigned int j = 0; j < numEvals; ++j) {
2824  (*cdfStaccVecs[j])[i] = cdfStaccs[j];
2825  }
2826  }
2827 
2828  return;
2829 }
2830 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2831 
2832 // --------------------------------------------------
2833 // --------------------------------------------------
2834 // --------------------------------------------------
2835 
2836 #ifdef UQ_CODE_HAS_MONITORS
2837 template <class V, class M>
2838 void
2839 SequenceOfVectors<V,M>::subMeanMonitorAlloc(unsigned int numberOfMonitorPositions)
2840 {
2841  m_subMeanMonitorPosSeq = new ScalarSequence<double>(m_env, numberOfMonitorPositions,(m_name+"_subMeanMonitorPosSeq").c_str());
2842  m_subMeanVecSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanVecSeq").c_str() );
2843  m_subMeanCltStdSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanCltStdSeq").c_str() );
2844 
2845  return;
2846 }
2847 // --------------------------------------------------
2848 template <class V, class M>
2849 void
2850 SequenceOfVectors<V,M>::subMeanInter0MonitorAlloc(unsigned int numberOfMonitorPositions)
2851 {
2852  m_subMeanInter0MonitorPosSeq = new ScalarSequence<double>(m_env, numberOfMonitorPositions,(m_name+"_subMeanInter0MonitorPosSeq").c_str() );
2853  m_subMeanInter0Mean = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0MeanSeq").c_str() );
2854  m_subMeanInter0Clt95 = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0Clt95Seq").c_str() );
2855  m_subMeanInter0Empirical90 = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0Empirical90Seq").c_str());
2856  m_subMeanInter0Min = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0MinSeq").c_str() );
2857  m_subMeanInter0Max = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0MaxSeq").c_str() );
2858 
2859  return;
2860 }
2861 // --------------------------------------------------
2862 template <class V, class M>
2863 void
2864 SequenceOfVectors<V,M>::unifiedMeanMonitorAlloc(unsigned int numberOfMonitorPositions)
2865 {
2866  m_unifiedMeanMonitorPosSeq = new ScalarSequence<double>(m_env, numberOfMonitorPositions,(m_name+"_unifiedMeanMonitorPosSeq").c_str());
2867  m_unifiedMeanVecSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_unifiedMeanVecSeq").c_str() );
2868  m_unifiedMeanCltStdSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_unifiedMeanCltStdSeq").c_str() );
2869 
2870  return;
2871 }
2872 // --------------------------------------------------
2873 template <class V, class M>
2874 void
2875 SequenceOfVectors<V,M>::subMeanMonitorRun(unsigned int monitorPosition,
2876  V& subMeanVec,
2877  V& subMeanCltStd)
2878 {
2879  this->subMeanExtra(0,
2880  monitorPosition,
2881  subMeanVec);
2882 
2883  this->subMeanCltStd(0,
2884  monitorPosition,
2885  subMeanVec,
2886  subMeanCltStd);
2887 
2888  return;
2889 }
2890 // --------------------------------------------------
2891 template <class V, class M>
2892 void
2893 SequenceOfVectors<V,M>::subMeanInter0MonitorRun(unsigned int monitorPosition,
2894  V& subMeanInter0Mean,
2895  V& subMeanInter0Clt95,
2896  V& subMeanInter0Empirical90,
2897  V& subMeanInter0Min,
2898  V& subMeanInter0Max)
2899 {
2900  V subMeanVec(m_vectorSpace.zeroVector());
2901  this->subMeanExtra(0,
2902  monitorPosition,
2903  subMeanVec);
2904 
2905  subMeanVec.mpiAllReduce(RawValue_MPI_SUM,m_env.inter0Comm(),subMeanInter0Mean);
2906  subMeanInter0Mean /= ((double) m_env.inter0Comm().NumProc());
2907 
2908  V subMeanInter0CltVariance = subMeanVec-subMeanInter0Mean;
2909  subMeanInter0CltVariance *= subMeanInter0CltVariance;
2910  subMeanInter0CltVariance.mpiAllReduce(RawValue_MPI_SUM,m_env.inter0Comm(),subMeanInter0Clt95);
2911  subMeanInter0Clt95 /= ((double) (m_env.inter0Comm().NumProc()-1));
2912  subMeanInter0Clt95 /= ((double) (m_env.inter0Comm().NumProc()-1));
2913  subMeanInter0Clt95.cwSqrt();
2914  subMeanInter0Clt95 *= 3.;
2915 
2916  V subMeanInter0Quantile5(m_vectorSpace.zeroVector());
2917  subMeanVec.mpiAllQuantile(.05,m_env.inter0Comm(),subMeanInter0Quantile5);
2918  V subMeanInter0Quantile95(m_vectorSpace.zeroVector());
2919  subMeanVec.mpiAllQuantile(.95,m_env.inter0Comm(),subMeanInter0Quantile95);
2920  subMeanInter0Empirical90 = subMeanInter0Quantile95 - subMeanInter0Quantile5;
2921 
2922  subMeanVec.mpiAllReduce(RawValue_MPI_MIN,m_env.inter0Comm(),subMeanInter0Min);
2923 
2924  subMeanVec.mpiAllReduce(RawValue_MPI_MAX,m_env.inter0Comm(),subMeanInter0Max);
2925 
2926  return;
2927 }
2928 // --------------------------------------------------
2929 template <class V, class M>
2930 void
2931 SequenceOfVectors<V,M>::unifiedMeanMonitorRun(unsigned int monitorPosition,
2932  V& unifiedMeanVec,
2933  V& unifiedMeanCltStd)
2934 {
2935  this->unifiedMeanExtra(0,
2936  monitorPosition,
2937  unifiedMeanVec);
2938 
2939  this->unifiedMeanCltStd(0,
2940  monitorPosition,
2941  unifiedMeanVec,
2942  unifiedMeanCltStd);
2943  return;
2944 }
2945 // --------------------------------------------------
2946 template <class V, class M>
2947 void
2948 SequenceOfVectors<V,M>::subMeanMonitorStore(unsigned int i,
2949  unsigned int monitorPosition,
2950  const V& subMeanVec,
2951  const V& subMeanCltStd)
2952 {
2953  (*m_subMeanMonitorPosSeq)[i] = monitorPosition;
2954  m_subMeanVecSeq->setPositionValues(i,subMeanVec);
2955  m_subMeanCltStdSeq->setPositionValues(i,subMeanCltStd);
2956 
2957  return;
2958 }
2959 // --------------------------------------------------
2960 template <class V, class M>
2961 void
2962 SequenceOfVectors<V,M>::subMeanInter0MonitorStore(unsigned int i,
2963  unsigned int monitorPosition,
2964  const V& subMeanInter0Mean,
2965  const V& subMeanInter0Clt95,
2966  const V& subMeanInter0Empirical90,
2967  const V& subMeanInter0Min,
2968  const V& subMeanInter0Max)
2969 {
2970  (*m_subMeanInter0MonitorPosSeq)[i] = monitorPosition;
2971  m_subMeanInter0Mean->setPositionValues(i,subMeanInter0Mean);
2972  m_subMeanInter0Clt95->setPositionValues(i,subMeanInter0Clt95);
2973  m_subMeanInter0Empirical90->setPositionValues(i,subMeanInter0Empirical90);
2974  m_subMeanInter0Min->setPositionValues(i,subMeanInter0Min);
2975  m_subMeanInter0Max->setPositionValues(i,subMeanInter0Max);
2976 
2977  return;
2978 }
2979 // --------------------------------------------------
2980 template <class V, class M>
2981 void
2982 SequenceOfVectors<V,M>::unifiedMeanMonitorStore(unsigned int i,
2983  unsigned int monitorPosition,
2984  V& unifiedMeanVec,
2985  V& unifiedMeanCltStd)
2986 {
2987  (*m_unifiedMeanMonitorPosSeq)[i] = monitorPosition;
2988  m_unifiedMeanVecSeq->setPositionValues(i,unifiedMeanVec);
2989  m_unifiedMeanCltStdSeq->setPositionValues(i,unifiedMeanCltStd);
2990 
2991  return;
2992 }
2993 // --------------------------------------------------
2994 template <class V, class M>
2995 void
2996 SequenceOfVectors<V,M>::subMeanMonitorWrite(std::ofstream& ofs)
2997 {
2998  m_subMeanMonitorPosSeq->subWriteContents(0,m_subMeanMonitorPosSeq->subSequenceSize(),ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2999  m_subMeanVecSeq->subWriteContents (0,m_subMeanVecSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3000  m_subMeanCltStdSeq->subWriteContents (0,m_subMeanCltStdSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3001 
3002  return;
3003 }
3004 // --------------------------------------------------
3005 template <class V, class M>
3006 void
3007 SequenceOfVectors<V,M>::subMeanInter0MonitorWrite(std::ofstream& ofs)
3008 {
3009  m_subMeanInter0MonitorPosSeq->subWriteContents(0,m_subMeanInter0MonitorPosSeq->subSequenceSize(),ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3010  m_subMeanInter0Mean->subWriteContents (0,m_subMeanInter0Mean->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3011  m_subMeanInter0Clt95->subWriteContents (0,m_subMeanInter0Clt95->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3012  m_subMeanInter0Empirical90->subWriteContents (0,m_subMeanInter0Empirical90->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3013  m_subMeanInter0Min->subWriteContents (0,m_subMeanInter0Min->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3014  m_subMeanInter0Max->subWriteContents (0,m_subMeanInter0Max->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
3015 
3016  return;
3017 }
3018 // --------------------------------------------------
3019 template <class V, class M>
3020 void
3021 SequenceOfVectors<V,M>::unifiedMeanMonitorWrite(std::ofstream& ofs)
3022 {
3023  // std::set<unsigned int> tmpSet;
3024  // tmpSet.insert(0);
3025  m_unifiedMeanMonitorPosSeq->subWriteContents(0,m_unifiedMeanMonitorPosSeq->subSequenceSize(),ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m" // Yes, 'subWriteContents()'
3026  m_unifiedMeanVecSeq->subWriteContents (0,m_unifiedMeanVecSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m" // Yes, 'subWriteContents()'
3027  m_unifiedMeanCltStdSeq->subWriteContents (0,m_unifiedMeanCltStdSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m" // Yes, 'subWriteContents()'
3028 
3029  return;
3030 }
3031 // --------------------------------------------------
3032 template <class V, class M>
3033 void
3034 SequenceOfVectors<V,M>::subMeanMonitorFree()
3035 {
3036  delete m_subMeanMonitorPosSeq;
3037  m_subMeanMonitorPosSeq = NULL;
3038  delete m_subMeanVecSeq;
3039  m_subMeanVecSeq = NULL;
3040  delete m_subMeanCltStdSeq;
3041  m_subMeanCltStdSeq = NULL;
3042 
3043  return;
3044 }
3045 // --------------------------------------------------
3046 template <class V, class M>
3047 void
3048 SequenceOfVectors<V,M>::subMeanInter0MonitorFree()
3049 {
3050  delete m_subMeanInter0MonitorPosSeq;
3051  m_subMeanInter0MonitorPosSeq = NULL;
3052  delete m_subMeanInter0Mean;
3053  m_subMeanInter0Mean = NULL;
3054  delete m_subMeanInter0Clt95;
3055  m_subMeanInter0Clt95 = NULL;
3056  delete m_subMeanInter0Empirical90;
3057  m_subMeanInter0Empirical90 = NULL;
3058  delete m_subMeanInter0Min;
3059  m_subMeanInter0Min = NULL;
3060  delete m_subMeanInter0Max;
3061  m_subMeanInter0Max = NULL;
3062 
3063  return;
3064 }
3065 // --------------------------------------------------
3066 template <class V, class M>
3067 void
3068 SequenceOfVectors<V,M>::unifiedMeanMonitorFree()
3069 {
3070  delete m_unifiedMeanMonitorPosSeq;
3071  m_unifiedMeanMonitorPosSeq = NULL;
3072  delete m_unifiedMeanVecSeq;
3073  m_unifiedMeanVecSeq = NULL;
3074  delete m_unifiedMeanCltStdSeq;
3075  m_unifiedMeanCltStdSeq = NULL;
3076 
3077  return;
3078 }
3079 #endif // #ifdef UQ_CODE_HAS_MONITORS
3080 
3081 } // End namespace QUESO
3082 
T subInterQuantileRange(unsigned int initialPos) const
Returns the interquartile range of the values in the sub-sequence.
double MiscGetEllapsedSeconds(struct timeval *timeval0)
void subUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &cdfGrids, ArrayOfOneDTables< V, M > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void deleteStoredVectors()
Deletes all the stored vectors.
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.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
Class to accommodate arrays of one-dimensional grid.
void unifiedMedianExtra(unsigned int initialPos, unsigned int numPos, V &unifiedMedianVec) const
Finds the median value of the unfed sequence, considering numPos positions starting at position initi...
void subHistogram(unsigned int initialPos, const V &minVec, const V &maxVec, std::vector< V * > &centersForAllBins, std::vector< V * > &quanttsForAllBins) const
Calculates the histogram of the sub-sequence.
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...
Struct for handling data input and output from files.
Definition: Environment.h:66
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
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.
void subMedianExtra(unsigned int initialPos, unsigned int numPos, V &medianVec) const
Finds the median value of the sub-sequence, considering numPos positions starting at position initial...
void unifiedInterQuantileRange(unsigned int initialPos, V &unifiedIqrVec) const
Returns the interquartile range of the values in the unified sequence.
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2411
void unifiedPopulationVariance(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedPopVec) const
Finds the population variance of the unified sequence, considering numPos positions starting at posit...
void subMeanExtra(unsigned int initialPos, unsigned int numPos, V &meanVec) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, ScalarSequence< double > &scalarSeq) const
Extracts a sequence of scalars.
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
Definition: Defines.h:207
void subSampleStd(unsigned int initialPos, unsigned int numPos, const V &meanVec, V &stdVec) const
Finds the sample standard deviation of the sub-sequence, considering numPos positions starting at pos...
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
SequenceOfVectors< V, M > & operator=(const SequenceOfVectors< V, M > &rhs)
Copies values from rhs to this.
const int UQ_OK_RC
Definition: Defines.h:76
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...
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:75
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:89
double estimateConvBrooksGelman(unsigned int initialPos, unsigned int numPos) const
Estimates convergence rate using Brooks &amp; Gelman method.
std::vector< const V * >::iterator seqVectorPositionIteratorTypedef
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...
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:223
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
void unifiedScalesForKde(unsigned int initialPos, const V &unifiedIqrVec, unsigned int kdeDimension, V &unifiedScaleVec) const
Selects the scales (bandwidth) for the kernel density estimation, considering the unified sequence...
void autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
void unifiedMinMaxExtra(unsigned int initialPos, unsigned int numPos, V &unifiedMinVec, V &unifiedMaxVec) const
Finds the minimum and the maximum values of the unified sequence, considering numPos positions starti...
void subGaussian1dKde(unsigned int initialPos, const V &scaleVec, const std::vector< V * > &evalParamVecs, std::vector< V * > &densityVecs) const
Gaussian kernel for the KDE estimate of the sub-sequence.
void copy(const BaseVectorSequence< V, M > &src)
Copies vector sequence src to this.
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
void subSampleVarianceExtra(unsigned int initialPos, unsigned int numPos, const V &meanVec, V &samVec) const
Finds the sample variance of the sub-sequence, considering numPos positions starting at position init...
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:78
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, V &minVec, V &maxVec) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
T unifiedInterQuantileRange(bool useOnlyInter0Comm, unsigned int initialPos) const
Returns the interquartile range of the values in the unified sequence.
void unifiedUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &unifiedCdfGrids, ArrayOfOneDTables< V, M > &unifiedCdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
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...
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates the autocorrelation via definition.
void subScalesForKde(unsigned int initialPos, const V &iqrVec, unsigned int kdeDimension, V &scaleVec) const
Selects the scales (bandwidth, scaleVec) for the kernel density estimation, considering only the sub-...
void unifiedHistogram(unsigned int initialPos, const V &unifiedMinVec, const V &unifiedMaxVec, std::vector< V * > &unifiedCentersForAllBins, std::vector< V * > &unifiedQuanttsForAllBins) const
Calculates the histogram of the unified sequence.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
void unifiedGaussian1dKde(unsigned int initialPos, const V &unifiedScaleVec, const std::vector< V * > &unifiedEvalParamVecs, std::vector< V * > &unifiedDensityVecs) const
Gaussian kernel for the KDE estimate of the unified sequence.
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.
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
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.
Class to accommodate arrays of one-dimensional tables.
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:90
void unifiedSampleVarianceExtra(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedSamVec) const
Finds the sample variance of the unified sequence, considering numPos positions starting at position ...
#define RawValue_MPI_MAX
Definition: MpiComm.h:51
Base class for handling vector and array samples (sequence of vectors or arrays). ...
#define RawValue_MPI_MIN
Definition: MpiComm.h:50
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.
void subInterQuantileRange(unsigned int initialPos, V &iqrVec) const
Returns the interquartile range of the values in the sub-sequence.
#define RawValue_MPI_SUM
Definition: MpiComm.h:52
void unifiedSampleStd(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedStdVec) const
Finds the sample standard deviation of the unified sequence, considering numPos positions starting at...
void copy(const SequenceOfVectors< V, M > &src)
Copies vector sequence src to this.
A class representing a vector space.
Definition: VectorSet.h:46
void select(const std::vector< unsigned int > &idsOfUniquePositions)
TODO: It shall select positions in the sequence of vectors.
Class for handling vector samples (sequence of vectors).
~SequenceOfVectors()
Destructor.
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
void resetValues(unsigned int initialPos, unsigned int numPos)
Resets a total of numPos values of the sequence starting at position initialPos.
void unifiedMeanExtra(unsigned int initialPos, unsigned int numPos, V &unifiedMeanVec) const
Finds the mean value of the unified sequence, considering numPos positions starting at position initi...
void subPopulationVariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, V &popVec) const
Finds the population variance of the sub-sequence, considering numPos positions starting at position ...
SequenceOfVectors(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
std::vector< const V * > m_seq
Sequence of vectors.

Generated on Thu Apr 23 2015 19:30:54 for queso-0.52.0 by  doxygen 1.8.5