queso-0.51.1
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,2009,2010,2011,2012,2013 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 
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 ...
T unifiedInterQuantileRange(bool useOnlyInter0Comm, unsigned int initialPos) const
Returns the interquartile range of the values in the unified sequence.
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...
Class to accommodate arrays of one-dimensional grid.
void subUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &cdfGrids, ArrayOfOneDTables< V, M > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void unifiedInterQuantileRange(unsigned int initialPos, V &unifiedIqrVec) const
Returns the interquartile range of the values in the unified sequence.
void deleteStoredVectors()
Deletes all the stored vectors.
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...
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 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 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.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:75
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.
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...
double estimateConvBrooksGelman(unsigned int initialPos, unsigned int numPos) const
Estimates convergence rate using Brooks &amp; Gelman method.
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 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 UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:89
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...
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 getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
void select(const std::vector< unsigned int > &idsOfUniquePositions)
TODO: It shall select positions in the sequence of vectors.
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 extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, ScalarSequence< double > &scalarSeq) const
Extracts a sequence of scalars.
void autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
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 setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
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...
void copy(const SequenceOfVectors< V, M > &src)
Copies vector sequence src to this.
#define RawValue_MPI_SUM
Definition: MpiComm.h:52
#define RawValue_MPI_MAX
Definition: MpiComm.h:51
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
Definition: Defines.h:207
void unifiedUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &unifiedCdfGrids, ArrayOfOneDTables< V, M > &unifiedCdfValues) const
Uniformly samples from the CDF from the sub-sequence.
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...
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:78
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 * >::iterator seqVectorPositionIteratorTypedef
Struct for handling data input and output from files.
Definition: Environment.h:66
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 extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
std::vector< const V * > m_seq
Sequence of vectors.
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.
Base class for handling vector and array samples (sequence of vectors or arrays). ...
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...
T subInterQuantileRange(unsigned int initialPos) const
Returns the interquartile range of the values in the sub-sequence.
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 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.
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
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...
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 subInterQuantileRange(unsigned int initialPos, V &iqrVec) const
Returns the interquartile range of the values in 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...
#define RawValue_MPI_MIN
Definition: MpiComm.h:50
Class to accommodate arrays of one-dimensional tables.
SequenceOfVectors< V, M > & operator=(const SequenceOfVectors< V, M > &rhs)
Copies values from rhs to this.
A class representing a vector space.
Definition: VectorSet.h:46
void copy(const BaseVectorSequence< V, M > &src)
Copies vector sequence src to this.
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:90
Class for handling vector samples (sequence of vectors).
SequenceOfVectors(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates the autocorrelation via definition.
const int UQ_OK_RC
Definition: Defines.h:76
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:223
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2409
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 resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
double MiscGetEllapsedSeconds(struct timeval *timeval0)
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...
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
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.
~SequenceOfVectors()
Destructor.

Generated on Thu Apr 23 2015 19:26:15 for queso-0.51.1 by  doxygen 1.8.5