queso-0.54.0
SequenceOfVectors.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2015 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/SequenceOfVectors.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Default constructor -----------------------------
32 template <class V, class M>
34  const VectorSpace<V,M>& vectorSpace,
35  unsigned int subSequenceSize,
36  const std::string& name)
37  :
38  BaseVectorSequence<V,M>(vectorSpace,subSequenceSize,name),
39  m_seq (subSequenceSize,NULL)
40 #ifdef UQ_CODE_HAS_MONITORS
41  ,
42  m_subMeanMonitorPosSeq (NULL),
43  m_subMeanVecSeq (NULL),
44  m_subMeanCltStdSeq (NULL),
45  m_subMeanInter0MonitorPosSeq (NULL),
46  m_subMeanInter0Mean (NULL),
47  m_subMeanInter0Clt95 (NULL),
48  m_subMeanInter0Empirical90 (NULL),
49  m_subMeanInter0Min (NULL),
50  m_subMeanInter0Max (NULL),
51  m_unifiedMeanMonitorPosSeq (NULL),
52  m_unifiedMeanVecSeq (NULL),
53  m_unifiedMeanCltStdSeq (NULL)
54 #endif
55 {
56  //if (m_env.subDisplayFile()) {
57  // *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::constructor()"
58  // << std::endl;
59  //}
60 
61  //if (m_env.subDisplayFile()) {
62  // *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::constructor()"
63  // << std::endl;
64  //}
65 }
66 // Destructor ---------------------------------------
67 template <class V, class M>
69 {
70 #ifdef UQ_CODE_HAS_MONITORS
71  if (m_subMeanMonitorPosSeq) delete m_subMeanMonitorPosSeq;
72  if (m_subMeanVecSeq ) delete m_subMeanVecSeq;
73  if (m_subMeanCltStdSeq ) delete m_subMeanCltStdSeq;
74 
75  if (m_subMeanInter0MonitorPosSeq) delete m_subMeanInter0MonitorPosSeq;
76  if (m_subMeanInter0Mean ) delete m_subMeanInter0Mean;
77  if (m_subMeanInter0Clt95 ) delete m_subMeanInter0Clt95;
78  if (m_subMeanInter0Empirical90 ) delete m_subMeanInter0Empirical90;
79  if (m_subMeanInter0Min ) delete m_subMeanInter0Min;
80  if (m_subMeanInter0Max ) delete m_subMeanInter0Max;
81 
82  if (m_unifiedMeanMonitorPosSeq) delete m_unifiedMeanMonitorPosSeq;
83  if (m_unifiedMeanVecSeq ) delete m_unifiedMeanVecSeq;
84  if (m_unifiedMeanCltStdSeq ) delete m_unifiedMeanCltStdSeq;
85 #endif
86 
87  for (unsigned int i = 0; i < (unsigned int) m_seq.size(); ++i) {
88  if (m_seq[i]) delete m_seq[i];
89  }
90 }
91 // Set methods --------------------------------------
92 template <class V, class M>
95 {
96  this->copy(rhs);
97  return *this;
98 }
99 
100 // Sequence methods ---------------------------------
101 template <class V, class M>
102 unsigned int
104 {
105  return m_seq.size();
106 }
107 //---------------------------------------------------
108 template <class V, class M>
109 void
110 SequenceOfVectors<V,M>::resizeSequence(unsigned int newSubSequenceSize)
111 {
112  if (newSubSequenceSize != this->subSequenceSize()) {
113  if (newSubSequenceSize < this->subSequenceSize()) {
114  this->resetValues(newSubSequenceSize,this->subSequenceSize()-newSubSequenceSize);
115  }
116  m_seq.resize(newSubSequenceSize,NULL);
117  std::vector<const V*>(m_seq).swap(m_seq);
119  }
120 
121  return;
122 }
123 //---------------------------------------------------
124 template <class V, class M>
125 void
126 SequenceOfVectors<V,M>::resetValues(unsigned int initialPos, unsigned int numPos)
127 {
128  bool bRC = ((initialPos < this->subSequenceSize()) &&
129  (0 < numPos ) &&
130  ((initialPos+numPos) <= this->subSequenceSize()));
131  if ((bRC == false) && (m_env.subDisplayFile())) {
132  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::resetValues()"
133  << ", initialPos = " << initialPos
134  << ", this->subSequenceSize() = " << this->subSequenceSize()
135  << ", numPos = " << numPos
136  << std::endl;
137  }
138  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 
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 
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 
221 
222  return;
223 }
224 //---------------------------------------------------
225 template <class V, class M>
226 void
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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  fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
1322  this->subWriteContents(initialPos,
1323  numPos,
1324  *filePtrSet.ofsVar,
1325  fileType);
1326  }
1327  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1328  queso_error_msg("hdf file type not supported yet");
1329  }
1330  else {
1331  queso_error_msg("invalid file type");
1332  }
1333 
1334  return;
1335 }
1336 //---------------------------------------------------
1337 template <class V, class M>
1338 void
1340  unsigned int initialPos,
1341  unsigned int numPos,
1342  std::ofstream& ofs,
1343  const std::string& fileType) const
1344 {
1345  queso_require_less_equal_msg((initialPos+numPos), this->subSequenceSize(), "invalid routine input parameters");
1346 
1347  if (initialPos == 0) {
1348  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1349  // Need sub matlab header here since this is subWriteContents
1350  this->writeSubMatlabHeader(ofs,
1351  this->subSequenceSize(),
1352  this->vectorSizeLocal());
1353  }
1354  else if (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT) {
1355  this->writeTxtHeader(ofs,
1356  this->subSequenceSize(),
1357  this->vectorSizeLocal());
1358  }
1359  }
1360 
1361  for (unsigned int j = initialPos; j < initialPos+numPos; ++j) {
1362  bool savedVectorPrintScientific = m_seq[j]->getPrintScientific();
1363  bool savedVectorPrintState = m_seq[j]->getPrintHorizontally();
1364  m_seq[j]->setPrintScientific (true);
1365  m_seq[j]->setPrintHorizontally(true);
1366 
1367  ofs << *(m_seq[j])
1368  << std::endl;
1369 
1370  m_seq[j]->setPrintHorizontally(savedVectorPrintState);
1371  m_seq[j]->setPrintScientific (savedVectorPrintScientific);
1372  }
1373 
1374  // Write Matlab-specific ending if desired
1375  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) &&
1376  ((initialPos + numPos) == this->subSequenceSize())) {
1377  ofs << "];\n";
1378  }
1379 }
1380 
1381 template <class V, class M>
1382 void
1384  double sequenceSize, double vectorSizeLocal) const
1385 {
1386  ofs << m_name << "_sub" << m_env.subIdString() << " = zeros(" << sequenceSize
1387  << "," << vectorSizeLocal
1388  << ");"
1389  << std::endl;
1390  ofs << m_name << "_sub" << m_env.subIdString() << " = [";
1391 }
1392 
1393 template <class V, class M>
1394 void
1396  double sequenceSize, double vectorSizeLocal) const
1397 {
1398  ofs << m_name << "_unified" << " = zeros(" << sequenceSize
1399  << "," << vectorSizeLocal
1400  << ");"
1401  << std::endl;
1402  ofs<< m_name << "_unified" << " = [";
1403 }
1404 
1405 template <class V, class M>
1406 void
1408  double sequenceSize, double vectorSizeLocal) const
1409 {
1410  ofs << sequenceSize << " " << vectorSizeLocal
1411  << std::endl;
1412 }
1413 
1414 //---------------------------------------------------
1415 template <class V, class M>
1416 void
1418  const std::string& fileName,
1419  const std::string& inputFileType) const
1420 {
1421  std::string fileType(inputFileType);
1422 #ifdef QUESO_HAS_HDF5
1423  // Do nothing
1424 #else
1425  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1426  if (m_env.subDisplayFile()) {
1427  *m_env.subDisplayFile() << "WARNING in SequenceOfVectors<V,M>::unifiedWriteContents()"
1428  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1429  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1430  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1431  << "' instead..."
1432  << std::endl;
1433  }
1434  if (m_env.subRank() == 0) {
1435  std::cerr << "WARNING in SequenceOfVectors<V,M>::unifiedWriteContents()"
1436  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1437  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1438  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1439  << "' instead..."
1440  << std::endl;
1441  }
1443  }
1444 #endif
1445 
1446  // All processors in 'fullComm' should call this routine...
1447 
1448  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ... // prudenci-2011-01-17
1449  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1450  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedWriteContents()"
1451  << ": worldRank " << m_env.worldRank()
1452  << ", fullRank " << m_env.fullRank()
1453  << ", subEnvironment " << m_env.subId()
1454  << ", subRank " << m_env.subRank()
1455  << ", inter0Rank " << m_env.inter0Rank()
1456  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1457  << ", fileName = " << fileName
1458  << std::endl;
1459  }
1460 
1461  if (m_env.inter0Rank() >= 0) {
1462  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) {
1463  if (m_env.inter0Rank() == (int) r) {
1464  // My turn
1465  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1466  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1467  << ": worldRank " << m_env.worldRank()
1468  << ", fullRank " << m_env.fullRank()
1469  << ", subEnvironment " << m_env.subId()
1470  << ", subRank " << m_env.subRank()
1471  << ", inter0Rank " << m_env.inter0Rank()
1472  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1473  << ", fileName = " << fileName
1474  << ", about to open file for r = " << r
1475  << std::endl;
1476  }
1477 
1478  // bool writeOver = (r == 0);
1479  bool writeOver = false; // A 'true' causes problems when the user chooses (via options
1480  // in the input file) to use just one file for all outputs.
1481  FilePtrSetStruct unifiedFilePtrSet;
1482  if (m_env.openUnifiedOutputFile(fileName,
1483  fileType, // "m or hdf"
1484  writeOver,
1485  unifiedFilePtrSet)) {
1486  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) { // 2013-02-23
1487  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1488  << ": worldRank " << m_env.worldRank()
1489  << ", fullRank " << m_env.fullRank()
1490  << ", subEnvironment " << m_env.subId()
1491  << ", subRank " << m_env.subRank()
1492  << ", inter0Rank " << m_env.inter0Rank()
1493  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1494  << ", fileName = " << fileName
1495  << ", just opened file for r = " << r
1496  << std::endl;
1497  }
1498 
1499  unsigned int chainSize = this->subSequenceSize();
1500  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
1501  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
1502  if (r == 0) {
1503  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1504  // Need unified matlab header here since this is unifiedWriteContents
1505  writeUnifiedMatlabHeader(*unifiedFilePtrSet.ofsVar,
1506  this->subSequenceSize()*m_env.inter0Comm().NumProc(),
1507  this->vectorSizeLocal());
1508  }
1509  else { // If we get here it's definitely a txt file not matlab
1510  writeTxtHeader(*unifiedFilePtrSet.ofsVar,
1511  this->subSequenceSize()*m_env.inter0Comm().NumProc(),
1512  this->vectorSizeLocal());
1513  }
1514  }
1515 
1516  for (unsigned int j = 0; j < chainSize; ++j) { // 2013-02-23
1517  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): m_seq[" << j << "] = " << m_seq[j]
1518  // << std::endl;
1519  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): &(m_seq[" << j << "].map()) = " << &(m_seq[j]->map())
1520  // << std::endl;
1521  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): (m_seq[" << j << "].map().NumMyElements = " << m_seq[j]->map().NumMyElements()
1522  // << std::endl;
1523  //V tmpVec(*(m_seq[j]));
1524  //std::cout << "*(m_seq[" << j << "]) = " << tmpVec
1525  // << std::endl;
1526  //std::cout << "*(m_seq[" << j << "]) = " << *(m_seq[j])
1527  // << std::endl;
1528 
1529  bool savedVectorPrintScientific = m_seq[j]->getPrintScientific();
1530  bool savedVectorPrintState = m_seq[j]->getPrintHorizontally();
1531  m_seq[j]->setPrintScientific (true);
1532  m_seq[j]->setPrintHorizontally(true);
1533 
1534  *unifiedFilePtrSet.ofsVar << *(m_seq[j])
1535  << std::endl;
1536 
1537  m_seq[j]->setPrintHorizontally(savedVectorPrintState);
1538  m_seq[j]->setPrintScientific (savedVectorPrintScientific);
1539  }
1540  }
1541 #ifdef QUESO_HAS_HDF5
1542  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1543  unsigned int numParams = m_vectorSpace.dimLocal();
1544  if (r == 0) {
1545  hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
1546  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data type created" << std::endl;
1547  hsize_t dimsf[2];
1548  dimsf[0] = numParams;
1549  dimsf[1] = chainSize;
1550  hid_t dataspace = H5Screate_simple(2, dimsf, NULL); // HDF5_rank = 2
1551  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data space created" << std::endl;
1552  hid_t dataset = H5Dcreate2(unifiedFilePtrSet.h5Var,
1553  "seq_of_vectors",
1554  datatype,
1555  dataspace,
1556  H5P_DEFAULT, // Link creation property list
1557  H5P_DEFAULT, // Dataset creation property list
1558  H5P_DEFAULT); // Dataset access property list
1559  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data set created" << std::endl;
1560 
1561  struct timeval timevalBegin;
1562  int iRC = UQ_OK_RC;
1563  iRC = gettimeofday(&timevalBegin,NULL);
1564  if (iRC) {}; // just to remover compiler warning
1565 
1566  //double* dataOut[numParams]; // avoid compiler warning
1567  std::vector<double*> dataOut((size_t) numParams,NULL);
1568  dataOut[0] = (double*) malloc(numParams*chainSize*sizeof(double));
1569  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
1570  dataOut[i] = dataOut[i-1] + chainSize; // Yes, just 'chainSize', not 'chainSize*sizeof(double)'
1571  }
1572  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, memory allocated" << std::endl;
1573  for (unsigned int j = 0; j < chainSize; ++j) {
1574  V tmpVec(*(m_seq[j]));
1575  for (unsigned int i = 0; i < numParams; ++i) {
1576  dataOut[i][j] = tmpVec[i];
1577  }
1578  }
1579  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, memory filled" << std::endl;
1580 
1581  herr_t status;
1582  status = H5Dwrite(dataset,
1583  H5T_NATIVE_DOUBLE,
1584  H5S_ALL,
1585  H5S_ALL,
1586  H5P_DEFAULT,
1587  (void*) dataOut[0]);
1588  if (status) {}; // just to remover compiler warning
1589  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data written" << std::endl;
1590 
1591  double writeTime = MiscGetEllapsedSeconds(&timevalBegin);
1592  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
1593  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1594  << ": worldRank " << m_env.worldRank()
1595  << ", fullRank " << m_env.fullRank()
1596  << ", subEnvironment " << m_env.subId()
1597  << ", subRank " << m_env.subRank()
1598  << ", inter0Rank " << m_env.inter0Rank()
1599  << ", fileName = " << fileName
1600  << ", numParams = " << numParams
1601  << ", chainSize = " << chainSize
1602  << ", writeTime = " << writeTime << " seconds"
1603  << std::endl;
1604  }
1605 
1606  H5Dclose(dataset);
1607  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data set closed" << std::endl;
1608  H5Sclose(dataspace);
1609  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data space closed" << std::endl;
1610  H5Tclose(datatype);
1611  //std::cout << "In SequenceOfVectors<V,M>::unifiedWriteContents(): h5 case, data type closed" << std::endl;
1612  //free(dataOut[0]); // related to the changes above for compiler warning
1613  for (unsigned int tmpIndex = 0; tmpIndex < dataOut.size(); tmpIndex++) {
1614  free (dataOut[tmpIndex]);
1615  }
1616  }
1617  else {
1618  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
1619  }
1620  }
1621 #endif
1622  else {
1623  queso_error_msg("invalid file type");
1624  }
1625 
1626  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1627  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1628  << ": worldRank " << m_env.worldRank()
1629  << ", fullRank " << m_env.fullRank()
1630  << ", subEnvironment " << m_env.subId()
1631  << ", subRank " << m_env.subRank()
1632  << ", inter0Rank " << m_env.inter0Rank()
1633  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1634  << ", fileName = " << fileName
1635  << ", about to close file for r = " << r
1636  << std::endl;
1637  }
1638 
1639  m_env.closeFile(unifiedFilePtrSet,fileType);
1640 
1641  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1642  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedWriteContents()"
1643  << ": worldRank " << m_env.worldRank()
1644  << ", fullRank " << m_env.fullRank()
1645  << ", subEnvironment " << m_env.subId()
1646  << ", subRank " << m_env.subRank()
1647  << ", inter0Rank " << m_env.inter0Rank()
1648  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1649  << ", fileName = " << fileName
1650  << ", just closed file for r = " << r
1651  << std::endl;
1652  }
1653  } // if (m_env.openUnifiedOutputFile())
1654  } // if (m_env.inter0Rank() == (int) r)
1655  m_env.inter0Comm().Barrier();
1656  } // for r
1657 
1658  if (m_env.inter0Rank() == 0) {
1659  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
1660  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
1661  FilePtrSetStruct unifiedFilePtrSet;
1662  if (m_env.openUnifiedOutputFile(fileName,
1663  fileType,
1664  false, // Yes, 'writeOver = false' in order to close the array for matlab
1665  unifiedFilePtrSet)) {
1666 
1667  if (fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) {
1668  *unifiedFilePtrSet.ofsVar << "];\n";
1669  }
1670 
1671  m_env.closeFile(unifiedFilePtrSet,fileType);
1672  }
1673  }
1674  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1675  // Do nothing
1676  }
1677  else {
1678  queso_error_msg("invalid file type");
1679  }
1680  }
1681  } // if (m_env.inter0Rank() >= 0)
1682 
1683  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 10)) {
1684  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedWriteContents()"
1685  << ", fileName = " << fileName
1686  << std::endl;
1687  }
1688  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ... // prudenci-2011-01-17
1689 
1690  return;
1691 }
1692 //---------------------------------------------------
1693 template <class V, class M>
1694 void
1696  const std::string& fileName,
1697  const std::string& inputFileType,
1698  const unsigned int subReadSize)
1699 {
1700  std::string fileType(inputFileType);
1701 #ifdef QUESO_HAS_HDF5
1702  // Do nothing
1703 #else
1704  if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1705  if (m_env.subDisplayFile()) {
1706  *m_env.subDisplayFile() << "WARNING in SequenceOfVectors<V,M>::unifiedReadContents()"
1707  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1708  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1709  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1710  << "' instead..."
1711  << std::endl;
1712  }
1713  if (m_env.subRank() == 0) {
1714  std::cerr << "WARNING in SequenceOfVectors<V,M>::unifiedReadContents()"
1715  << ": file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1716  << "' has been requested, but this QUESO library has not been built with 'hdf5'"
1717  << ". Code will therefore process the file format '" << UQ_FILE_EXTENSION_FOR_HDF_FORMAT
1718  << "' instead..."
1719  << std::endl;
1720  }
1722  }
1723 #endif
1724 
1725  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1726  if (m_env.subDisplayFile()) {
1727  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::unifiedReadContents()"
1728  << ": worldRank " << m_env.worldRank()
1729  << ", fullRank " << m_env.fullRank()
1730  << ", subEnvironment " << m_env.subId()
1731  << ", subRank " << m_env.subRank()
1732  << ", inter0Rank " << m_env.inter0Rank()
1733  //<< ", m_env.inter0Comm().NumProc() = " << m_env.inter0Comm().NumProc()
1734  << ", fileName = " << fileName
1735  << ", subReadSize = " << subReadSize
1736  //<< ", unifiedReadSize = " << unifiedReadSize
1737  << std::endl;
1738  }
1739 
1740  this->resizeSequence(subReadSize);
1741 
1742  if (m_env.inter0Rank() >= 0) {
1743  double unifiedReadSize = subReadSize*m_env.inter0Comm().NumProc();
1744 
1745  // In the logic below, the id of a line' begins with value 0 (zero)
1746  unsigned int idOfMyFirstLine = 1 + m_env.inter0Rank()*subReadSize;
1747  unsigned int idOfMyLastLine = (1 + m_env.inter0Rank())*subReadSize;
1748  unsigned int numParams = this->vectorSizeLocal();
1749 
1750  for (unsigned int r = 0; r < (unsigned int) m_env.inter0Comm().NumProc(); ++r) { // "m or hdf"
1751  if (m_env.inter0Rank() == (int) r) {
1752  // My turn
1753  FilePtrSetStruct unifiedFilePtrSet;
1754  if (m_env.openUnifiedInputFile(fileName,
1755  fileType,
1756  unifiedFilePtrSet)) {
1757  if ((fileType == UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT) ||
1758  (fileType == UQ_FILE_EXTENSION_FOR_TXT_FORMAT)) {
1759  if (r == 0) {
1760  // Read number of chain positions in the file by taking care of the first line,
1761  // which resembles something like 'variable_name = zeros(n_positions,m_params);'
1762  std::string tmpString;
1763 
1764  // Read 'variable name' string
1765  *unifiedFilePtrSet.ifsVar >> tmpString;
1766  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1767 
1768  // Read '=' sign
1769  *unifiedFilePtrSet.ifsVar >> tmpString;
1770  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1771  queso_require_equal_to_msg(tmpString, "=", "string should be the '=' sign");
1772 
1773  // Read 'zeros(n_positions,n_params)' string
1774  *unifiedFilePtrSet.ifsVar >> tmpString;
1775  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1776  unsigned int posInTmpString = 6;
1777 
1778  // Isolate 'n_positions' in a string
1779  //char nPositionsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
1780  std::string nPositionsString((size_t) (tmpString.size()-posInTmpString+1),' ');
1781  unsigned int posInPositionsString = 0;
1782  do {
1783  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
1784  nPositionsString[posInPositionsString++] = tmpString[posInTmpString++];
1785  } while (tmpString[posInTmpString] != ',');
1786  nPositionsString[posInPositionsString] = '\0';
1787 
1788  // Isolate 'n_params' in a string
1789  posInTmpString++; // Avoid reading ',' char
1790  //char nParamsString[tmpString.size()-posInTmpString+1]; // avoid compiler warning
1791  std::string nParamsString((size_t) (tmpString.size()-posInTmpString+1),' ');
1792  unsigned int posInParamsString = 0;
1793  do {
1794  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
1795  nParamsString[posInParamsString++] = tmpString[posInTmpString++];
1796  } while (tmpString[posInTmpString] != ')');
1797  nParamsString[posInParamsString] = '\0';
1798 
1799  // Convert 'n_positions' and 'n_params' strings to numbers
1800  unsigned int sizeOfChainInFile = (unsigned int) strtod(nPositionsString.c_str(),NULL);
1801  unsigned int numParamsInFile = (unsigned int) strtod(nParamsString.c_str(), NULL);
1802  if (m_env.subDisplayFile()) {
1803  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedReadContents()"
1804  << ": worldRank " << m_env.worldRank()
1805  << ", fullRank " << m_env.fullRank()
1806  << ", sizeOfChainInFile = " << sizeOfChainInFile
1807  << ", numParamsInFile = " << numParamsInFile
1808  << std::endl;
1809  }
1810 
1811  // Check if [size of chain in file] >= [requested unified sequence size]
1812  queso_require_greater_equal_msg(sizeOfChainInFile, unifiedReadSize, "size of chain in file is not big enough");
1813 
1814  // Check if [num params in file] == [num params in current chain]
1815  queso_require_equal_to_msg(numParamsInFile, numParams, "number of parameters of chain in file is different than number of parameters in this chain object");
1816  } // if (r == 0)
1817 
1818  // Code common to any core in 'inter0Comm', including core of rank 0
1819  unsigned int maxCharsPerLine = 64*numParams; // Up to about 60 characters to represent each parameter value
1820 
1821  unsigned int lineId = 0;
1822  while (lineId < idOfMyFirstLine) {
1823  unifiedFilePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
1824  lineId++;
1825  };
1826 
1827  if (r == 0) {
1828  // Take care of initial part of the first data line,
1829  // which resembles something like 'variable_name = [value1 value2 ...'
1830  std::string tmpString;
1831 
1832  // Read 'variable name' string
1833  *unifiedFilePtrSet.ifsVar >> tmpString;
1834  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1835 
1836  // Read '=' sign
1837  *unifiedFilePtrSet.ifsVar >> tmpString;
1838  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1839  queso_require_equal_to_msg(tmpString, "=", "in core 0, string should be the '=' sign");
1840 
1841  // Take into account the ' [' portion
1842  std::streampos tmpPos = unifiedFilePtrSet.ifsVar->tellg();
1843  unifiedFilePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
1844  }
1845 
1846  V tmpVec(m_vectorSpace.zeroVector());
1847  while (lineId <= idOfMyLastLine) {
1848  for (unsigned int i = 0; i < numParams; ++i) {
1849  *unifiedFilePtrSet.ifsVar >> tmpVec[i];
1850  }
1851  this->setPositionValues(lineId - idOfMyFirstLine, tmpVec);
1852  lineId++;
1853  };
1854  }
1855 #ifdef QUESO_HAS_HDF5
1856  else if (fileType == UQ_FILE_EXTENSION_FOR_HDF_FORMAT) {
1857  if (r == 0) {
1858  hid_t dataset = H5Dopen2(unifiedFilePtrSet.h5Var,
1859  "seq_of_vectors",
1860  H5P_DEFAULT); // Dataset access property list
1861  hid_t datatype = H5Dget_type(dataset);
1862  H5T_class_t t_class = H5Tget_class(datatype);
1863  queso_require_equal_to_msg(t_class, H5T_FLOAT, "t_class is not H5T_DOUBLE");
1864  hid_t dataspace = H5Dget_space(dataset);
1865  int rank = H5Sget_simple_extent_ndims(dataspace);
1866  queso_require_equal_to_msg(rank, 2, "hdf rank is not 2");
1867  hsize_t dims_in[2];
1868  int status_n;
1869  status_n = H5Sget_simple_extent_dims(dataspace, dims_in, NULL);
1870  if (status_n) {}; // just to remover compiler warning
1871  //std::cout << "In SequenceOfVectors<V,M>::unifiedReadContents()"
1872  // << ": dims_in[0] = " << dims_in[0]
1873  // << ", dims_in[1] = " << dims_in[1]
1874  // << std::endl;
1875  queso_require_equal_to_msg(dims_in[0], numParams, "dims_in[0] is not equal to 'numParams'");
1876  queso_require_greater_equal_msg(dims_in[1], subReadSize, "dims_in[1] is smaller that requested 'subReadSize'");
1877 
1878  struct timeval timevalBegin;
1879  int iRC = UQ_OK_RC;
1880  iRC = gettimeofday(&timevalBegin,NULL);
1881  if (iRC) {}; // just to remover compiler warning
1882 
1883  unsigned int chainSizeIn = (unsigned int) dims_in[1];
1884  //double* dataIn[numParams]; // avoid compiler warning
1885  std::vector<double*> dataIn((size_t) numParams,NULL);
1886  dataIn[0] = (double*) malloc(numParams*chainSizeIn*sizeof(double));
1887  for (unsigned int i = 1; i < numParams; ++i) { // Yes, from '1'
1888  dataIn[i] = dataIn[i-1] + chainSizeIn; // Yes, just 'chainSizeIn', not 'chainSizeIn*sizeof(double)'
1889  }
1890  //std::cout << "In SequenceOfVectors<V,M>::unifiedReadContents(): h5 case, memory allocated" << std::endl;
1891  herr_t status;
1892  status = H5Dread(dataset,
1893  H5T_NATIVE_DOUBLE,
1894  H5S_ALL,
1895  dataspace,
1896  H5P_DEFAULT,
1897  dataIn[0]);
1898  if (status) {}; // just to remover compiler warning
1899  //std::cout << "In SequenceOfVectors<V,M>::unifiedReadContents(): h5 case, data read" << std::endl;
1900  V tmpVec(m_vectorSpace.zeroVector());
1901  for (unsigned int j = 0; j < subReadSize; ++j) { // Yes, 'subReadSize', not 'chainSizeIn'
1902  for (unsigned int i = 0; i < numParams; ++i) {
1903  tmpVec[i] = dataIn[i][j];
1904  }
1905  this->setPositionValues(j, tmpVec);
1906  }
1907 
1908  double readTime = MiscGetEllapsedSeconds(&timevalBegin);
1909  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
1910  *m_env.subDisplayFile() << "In SequenceOfVectors<V,M>::unifiedReadContents()"
1911  << ": worldRank " << m_env.worldRank()
1912  << ", fullRank " << m_env.fullRank()
1913  << ", subEnvironment " << m_env.subId()
1914  << ", subRank " << m_env.subRank()
1915  << ", inter0Rank " << m_env.inter0Rank()
1916  << ", fileName = " << fileName
1917  << ", numParams = " << numParams
1918  << ", chainSizeIn = " << chainSizeIn
1919  << ", subReadSize = " << subReadSize
1920  << ", readTime = " << readTime << " seconds"
1921  << std::endl;
1922  }
1923 
1924  H5Sclose(dataspace);
1925  H5Tclose(datatype);
1926  H5Dclose(dataset);
1927  //free(dataIn[0]); // related to the changes above for compiler warning
1928  for (unsigned int tmpIndex = 0; tmpIndex < dataIn.size(); tmpIndex++) {
1929  free (dataIn[tmpIndex]);
1930  }
1931  }
1932  else {
1933  queso_error_msg("hdf file type not supported for multiple sub-environments yet");
1934  }
1935  }
1936 #endif
1937  else {
1938  queso_error_msg("invalid file type");
1939  }
1940  m_env.closeFile(unifiedFilePtrSet,fileType);
1941  } // if (m_env.openUnifiedInputFile())
1942  } // if (m_env.inter0Rank() == (int) r)
1943  m_env.inter0Comm().Barrier();
1944  } // for r
1945  } // if (m_env.inter0Rank() >= 0)
1946  else {
1947  V tmpVec(m_vectorSpace.zeroVector());
1948  for (unsigned int i = 1; i < subReadSize; ++i) {
1949  this->setPositionValues(i,tmpVec);
1950  }
1951  }
1952 
1953  if (m_env.subDisplayFile()) {
1954  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::unifiedReadContents()"
1955  << ", fileName = " << fileName
1956  << std::endl;
1957  }
1958  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
1959 
1960  return;
1961 }
1962 //---------------------------------------------------
1963 template <class V, class M>
1964 void
1965 SequenceOfVectors<V,M>::select(const std::vector<unsigned int>& /* idsOfUniquePositions */)
1966 {
1967  queso_error_msg("Code is not complete yet");
1968 
1969  return;
1970 }
1971 //---------------------------------------------------
1972 template <class V, class M>
1973 void
1975  unsigned int initialPos,
1976  unsigned int spacing)
1977 {
1978  if (m_env.subDisplayFile()) {
1979  *m_env.subDisplayFile() << "Entering SequenceOfVectors<V,M>::filter()"
1980  << ": initialPos = " << initialPos
1981  << ", spacing = " << spacing
1982  << ", subSequenceSize = " << this->subSequenceSize()
1983  << std::endl;
1984  }
1985 
1986  unsigned int i = 0;
1987  unsigned int j = initialPos;
1988  unsigned int originalSubSequenceSize = this->subSequenceSize();
1989  while (j < originalSubSequenceSize) {
1990  if (i != j) {
1991  //*m_env.subDisplayFile() << i << "--" << j << " ";
1992  delete m_seq[i];
1993  m_seq[i] = new V(*(m_seq[j]));
1994  }
1995  i++;
1996  j += spacing;
1997  }
1998 
1999  this->resetValues(i,originalSubSequenceSize-i);
2000  this->resizeSequence(i);
2001 
2002  if (m_env.subDisplayFile()) {
2003  *m_env.subDisplayFile() << "Leaving SequenceOfVectors<V,M>::filter()"
2004  << ": initialPos = " << initialPos
2005  << ", spacing = " << spacing
2006  << ", subSequenceSize = " << this->subSequenceSize()
2007  << std::endl;
2008  }
2009 
2010  return;
2011 }
2012 //---------------------------------------------------
2013 template <class V, class M>
2014 double
2016  unsigned int initialPos,
2017  unsigned int numPos) const
2018 {
2019  // This method requires *at least* two sequences. Error if there is only one.
2020  queso_require_greater_equal_msg(m_env.numSubEnvironments(), 2, "At least two sequences required for Brooks-Gelman convergence test.");
2021 
2022  // TODO: Need special case for 1-dimensional parameter space.
2023 
2024  // Initialize with garbage to give the user a clue something is funky.
2025  double convMeasure = -1.0;
2026 
2027  // We only do the work on the subenvironment where the sequence data
2028  // is stashed.
2029  if( m_env.inter0Rank() >= 0 )
2030  {
2031  // Sanity Checks
2032 
2033  // REMEMBER: \psi is a *vector* of parameters
2034  // Get quantities we will use several times
2035  V psi_j_dot = m_vectorSpace.zeroVector();
2036  V psi_dot_dot = m_vectorSpace.zeroVector();
2037  V work = m_vectorSpace.zeroVector();
2038 
2039  // m = number of chains > 1
2040  // n = number of steps for which we are computing the metric
2041  int m = m_env.numSubEnvironments();
2042  int n = numPos;
2043 
2044  this->subMeanExtra ( initialPos, numPos, psi_j_dot );
2045  this->unifiedMeanExtra( initialPos, numPos, psi_dot_dot );
2046 
2047 #if 0
2048  std::cout << "psi_j_dot = " << psi_j_dot << std::endl;
2049  std::cout << "psi_dot_dot = " << psi_dot_dot << std::endl;
2050 #endif
2051 
2052  /* W = \frac{1}{m*(n-1)}*\sum_{j=1}^m \sum{t=1}^n
2053  (\psi_{jt} - \overline{\psi_{j\cdot}})*(\psi_{jt} - \overline{\psi_{j\cdot}})^T
2054  This corresponds to the "within-sequence" covariance matrix. */
2055  M* W_local = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2056  M* W = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2057  V psi_j_t = m_vectorSpace.zeroVector();
2058 
2059  // Sum within the chain
2060  for( unsigned int t = initialPos; t < initialPos+numPos; ++t )
2061  {
2062  psi_j_t = *(m_seq[t]);
2063 
2064  work = psi_j_t - psi_j_dot;
2065 
2066  (*W_local) += matrixProduct( work, work );
2067  }
2068 
2069  // Now do the sum over the chains
2070  // W will be available on all inter0 processors
2071  W_local->mpiSum( m_env.inter0Comm(), (*W) );
2072 
2073  (*W) = 1.0/(double(m)*(double(n)-1.0)) * (*W);
2074 
2075 #if 0
2076  std::cout << "n, m = " << n << ", " << m << std::endl;
2077  std::cout << "W_local = " << *W_local << std::endl;
2078  std::cout << "W = " << *W << std::endl;
2079 #endif
2080 
2081  // Need to delete pointers to temporary covariance matrices
2082  delete W_local;
2083 
2084  /* B/n = \frac{1}{m-1}\sum_{j=1}^m
2085  (\overline{\psi_{j\cdot}} - \overline{\psi_{\cdot \cdot}})*
2086  (\overline{\psi_{j\cdot}} - \overline{\psi_{\cdot \cdot}})^T
2087  This corresponds to the "between-sequence" covariance matrix. */
2088  M* B_over_n_local = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2089  M* B_over_n = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2090 
2091  work = psi_j_dot - psi_dot_dot;
2092  (*B_over_n_local) = matrixProduct( work, work );
2093 
2094  B_over_n_local->mpiSum( m_env.inter0Comm(), (*B_over_n) );
2095 
2096  // Need to delete pointers to temporary covariance matrices
2097  delete B_over_n_local;
2098 
2099  (*B_over_n) = 1.0/(double(m)-1.0) * (*B_over_n);
2100 
2101 #if 0
2102  std::cout << "B_over_n = " << *B_over_n << std::endl;
2103 #endif
2104 
2105 
2106  /* R_p = (n-1)/n + (m+1)/m * \lambda
2107  \lambda = largest eigenvalue of W^{-1}*B/n */
2108  M* A = m_vectorSpace.newDiagMatrix( m_vectorSpace.zeroVector() );
2109 
2110  W->invertMultiply( *B_over_n, *A );
2111 
2112 #if 0
2113  std::cout << "A = " << *A << std::endl;
2114  std::cout.flush();
2115 #endif
2116  // Need to delete pointers to temporary covariance matrices
2117  delete W;
2118  delete B_over_n;
2119 
2120  double eigenValue;
2121  V eigenVector = m_vectorSpace.zeroVector();
2122 
2123  A->largestEigen( eigenValue, eigenVector );
2124 
2125  // Need to delete pointers to temporary covariance matrices
2126  delete A;
2127 
2128  // Now, finally compute the final convMeasure
2129  convMeasure = (double(n)-1.0)/double(n) + (double(m)+1.0)/double(m)*eigenValue;
2130 
2131  } // End of check on inter0Rank
2132 
2133  //TODO: Do we need a Barrier here?
2134  //TODO: Error checking on MPI.
2135  //m_env.fullComm().Barrier(); // Dangerous to barrier on fullComm ...
2136 
2137  return convMeasure;
2138 }
2139 //---------------------------------------------------
2140 template <class V, class M>
2141 void
2143  unsigned int initialPos,
2144  unsigned int spacing,
2145  unsigned int numPos,
2146  unsigned int paramId,
2147  ScalarSequence<double>& scalarSeq) const
2148 {
2149  scalarSeq.resizeSequence(numPos);
2150  if (spacing == 1) {
2151  for (unsigned int j = 0; j < numPos; ++j) {
2152  scalarSeq[j] = (*(m_seq[initialPos+j ]))[paramId];
2153  }
2154  }
2155  else {
2156  for (unsigned int j = 0; j < numPos; ++j) {
2157  scalarSeq[j] = (*(m_seq[initialPos+j*spacing]))[paramId];
2158  }
2159  }
2160 
2161  return;
2162 }
2163 // Private methods ------------------------------------
2164 template <class V, class M>
2165 void
2167 {
2169  for (unsigned int i = 0; i < (unsigned int) m_seq.size(); ++i) {
2170  if (m_seq[i]) {
2171  delete m_seq[i];
2172  m_seq[i] = NULL;
2173  }
2174  }
2175  m_seq.resize(src.subSequenceSize(),NULL);
2176  for (unsigned int i = 0; i < m_seq.size(); ++i) {
2177  m_seq[i] = new V(*(src.m_seq[i]));
2178  }
2179 
2180  return;
2181 }
2182 //---------------------------------------------------
2183 template <class V, class M>
2184 void
2186  unsigned int initialPos,
2187  unsigned int spacing,
2188  unsigned int numPos,
2189  unsigned int paramId,
2190  std::vector<double>& rawData) const
2191 {
2192  rawData.resize(numPos);
2193  if (spacing == 1) {
2194  for (unsigned int j = 0; j < numPos; ++j) {
2195  rawData[j] = (*(m_seq[initialPos+j ]))[paramId];
2196  }
2197  }
2198  else {
2199  for (unsigned int j = 0; j < numPos; ++j) {
2200  rawData[j] = (*(m_seq[initialPos+j*spacing]))[paramId];
2201  }
2202  }
2203 
2204  return;
2205 }
2206 
2207 // --------------------------------------------------
2208 // Methods conditionally available ------------------
2209 // --------------------------------------------------
2210 // --------------------------------------------------
2211 #ifdef UQ_SEQ_VEC_USES_OPERATOR
2212 template <class V, class M>
2213 const V*
2214 SequenceOfVectors<V,M>::operator[](unsigned int posId) const
2215 {
2216  queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
2217 
2218  return (const V*) (m_seq[posId]);
2219 }
2220 // --------------------------------------------------
2221 template <class V, class M>
2222 const V*&
2223 SequenceOfVectors<V,M>::operator[](unsigned int posId)
2224 {
2225  queso_require_less_msg(posId, this->subSequenceSize(), "posId > subSequenceSize()");
2226 
2227  return m_seq[posId];
2228 }
2229 #endif
2230 
2231 // --------------------------------------------------
2232 // --------------------------------------------------
2233 // --------------------------------------------------
2234 
2235 #ifdef UQ_ALSO_COMPUTE_MDFS_WITHOUT_KDE
2236 template <class V, class M>
2237 void
2238 SequenceOfVectors<V,M>::subUniformlySampledMdf(
2239  const V& numEvaluationPointsVec,
2240  ArrayOfOneDGrids <V,M>& mdfGrids,
2241  ArrayOfOneDTables<V,M>& mdfValues) const
2242 {
2243  V minDomainValues(m_vectorSpace.zeroVector());
2244  V maxDomainValues(m_vectorSpace.zeroVector());
2245 
2246  ScalarSequence<double> data(m_env,0,"");
2247 
2248  unsigned int numParams = this->vectorSizeLocal();
2249  for (unsigned int i = 0; i < numParams; ++i) {
2250  this->extractScalarSeq(0, // initialPos
2251  1, // spacing
2252  subSequenceSize(), // numPos
2253  i,
2254  data);
2255 
2256  std::vector<double> aMdf(0);
2257  data.subUniformlySampledMdf((unsigned int) numEvaluationPointsVec[i],
2258  minDomainValues[i],
2259  maxDomainValues[i],
2260  aMdf);
2261  mdfValues.setOneDTable(i,aMdf);
2262  }
2263 
2264  mdfGrids.setUniformGrids(numEvaluationPointsVec,
2265  minDomainValues,
2266  maxDomainValues);
2267 
2268  return;
2269 }
2270 #endif
2271 
2272 // --------------------------------------------------
2273 // --------------------------------------------------
2274 // --------------------------------------------------
2275 
2276 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2277 template <class V, class M>
2278 void
2279 SequenceOfVectors<V,M>::subMeanCltStd(
2280  unsigned int initialPos,
2281  unsigned int numPos,
2282  const V& meanVec,
2283  V& stdVec) const
2284 {
2285  bool bRC = ((initialPos < this->subSequenceSize()) &&
2286  (0 < numPos ) &&
2287  ((initialPos+numPos) <= this->subSequenceSize()) &&
2288  (this->vectorSizeLocal() == meanVec.sizeLocal() ) &&
2289  (this->vectorSizeLocal() == stdVec.sizeLocal() ));
2290  queso_require_msg(bRC, "invalid input data");
2291 
2292  ScalarSequence<double> data(m_env,0,"");
2293 
2294  unsigned int numParams = this->vectorSizeLocal();
2295  for (unsigned int i = 0; i < numParams; ++i) {
2296  this->extractScalarSeq(initialPos,
2297  1, // spacing
2298  numPos,
2299  i,
2300  data);
2301  stdVec[i] = data.subMeanCltStd(0,
2302  numPos,
2303  meanVec[i]);
2304  }
2305 
2306  return;
2307 }
2308 
2309 template <class V, class M>
2310 void
2311 SequenceOfVectors<V,M>::unifiedMeanCltStd(
2312  unsigned int initialPos,
2313  unsigned int numPos,
2314  const V& unifiedMeanVec,
2315  V& unifiedSamVec) const
2316 {
2317  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2318  (0 < numPos ) &&
2319  ((initialPos+numPos) <= this->subSequenceSize() ) &&
2320  (this->vectorSizeLocal() == unifiedMeanVec.sizeLocal()) &&
2321  (this->vectorSizeLocal() == unifiedSamVec.sizeLocal() ));
2322  queso_require_msg(bRC, "invalid input data");
2323 
2324  ScalarSequence<double> data(m_env,0,"");
2325 
2326  unsigned int numParams = this->vectorSizeLocal();
2327  for (unsigned int i = 0; i < numParams; ++i) {
2328  this->extractScalarSeq(initialPos,
2329  1, // spacing
2330  numPos,
2331  i,
2332  data);
2333  unifiedSamVec[i] = data.unifiedMeanCltStd(m_vectorSpace.numOfProcsForStorage() == 1,
2334  0,
2335  numPos,
2336  unifiedMeanVec[i]);
2337  }
2338 
2339  return;
2340 }
2341 //---------------------------------------------------
2342 template <class V, class M>
2343 void
2344 SequenceOfVectors<V,M>::bmm(
2345  unsigned int initialPos,
2346  unsigned int batchLength,
2347  V& bmmVec) const
2348 {
2349  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2350  (batchLength < (this->subSequenceSize()-initialPos)) &&
2351  (this->vectorSizeLocal() == bmmVec.sizeLocal() ));
2352  queso_require_msg(bRC, "invalid input data");
2353 
2354  ScalarSequence<double> data(m_env,0,"");
2355 
2356  unsigned int numParams = this->vectorSizeLocal();
2357  for (unsigned int i = 0; i < numParams; ++i) {
2358  this->extractScalarSeq(initialPos,
2359  1, // spacing
2360  this->subSequenceSize()-initialPos,
2361  i,
2362  data);
2363  bmmVec[i] = data.bmm(0,
2364  batchLength);
2365  }
2366 
2367  return;
2368 }
2369 //---------------------------------------------------
2370 template <class V, class M>
2371 void
2372 SequenceOfVectors<V,M>::fftForward(
2373  unsigned int initialPos,
2374  unsigned int fftSize,
2375  unsigned int paramId,
2376  std::vector<std::complex<double> >& fftResult) const
2377 {
2378  bool bRC = ((initialPos < this->subSequenceSize()) &&
2379  (paramId < this->vectorSizeLocal()) &&
2380  (0 < fftSize ) &&
2381  ((initialPos+fftSize) <= this->subSequenceSize()) &&
2382  (fftSize < this->subSequenceSize()));
2383  queso_require_msg(bRC, "invalid input data");
2384 
2385  std::vector<double> rawData(fftSize,0.);
2386  this->extractRawData(initialPos,
2387  1, // spacing
2388  fftSize,
2389  paramId,
2390  rawData);
2391 
2392  m_fftObj->forward(rawData,fftSize,fftResult);
2393 
2394  return;
2395 }
2396 //---------------------------------------------------
2397 template <class V, class M>
2398 void
2399 SequenceOfVectors<V,M>::psd(
2400  unsigned int initialPos,
2401  unsigned int numBlocks,
2402  double hopSizeRatio,
2403  unsigned int paramId,
2404  std::vector<double>& psdResult) const
2405 {
2406  bool bRC = ((initialPos < this->subSequenceSize()) &&
2407  (paramId < this->vectorSizeLocal()));
2408  queso_require_msg(bRC, "invalid input data");
2409 
2410  ScalarSequence<double> data(m_env,0,"");
2411 
2412  this->extractScalarSeq(initialPos,
2413  1, // spacing
2414  this->subSequenceSize()-initialPos,
2415  paramId,
2416  data);
2417  data.psd(0,
2418  numBlocks,
2419  hopSizeRatio,
2420  psdResult);
2421 
2422  return;
2423 }
2424 //---------------------------------------------------
2425 template <class V, class M>
2426 void
2427 SequenceOfVectors<V,M>::psdAtZero(
2428  unsigned int initialPos,
2429  unsigned int numBlocks,
2430  double hopSizeRatio,
2431  V& psdVec) const
2432 {
2433  bool bRC = ((initialPos < this->subSequenceSize()) &&
2434  (this->vectorSizeLocal() == psdVec.sizeLocal()));
2435  queso_require_msg(bRC, "invalid input data");
2436 
2437  ScalarSequence<double> data(m_env,0,"");
2438  std::vector<double> psdResult(0,0.); // size will be determined by 'data.psd()'
2439 
2440  unsigned int numParams = this->vectorSizeLocal();
2441  for (unsigned int i = 0; i < numParams; ++i) {
2442  this->extractScalarSeq(initialPos,
2443  1, // spacing
2444  this->subSequenceSize()-initialPos,
2445  i,
2446  data);
2447  data.psd(0,
2448  numBlocks,
2449  hopSizeRatio,
2450  psdResult);
2451  psdVec[i] = psdResult[0];
2452  //*m_env.subDisplayFile() << "psdResult[0] = " << psdResult[0] << std::endl;
2453  }
2454 
2455  return;
2456 }
2457 //---------------------------------------------------
2458 template <class V, class M>
2459 void
2460 SequenceOfVectors<V,M>::geweke(
2461  unsigned int initialPos,
2462  double ratioNa,
2463  double ratioNb,
2464  V& gewVec) const
2465 {
2466  bool bRC = ((initialPos < this->subSequenceSize()) &&
2467  (this->vectorSizeLocal() == gewVec.sizeLocal() ));
2468  queso_require_msg(bRC, "invalid input data");
2469 
2470  unsigned int numPos = this->subSequenceSize() - initialPos;
2471  ScalarSequence<double> data(m_env,0,"");
2472 
2473  unsigned int numParams = this->vectorSizeLocal();
2474  for (unsigned int i = 0; i < numParams; ++i) {
2475  this->extractScalarSeq(initialPos,
2476  1, // spacing
2477  numPos,
2478  i,
2479  data);
2480  gewVec[i] = data.geweke(0,
2481  ratioNa,
2482  ratioNb);
2483  }
2484 
2485  return;
2486 }
2487 //---------------------------------------------------
2488 template <class V, class M>
2489 void
2490 SequenceOfVectors<V,M>::meanStacc(
2491  unsigned int initialPos,
2492  V& meanStaccVec) const
2493 {
2494  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2495  (this->vectorSizeLocal() == meanStaccVec.sizeLocal()));
2496  queso_require_msg(bRC, "invalid input data");
2497 
2498  unsigned int numPos = this->subSequenceSize() - initialPos;
2499  ScalarSequence<double> data(m_env,0,"");
2500 
2501  unsigned int numParams = this->vectorSizeLocal();
2502  for (unsigned int i = 0; i < numParams; ++i) {
2503  this->extractScalarSeq(initialPos,
2504  1, // spacing
2505  numPos,
2506  i,
2507  data);
2508  meanStaccVec[i] = data.meanStacc(0);
2509  }
2510 
2511  return;
2512 }
2513 //---------------------------------------------------
2514 template <class V, class M>
2515 void
2516 SequenceOfVectors<V,M>::subCdfPercentageRange(
2517  unsigned int initialPos,
2518  unsigned int numPos,
2519  double range, // \in [0,1]
2520  V& lowerVec,
2521  V& upperVec) const
2522 {
2523  bool bRC = ((0 < numPos ) &&
2524  ((initialPos+numPos) <= this->subSequenceSize()) &&
2525  (this->vectorSizeLocal() == lowerVec.sizeLocal() ) &&
2526  (this->vectorSizeLocal() == upperVec.sizeLocal() ));
2527  queso_require_msg(bRC, "invalid input data");
2528 
2529  unsigned int numParams = this->vectorSizeLocal();
2530  ScalarSequence<double> data(m_env,0,"");
2531 
2532  for (unsigned int i = 0; i < numParams; ++i) {
2533  this->extractScalarSeq(initialPos,
2534  1, // spacing
2535  numPos,
2536  i,
2537  data);
2538  data.subCdfPercentageRange(0,
2539  numPos,
2540  range,
2541  lowerVec[i],
2542  upperVec[i]);
2543  }
2544 
2545  return;
2546 }
2547 //---------------------------------------------------
2548 template <class V, class M>
2549 void
2550 SequenceOfVectors<V,M>::unifiedCdfPercentageRange(
2551  unsigned int initialPos,
2552  unsigned int numPos,
2553  double range, // \in [0,1]
2554  V& lowerVec,
2555  V& upperVec) const
2556 {
2557  bool bRC = ((0 < numPos ) &&
2558  ((initialPos+numPos) <= this->subSequenceSize()) &&
2559  (this->vectorSizeLocal() == lowerVec.sizeLocal() ) &&
2560  (this->vectorSizeLocal() == upperVec.sizeLocal() ));
2561  queso_require_msg(bRC, "invalid input data");
2562 
2563  unsigned int numParams = this->vectorSizeLocal();
2564  ScalarSequence<double> data(m_env,0,"");
2565 
2566  for (unsigned int i = 0; i < numParams; ++i) {
2567  this->extractScalarSeq(initialPos,
2568  1, // spacing
2569  numPos,
2570  i,
2571  data);
2572  data.unifiedCdfPercentageRange(m_vectorSpace.numOfProcsForStorage() == 1,
2573  0,
2574  numPos,
2575  range,
2576  lowerVec[i],
2577  upperVec[i]);
2578  }
2579 
2580  return;
2581 }
2582 //---------------------------------------------------
2583 template <class V, class M>
2584 void
2585 SequenceOfVectors<V,M>::subCdfStacc(
2586  unsigned int initialPos,
2587  std::vector<V*>& cdfStaccVecs,
2588  std::vector<V*>& cdfStaccVecsUp,
2589  std::vector<V*>& cdfStaccVecsLow,
2590  std::vector<V*>& sortedDataVecs) const
2591 {
2592  bool bRC = (initialPos < this->subSequenceSize());
2593  queso_require_msg(bRC, "invalid input data");
2594 
2595  unsigned int numPos = this->subSequenceSize() - initialPos;
2596  unsigned int numEvals = numPos;
2597  for (unsigned int j = 0; j < numEvals; ++j) {
2598  cdfStaccVecs [j] = new V(m_vectorSpace.zeroVector());
2599  cdfStaccVecsUp [j] = new V(m_vectorSpace.zeroVector());
2600  cdfStaccVecsLow[j] = new V(m_vectorSpace.zeroVector());
2601  sortedDataVecs [j] = new V(m_vectorSpace.zeroVector());
2602  }
2603  std::vector<double> cdfStaccs (numEvals,0.);
2604  std::vector<double> cdfStaccsup (numEvals,0.);
2605  std::vector<double> cdfStaccslow(numEvals,0.);
2606 
2607  ScalarSequence<double> data (m_env,0,"");
2608  ScalarSequence<double> sortedData(m_env,0,"");
2609  unsigned int numParams = this->vectorSizeLocal();
2610  for (unsigned int i = 0; i < numParams; ++i) {
2611  this->extractScalarSeq(initialPos,
2612  1, // spacing
2613  numPos,
2614  i,
2615  data);
2616  //std::cout << "x-data" << data<< std::endl;
2617  data.subSort(initialPos,sortedData);
2618  data.subCdfStacc(initialPos,
2619  cdfStaccs,
2620  cdfStaccsup,
2621  cdfStaccslow,
2622  sortedData);
2623 
2624  for (unsigned int j = 0; j < numEvals; ++j) {
2625  (*sortedDataVecs [j])[i] = sortedData [j];
2626  (*cdfStaccVecs [j])[i] = cdfStaccs [j];
2627  (*cdfStaccVecsUp [j])[i] = cdfStaccsup [j];
2628  (*cdfStaccVecsLow[j])[i] = cdfStaccslow[j];
2629  }
2630  }
2631 
2632  return;
2633 }
2634 //---------------------------------------------------
2635 template <class V, class M>
2636 void
2637 SequenceOfVectors<V,M>::subCdfStacc(
2638  unsigned int initialPos,
2639  const std::vector<V*>& evalPositionsVecs,
2640  std::vector<V*>& cdfStaccVecs) const
2641 {
2642  bool bRC = ((initialPos < this->subSequenceSize() ) &&
2643  (0 < evalPositionsVecs.size()) &&
2644  (evalPositionsVecs.size() == cdfStaccVecs.size() ));
2645  queso_require_msg(bRC, "invalid input data");
2646 
2647  unsigned int numPos = this->subSequenceSize() - initialPos;
2648  ScalarSequence<double> data(m_env,0,"");
2649 
2650  unsigned int numEvals = evalPositionsVecs.size();
2651  for (unsigned int j = 0; j < numEvals; ++j) {
2652  cdfStaccVecs[j] = new V(m_vectorSpace.zeroVector());
2653  }
2654  std::vector<double> evalPositions(numEvals,0.);
2655  std::vector<double> cdfStaccs (numEvals,0.);
2656 
2657  unsigned int numParams = this->vectorSizeLocal();
2658  for (unsigned int i = 0; i < numParams; ++i) {
2659  this->extractScalarSeq(initialPos,
2660  1, // spacing
2661  numPos,
2662  i,
2663  data);
2664 
2665  for (unsigned int j = 0; j < numEvals; ++j) {
2666  evalPositions[j] = (*evalPositionsVecs[j])[i];
2667  }
2668 
2669  data.subCdfStacc(0,
2670  evalPositions,
2671  cdfStaccs);
2672 
2673  for (unsigned int j = 0; j < numEvals; ++j) {
2674  (*cdfStaccVecs[j])[i] = cdfStaccs[j];
2675  }
2676  }
2677 
2678  return;
2679 }
2680 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
2681 
2682 // --------------------------------------------------
2683 // --------------------------------------------------
2684 // --------------------------------------------------
2685 
2686 #ifdef UQ_CODE_HAS_MONITORS
2687 template <class V, class M>
2688 void
2689 SequenceOfVectors<V,M>::subMeanMonitorAlloc(unsigned int numberOfMonitorPositions)
2690 {
2691  m_subMeanMonitorPosSeq = new ScalarSequence<double>(m_env, numberOfMonitorPositions,(m_name+"_subMeanMonitorPosSeq").c_str());
2692  m_subMeanVecSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanVecSeq").c_str() );
2693  m_subMeanCltStdSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanCltStdSeq").c_str() );
2694 
2695  return;
2696 }
2697 // --------------------------------------------------
2698 template <class V, class M>
2699 void
2700 SequenceOfVectors<V,M>::subMeanInter0MonitorAlloc(unsigned int numberOfMonitorPositions)
2701 {
2702  m_subMeanInter0MonitorPosSeq = new ScalarSequence<double>(m_env, numberOfMonitorPositions,(m_name+"_subMeanInter0MonitorPosSeq").c_str() );
2703  m_subMeanInter0Mean = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0MeanSeq").c_str() );
2704  m_subMeanInter0Clt95 = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0Clt95Seq").c_str() );
2705  m_subMeanInter0Empirical90 = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0Empirical90Seq").c_str());
2706  m_subMeanInter0Min = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0MinSeq").c_str() );
2707  m_subMeanInter0Max = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_subMeanInter0MaxSeq").c_str() );
2708 
2709  return;
2710 }
2711 // --------------------------------------------------
2712 template <class V, class M>
2713 void
2714 SequenceOfVectors<V,M>::unifiedMeanMonitorAlloc(unsigned int numberOfMonitorPositions)
2715 {
2716  m_unifiedMeanMonitorPosSeq = new ScalarSequence<double>(m_env, numberOfMonitorPositions,(m_name+"_unifiedMeanMonitorPosSeq").c_str());
2717  m_unifiedMeanVecSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_unifiedMeanVecSeq").c_str() );
2718  m_unifiedMeanCltStdSeq = new SequenceOfVectors<V,M>(m_vectorSpace,numberOfMonitorPositions,(m_name+"_unifiedMeanCltStdSeq").c_str() );
2719 
2720  return;
2721 }
2722 // --------------------------------------------------
2723 template <class V, class M>
2724 void
2725 SequenceOfVectors<V,M>::subMeanMonitorRun(unsigned int monitorPosition,
2726  V& subMeanVec,
2727  V& subMeanCltStd)
2728 {
2729  this->subMeanExtra(0,
2730  monitorPosition,
2731  subMeanVec);
2732 
2733  this->subMeanCltStd(0,
2734  monitorPosition,
2735  subMeanVec,
2736  subMeanCltStd);
2737 
2738  return;
2739 }
2740 // --------------------------------------------------
2741 template <class V, class M>
2742 void
2743 SequenceOfVectors<V,M>::subMeanInter0MonitorRun(unsigned int monitorPosition,
2744  V& subMeanInter0Mean,
2745  V& subMeanInter0Clt95,
2746  V& subMeanInter0Empirical90,
2747  V& subMeanInter0Min,
2748  V& subMeanInter0Max)
2749 {
2750  V subMeanVec(m_vectorSpace.zeroVector());
2751  this->subMeanExtra(0,
2752  monitorPosition,
2753  subMeanVec);
2754 
2755  subMeanVec.mpiAllReduce(RawValue_MPI_SUM,m_env.inter0Comm(),subMeanInter0Mean);
2756  subMeanInter0Mean /= ((double) m_env.inter0Comm().NumProc());
2757 
2758  V subMeanInter0CltVariance = subMeanVec-subMeanInter0Mean;
2759  subMeanInter0CltVariance *= subMeanInter0CltVariance;
2760  subMeanInter0CltVariance.mpiAllReduce(RawValue_MPI_SUM,m_env.inter0Comm(),subMeanInter0Clt95);
2761  subMeanInter0Clt95 /= ((double) (m_env.inter0Comm().NumProc()-1));
2762  subMeanInter0Clt95 /= ((double) (m_env.inter0Comm().NumProc()-1));
2763  subMeanInter0Clt95.cwSqrt();
2764  subMeanInter0Clt95 *= 3.;
2765 
2766  V subMeanInter0Quantile5(m_vectorSpace.zeroVector());
2767  subMeanVec.mpiAllQuantile(.05,m_env.inter0Comm(),subMeanInter0Quantile5);
2768  V subMeanInter0Quantile95(m_vectorSpace.zeroVector());
2769  subMeanVec.mpiAllQuantile(.95,m_env.inter0Comm(),subMeanInter0Quantile95);
2770  subMeanInter0Empirical90 = subMeanInter0Quantile95 - subMeanInter0Quantile5;
2771 
2772  subMeanVec.mpiAllReduce(RawValue_MPI_MIN,m_env.inter0Comm(),subMeanInter0Min);
2773 
2774  subMeanVec.mpiAllReduce(RawValue_MPI_MAX,m_env.inter0Comm(),subMeanInter0Max);
2775 
2776  return;
2777 }
2778 // --------------------------------------------------
2779 template <class V, class M>
2780 void
2781 SequenceOfVectors<V,M>::unifiedMeanMonitorRun(unsigned int monitorPosition,
2782  V& unifiedMeanVec,
2783  V& unifiedMeanCltStd)
2784 {
2785  this->unifiedMeanExtra(0,
2786  monitorPosition,
2787  unifiedMeanVec);
2788 
2789  this->unifiedMeanCltStd(0,
2790  monitorPosition,
2791  unifiedMeanVec,
2792  unifiedMeanCltStd);
2793  return;
2794 }
2795 // --------------------------------------------------
2796 template <class V, class M>
2797 void
2798 SequenceOfVectors<V,M>::subMeanMonitorStore(unsigned int i,
2799  unsigned int monitorPosition,
2800  const V& subMeanVec,
2801  const V& subMeanCltStd)
2802 {
2803  (*m_subMeanMonitorPosSeq)[i] = monitorPosition;
2804  m_subMeanVecSeq->setPositionValues(i,subMeanVec);
2805  m_subMeanCltStdSeq->setPositionValues(i,subMeanCltStd);
2806 
2807  return;
2808 }
2809 // --------------------------------------------------
2810 template <class V, class M>
2811 void
2812 SequenceOfVectors<V,M>::subMeanInter0MonitorStore(unsigned int i,
2813  unsigned int monitorPosition,
2814  const V& subMeanInter0Mean,
2815  const V& subMeanInter0Clt95,
2816  const V& subMeanInter0Empirical90,
2817  const V& subMeanInter0Min,
2818  const V& subMeanInter0Max)
2819 {
2820  (*m_subMeanInter0MonitorPosSeq)[i] = monitorPosition;
2821  m_subMeanInter0Mean->setPositionValues(i,subMeanInter0Mean);
2822  m_subMeanInter0Clt95->setPositionValues(i,subMeanInter0Clt95);
2823  m_subMeanInter0Empirical90->setPositionValues(i,subMeanInter0Empirical90);
2824  m_subMeanInter0Min->setPositionValues(i,subMeanInter0Min);
2825  m_subMeanInter0Max->setPositionValues(i,subMeanInter0Max);
2826 
2827  return;
2828 }
2829 // --------------------------------------------------
2830 template <class V, class M>
2831 void
2832 SequenceOfVectors<V,M>::unifiedMeanMonitorStore(unsigned int i,
2833  unsigned int monitorPosition,
2834  V& unifiedMeanVec,
2835  V& unifiedMeanCltStd)
2836 {
2837  (*m_unifiedMeanMonitorPosSeq)[i] = monitorPosition;
2838  m_unifiedMeanVecSeq->setPositionValues(i,unifiedMeanVec);
2839  m_unifiedMeanCltStdSeq->setPositionValues(i,unifiedMeanCltStd);
2840 
2841  return;
2842 }
2843 // --------------------------------------------------
2844 template <class V, class M>
2845 void
2846 SequenceOfVectors<V,M>::subMeanMonitorWrite(std::ofstream& ofs)
2847 {
2848  m_subMeanMonitorPosSeq->subWriteContents(0,m_subMeanMonitorPosSeq->subSequenceSize(),ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2849  m_subMeanVecSeq->subWriteContents (0,m_subMeanVecSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2850  m_subMeanCltStdSeq->subWriteContents (0,m_subMeanCltStdSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2851 
2852  return;
2853 }
2854 // --------------------------------------------------
2855 template <class V, class M>
2856 void
2857 SequenceOfVectors<V,M>::subMeanInter0MonitorWrite(std::ofstream& ofs)
2858 {
2859  m_subMeanInter0MonitorPosSeq->subWriteContents(0,m_subMeanInter0MonitorPosSeq->subSequenceSize(),ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2860  m_subMeanInter0Mean->subWriteContents (0,m_subMeanInter0Mean->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2861  m_subMeanInter0Clt95->subWriteContents (0,m_subMeanInter0Clt95->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2862  m_subMeanInter0Empirical90->subWriteContents (0,m_subMeanInter0Empirical90->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2863  m_subMeanInter0Min->subWriteContents (0,m_subMeanInter0Min->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2864  m_subMeanInter0Max->subWriteContents (0,m_subMeanInter0Max->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m"
2865 
2866  return;
2867 }
2868 // --------------------------------------------------
2869 template <class V, class M>
2870 void
2871 SequenceOfVectors<V,M>::unifiedMeanMonitorWrite(std::ofstream& ofs)
2872 {
2873  // std::set<unsigned int> tmpSet;
2874  // tmpSet.insert(0);
2875  m_unifiedMeanMonitorPosSeq->subWriteContents(0,m_unifiedMeanMonitorPosSeq->subSequenceSize(),ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m" // Yes, 'subWriteContents()'
2876  m_unifiedMeanVecSeq->subWriteContents (0,m_unifiedMeanVecSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m" // Yes, 'subWriteContents()'
2877  m_unifiedMeanCltStdSeq->subWriteContents (0,m_unifiedMeanCltStdSeq->subSequenceSize(), ofs,UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT); // Yes, always ".m" // Yes, 'subWriteContents()'
2878 
2879  return;
2880 }
2881 // --------------------------------------------------
2882 template <class V, class M>
2883 void
2884 SequenceOfVectors<V,M>::subMeanMonitorFree()
2885 {
2886  delete m_subMeanMonitorPosSeq;
2887  m_subMeanMonitorPosSeq = NULL;
2888  delete m_subMeanVecSeq;
2889  m_subMeanVecSeq = NULL;
2890  delete m_subMeanCltStdSeq;
2891  m_subMeanCltStdSeq = NULL;
2892 
2893  return;
2894 }
2895 // --------------------------------------------------
2896 template <class V, class M>
2897 void
2898 SequenceOfVectors<V,M>::subMeanInter0MonitorFree()
2899 {
2900  delete m_subMeanInter0MonitorPosSeq;
2901  m_subMeanInter0MonitorPosSeq = NULL;
2902  delete m_subMeanInter0Mean;
2903  m_subMeanInter0Mean = NULL;
2904  delete m_subMeanInter0Clt95;
2905  m_subMeanInter0Clt95 = NULL;
2906  delete m_subMeanInter0Empirical90;
2907  m_subMeanInter0Empirical90 = NULL;
2908  delete m_subMeanInter0Min;
2909  m_subMeanInter0Min = NULL;
2910  delete m_subMeanInter0Max;
2911  m_subMeanInter0Max = NULL;
2912 
2913  return;
2914 }
2915 // --------------------------------------------------
2916 template <class V, class M>
2917 void
2918 SequenceOfVectors<V,M>::unifiedMeanMonitorFree()
2919 {
2920  delete m_unifiedMeanMonitorPosSeq;
2921  m_unifiedMeanMonitorPosSeq = NULL;
2922  delete m_unifiedMeanVecSeq;
2923  m_unifiedMeanVecSeq = NULL;
2924  delete m_unifiedMeanCltStdSeq;
2925  m_unifiedMeanCltStdSeq = NULL;
2926 
2927  return;
2928 }
2929 #endif // #ifdef UQ_CODE_HAS_MONITORS
2930 
2931 } // End namespace QUESO
2932 
Class to accommodate arrays of one-dimensional grid.
#define UQ_FILE_EXTENSION_FOR_MATLAB_FORMAT
Definition: Defines.h:102
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
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...
std::ofstream * ofsVar
Provides a stream interface to write data to files.
Definition: Environment.h:81
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.
#define UQ_FILE_EXTENSION_FOR_TXT_FORMAT
Definition: Defines.h:103
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
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...
#define queso_error_msg(msg)
Definition: asserts.h:47
void subGaussian1dKde(unsigned int initialPos, const V &scaleVec, const std::vector< V * > &evalParamVecs, std::vector< V * > &densityVecs) const
Gaussian kernel for the KDE estimate of the sub-sequence.
void unifiedWriteContents(const std::string &fileName, const std::string &fileType) const
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...
istream * dataIn
Definition: ann_sample.cpp:58
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
void subInterQuantileRange(unsigned int initialPos, V &iqrVec) const
Returns the interquartile range of the values in the sub-sequence.
T subMedianExtra(unsigned int initialPos, unsigned int numPos) const
Finds the median value of the sub-sequence, considering numPos positions starting at position initial...
#define RawValue_MPI_SUM
Definition: MpiComm.h:71
void deleteStoredVectors()
Deletes all the stored vectors.
void unifiedUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &unifiedCdfGrids, ArrayOfOneDTables< V, M > &unifiedCdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void unifiedPopulationVariance(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedPopVec) const
Finds the population variance of the unified sequence, considering numPos positions starting at posit...
void 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...
#define queso_require_less_msg(expr1, expr2, msg)
Definition: asserts.h:87
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...
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
T unifiedInterQuantileRange(bool useOnlyInter0Comm, unsigned int initialPos) const
Returns the interquartile range of the values in the unified sequence.
void unifiedGaussian1dKde(unsigned int initialPos, const V &unifiedScaleVec, const std::vector< V * > &unifiedEvalParamVecs, std::vector< V * > &unifiedDensityVecs) const
Gaussian kernel for the KDE estimate of the unified sequence.
void 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.
SequenceOfVectors< V, M > & operator=(const SequenceOfVectors< V, M > &rhs)
Copies values from rhs to this.
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
Definition: GslMatrix.C:2001
void writeTxtHeader(std::ofstream &ofs, double sequenceSize, double vectorSizeLocal) const
Helper function to write plain txt info for vectors.
void setPositionValues(unsigned int posId, const V &vec)
Set the values in vec at position posId of the sequence.
SequenceOfVectors(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void writeUnifiedMatlabHeader(std::ofstream &ofs, double sequenceSize, double vectorSizeLocal) const
double MiscGetEllapsedSeconds(struct timeval *timeval0)
#define queso_require_less_equal_msg(expr1, expr2, msg)
Definition: asserts.h:89
double estimateConvBrooksGelman(unsigned int initialPos, unsigned int numPos) const
Estimates convergence rate using Brooks &amp; Gelman method.
Class for handling vector samples (sequence of vectors).
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors.
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90
void unifiedMinMaxExtra(bool useOnlyInter0Comm, unsigned int initialPos, unsigned int numPos, T &unifiedMinValue, T &unifiedMaxValue) const
Finds the minimum and the maximum values of the unified sequence, considering numPos positions starti...
void unifiedGaussian1dKde(bool useOnlyInter0Comm, unsigned int initialPos, double unifiedScaleValue, const std::vector< T > &unifiedEvaluationPositions, std::vector< double > &unifiedDensityValues) const
Gaussian kernel for the KDE estimate of the unified sequence.
void unifiedScalesForKde(unsigned int initialPos, const V &unifiedIqrVec, unsigned int kdeDimension, V &unifiedScaleVec) const
Selects the scales (bandwidth) for the kernel density estimation, considering the unified sequence...
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...
const int UQ_OK_RC
Definition: Defines.h:89
#define RawValue_MPI_MIN
Definition: MpiComm.h:69
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...
unsigned int subSequenceSize() const
Size of the sub-sequence of vectors.
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...
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.
Struct for handling data input and output from files.
Definition: Environment.h:72
~SequenceOfVectors()
Destructor.
T autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag) const
Calculates the autocorrelation via definition.
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...
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...
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
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...
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
void select(const std::vector< unsigned int > &idsOfUniquePositions)
TODO: It shall select positions in the sequence of vectors.
void unifiedInterQuantileRange(unsigned int initialPos, V &unifiedIqrVec) const
Returns the interquartile range of the values in the unified sequence.
void resetValues(unsigned int initialPos, unsigned int numPos)
Resets a total of numPos values of the sequence starting at position initialPos.
#define UQ_FILE_EXTENSION_FOR_HDF_FORMAT
Definition: Defines.h:104
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 ...
#define RawValue_MPI_MAX
Definition: MpiComm.h:70
void unifiedHistogram(unsigned int initialPos, const V &unifiedMinVec, const V &unifiedMaxVec, std::vector< V * > &unifiedCentersForAllBins, std::vector< V * > &unifiedQuanttsForAllBins) const
Calculates the histogram of the unified sequence.
void subUniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &cdfGrids, ArrayOfOneDTables< V, M > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates the autocorrelation via definition.
Class to accommodate arrays of one-dimensional tables.
T subInterQuantileRange(unsigned int initialPos) const
Returns the interquartile range of the values in the sub-sequence.
void subWriteContents(unsigned int initialPos, unsigned int numPos, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Writes the sub-sequence to a file.
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...
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
std::ifstream * ifsVar
Provides a stream interface to read data from files.
Definition: Environment.h:84
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
std::vector< const V * >::iterator seqVectorPositionIteratorTypedef
void unifiedMedianExtra(unsigned int initialPos, unsigned int numPos, V &unifiedMedianVec) const
Finds the median value of the unfed sequence, considering numPos positions starting at position initi...
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
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 ...
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...
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 ...
void unifiedSampleStd(unsigned int initialPos, unsigned int numPos, const V &unifiedMeanVec, V &unifiedStdVec) const
Finds the sample standard deviation of the unified sequence, considering numPos positions starting at...
void copy(const SequenceOfVectors< V, M > &src)
Copies vector sequence src to this.
A class representing a vector space.
Definition: VectorSet.h:49
T autoCovariance(unsigned int initialPos, unsigned int numPos, const T &meanValue, unsigned int lag) const
Calculates the autocovariance.
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
void autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
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...
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 ...
void copy(const BaseVectorSequence< V, M > &src)
Copies vector sequence src to this.
void subMeanExtra(unsigned int initialPos, unsigned int numPos, V &meanVec) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, ScalarSequence< double > &scalarSeq) const
Extracts a sequence of scalars.
void subScalesForKde(unsigned int initialPos, const V &iqrVec, unsigned int kdeDimension, V &scaleVec) const
Selects the scales (bandwidth, scaleVec) for the kernel density estimation, considering only the sub-...
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, V &minVec, V &maxVec) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
void subGaussian1dKde(unsigned int initialPos, double scaleValue, const std::vector< T > &evaluationPositions, std::vector< double > &densityValues) const
Gaussian kernel for the KDE estimate of the sub-sequence.
void writeSubMatlabHeader(std::ofstream &ofs, double sequenceSize, double vectorSizeLocal) const
Helper function to write matlab-specific header info for vectors.
Base class for handling vector and array samples (sequence of vectors or arrays). ...
std::vector< const V * > m_seq
Sequence of vectors.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.

Generated on Fri Jun 17 2016 14:10:48 for queso-0.54.0 by  doxygen 1.8.5