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