queso-0.52.0
ArrayOfSequences.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/ArrayOfSequences.h>
26 #include <queso/VectorSequence.h>
27 
28 namespace QUESO {
29 
30 // Default constructor
31 template <class V, class M>
32 ArrayOfSequences<V,M>::ArrayOfSequences(const VectorSpace<V,M>& vectorSpace,
33  unsigned int subSequenceSize,
34  const std::string& name)
35  : BaseVectorSequence<V,M>(vectorSpace,subSequenceSize,name),
36  m_scalarSequences(m_vectorSpace.map(),1)
37 {
38 
39  //if (m_env.subDisplayFile()) {
40  // *m_env.subDisplayFile() << "Entering ArrayOfSequences<V,M>::constructor()"
41  // << std::endl;
42  //}
43 
44  //if (m_env.subDisplayFile()) {
45  // *m_env.subDisplayFile() << "In ArrayOfSequences<V,M>::constructor()"
46  // << "\n subSequenceSize = " << subSequenceSize
47  // << "\n m_scalarSequences.MyLength() = " << m_scalarSequences.MyLength()
48  // << std::endl;
49  //}
50 
51  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
52  m_scalarSequences(i,0) = new ScalarSequence<double>(m_env,subSequenceSize);
53  }
54 
55  //if (m_env.subDisplayFile()) {
56  // *m_env.subDisplayFile() << "Leaving ArrayOfSequences<V,M>::constructor()"
57  // << std::endl;
58  //}
59 }
60 
61 // Destructor
62 template <class V, class M>
63 ArrayOfSequences<V,M>::~ArrayOfSequences()
64 {
65  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
66  if (m_scalarSequences(i,0)) delete m_scalarSequences(i,0);
67  }
68 }
69 
70 // Math methods
71 template <class V, class M>
72 unsigned int ArrayOfSequences<V,M>::subSequenceSize() const
73 {
74  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
75 
76  return tmp->m_scalarSequences(0,0)->subSequenceSize();
77 }
78 
79 template <class V, class M>
80 void ArrayOfSequences<V,M>::resizeSequence(unsigned int newSubSequenceSize)
81 {
82  if (newSubSequenceSize != this->subSequenceSize()) {
83  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
84  m_scalarSequences(i,0)->resizeSequence(newSubSequenceSize);
85  }
86  }
87 
88  BaseVectorSequence<V,M>::deleteStoredVectors();
89  return;
90 }
91 
92 template <class V, class M>
93 void ArrayOfSequences<V,M>::resetValues(unsigned int initialPos,
94  unsigned int numPos)
95 {
96  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
97  m_scalarSequences(i,0)->resetValues(initialPos,numPos);
98  }
99 
100  BaseVectorSequence<V,M>::deleteStoredVectors();
101 
102  return;
103 }
104 
105 template <class V, class M>
106 void ArrayOfSequences<V,M>::erasePositions(unsigned int initialPos,
107  unsigned int numPos)
108 {
109  if (initialPos < this->subSequenceSize()) {
110  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
111  ScalarSequence<double>& seq = *(m_scalarSequences(i,0));
112  seq.erasePositions(initialPos,numPos);
113  }
114  }
115 
116  BaseVectorSequence<V,M>::deleteStoredVectors();
117 
118  return;
119 }
120 
121 template <class V, class M>
122 void ArrayOfSequences<V,M>::getPositionValues(unsigned int posId, V& vec) const
123 {
124  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
125  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
126  vec[i] = (*(tmp->m_scalarSequences(i,0)))[posId];
127  }
128 
129  return;
130 }
131 
132 template <class V, class M>
133 void ArrayOfSequences<V,M>::setPositionValues(unsigned int posId, const V& vec)
134 {
135  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
136  ScalarSequence<double>& seq = *(m_scalarSequences(i,0));
137  seq[posId] = vec[i];
138  }
139 
140  BaseVectorSequence<V,M>::deleteStoredVectors();
141 
142  return;
143 }
144 
145 template <class V, class M>
146 void ArrayOfSequences<V,M>::setGaussian(const V& meanVec, const V& stdDevVec)
147 {
148  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
149  ScalarSequence<double>& seq = *(m_scalarSequences(i,0));
150  seq.setGaussian(meanVec[i],stdDevVec[i]);
151  }
152 
153  BaseVectorSequence<V,M>::deleteStoredVectors();
154 
155  return;
156 }
157 
158 template <class V, class M>
159 void ArrayOfSequences<V,M>::setUniform(const V& aVec, const V& bVec)
160 {
161  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
162  ScalarSequence<double>& seq = *(m_scalarSequences(i,0));
163  seq.setUniform(aVec[i],bVec[i]);
164  }
165 
166  BaseVectorSequence<V,M>::deleteStoredVectors();
167 
168  return;
169 }
170 
171 template <class V, class M>
172 void ArrayOfSequences<V,M>::mean(unsigned int initialPos, unsigned int numPos,
173  V& meanVec) const
174 {
175  bool bRC = ((0 <= initialPos ) &&
176  (0 < numPos ) &&
177  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
178  UQ_FATAL_TEST_MACRO(bRC == false,
179  m_env.worldRank(),
180  "ArrayOfSequences<V,M>::mean()",
181  "invalid initial position or number of positions");
182 
183  bRC = (this->vectorSize() == meanVec.size());
184  UQ_FATAL_TEST_MACRO(bRC == false,
185  m_env.worldRank(),
186  "ArrayOfSequences<V,M>::mean()",
187  "incompatible sizes between meanVec vector and vectors in sequence");
188 
189  meanVec.cwSet(0.);
190  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
191  for (unsigned int i = 0; i < meanVec.size(); ++i) {
192  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
193  meanVec[i] = seq.subMeanExtra(initialPos, numPos);
194  }
195 
196  return;
197 }
198 
199 template <class V, class M>
200 void ArrayOfSequences<V,M>::sampleVariance(unsigned int initialPos,
201  unsigned int numPos,
202  const V& meanVec,
203  V& samVec) const
204 {
205  bool bRC = ((0 <= initialPos ) &&
206  (0 < numPos ) &&
207  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
208  UQ_FATAL_TEST_MACRO(bRC == false,
209  m_env.worldRank(),
210  "ArrayOfSequences<V,M>::sampleVariance()",
211  "invalid initial position or number of positions");
212 
213  bRC = (this->vectorSize() == samVec.size());
214  UQ_FATAL_TEST_MACRO(bRC == false,
215  m_env.worldRank(),
216  "ArrayOfSequences<V,M>::sampleVariance()",
217  "incompatible sizes between samVec vector and vectors in sequence");
218 
219  bRC = (this->vectorSize() == meanVec.size());
220  UQ_FATAL_TEST_MACRO(bRC == false,
221  m_env.worldRank(),
222  "ArrayOfSequences<V,M>::sampleVariance()",
223  "incompatible sizes between meanVec vector and vectors in sequence");
224 
225  unsigned int loopSize = numPos;
226  unsigned int finalPosPlus1 = initialPos + loopSize;
227  double doubleLoopSize = (double) loopSize;
228  samVec.cwSet(0.);
229 
230  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
231  for (unsigned int i = 0; i < samVec.size(); ++i) {
232  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
233  double tmpMean = meanVec[i];
234  double result = 0.;
235  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
236  double diff = seq[j] - tmpMean;
237  result += diff*diff;
238  }
239  samVec[i] = result/(doubleLoopSize - 1.);
240  }
241 
242  return;
243 }
244 
245 template <class V, class M>
246 void ArrayOfSequences<V,M>::populationVariance(unsigned int initialPos,
247  unsigned int numPos,
248  const V& meanVec,
249  V& popVec) const
250 {
251  bool bRC = ((0 <= initialPos ) &&
252  (0 < numPos ) &&
253  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
254  UQ_FATAL_TEST_MACRO(bRC == false,
255  m_env.worldRank(),
256  "ArrayOfSequences<V,M>::populationVariance()",
257  "invalid initial position or number of positions");
258 
259  bRC = (this->vectorSize() == popVec.size());
260  UQ_FATAL_TEST_MACRO(bRC == false,
261  m_env.worldRank(),
262  "ArrayOfSequences<V,M>::populationVariance()",
263  "incompatible sizes between popVec vector and vectors in sequence");
264 
265  bRC = (this->vectorSize() == meanVec.size());
266  UQ_FATAL_TEST_MACRO(bRC == false,
267  m_env.worldRank(),
268  "ArrayOfSequences<V,M>::populationVariance()",
269  "incompatible sizes between meanVec vector and vectors in sequence");
270 
271  unsigned int loopSize = numPos;
272  unsigned int finalPosPlus1 = initialPos + loopSize;
273  double doubleLoopSize = (double) loopSize;
274  popVec.cwSet(0.);
275 
276  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
277  for (unsigned int i = 0; i < popVec.size(); ++i) {
278  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
279  double tmpMean = meanVec[i];
280  double result = 0.;
281  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
282  double diff = seq[j] - tmpMean;
283  result += diff*diff;
284  }
285  popVec[i] = result/doubleLoopSize;
286  }
287 
288  return;
289 }
290 
291 template <class V, class M>
292 void ArrayOfSequences<V,M>::autoCovariance(unsigned int initialPos,
293  unsigned int numPos,
294  const V& meanVec,
295  unsigned int lag,
296  V& covVec) const
297 {
298  bool bRC = ((0 <= initialPos ) &&
299  (0 < numPos ) &&
300  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
301  UQ_FATAL_TEST_MACRO(bRC == false,
302  m_env.worldRank(),
303  "VectorSequenceAutoCovariance<V,M>()",
304  "invalid initial position or number of positions");
305 
306  bRC = (numPos > lag);
307  UQ_FATAL_TEST_MACRO(bRC == false,
308  m_env.worldRank(),
309  "VectorSequenceAutoCovariance<V,M>()",
310  "lag is too large");
311 
312  bRC = (this->vectorSize() == covVec.size());
313  UQ_FATAL_TEST_MACRO(bRC == false,
314  m_env.worldRank(),
315  "VectorSequenceAutoCovariance<V,M>()",
316  "incompatible sizes between covVec vector and vectors in sequence");
317 
318  bRC = (this->vectorSize() == meanVec.size());
319  UQ_FATAL_TEST_MACRO(bRC == false,
320  m_env.worldRank(),
321  "VectorSequenceAutoCovariance<V,M>()",
322  "incompatible sizes between meanVec vector and vectors in sequence");
323 
324  unsigned int loopSize = numPos - lag;
325  unsigned int finalPosPlus1 = initialPos + loopSize;
326  double doubleLoopSize = (double) loopSize;
327  covVec.cwSet(0.);
328 
329  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
330  for (unsigned int i = 0; i < covVec.size(); ++i) {
331  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
332  double meanValue = meanVec[i];
333  double result = 0;
334  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
335  double diff1 = seq[j] - meanValue;
336  double diff2 = seq[j+lag] - meanValue;
337  result += diff1*diff2;
338  }
339  covVec[i] = result/doubleLoopSize;
340  }
341 
342  return;
343 }
344 
345 template <class V, class M>
346 void ArrayOfSequences<V,M>::autoCorrViaDef(unsigned int initialPos,
347  unsigned int numPos,
348  unsigned int lag,
349  V& corrVec) const
350 {
351  V subChainMean (m_vectorSpace.zeroVector());
352  V subChainAutoCovarianceLag0(m_vectorSpace.zeroVector());
353 
354  this->mean(initialPos,
355  numPos,
356  subChainMean);
357  this->autoCovariance(initialPos,
358  numPos,
359  subChainMean,
360  0, // lag
361  subChainAutoCovarianceLag0);
362 
363  this->autoCovariance(initialPos,
364  numPos,
365  subChainMean,
366  lag,
367  corrVec);
368  corrVec /= subChainAutoCovarianceLag0;
369 
370  return;
371 }
372 
373 template <class V, class M>
374 void ArrayOfSequences<V,M>::autoCorrViaFft(unsigned int initialPos,
375  unsigned int numPos,
376  const std::vector<unsigned int>& lags,
377  std::vector<V*>& corrVecs) const
378 {
379  return;
380 }
381 
382 template <class V, class M>
383 void ArrayOfSequences<V,M>::autoCorrViaFft(unsigned int initialPos,
384  unsigned int numPos,
385  unsigned int numSum,
386  V& autoCorrsSumVec) const
387 {
388 #if 0
389  bool bRC = ((initialPos < this->subSequenceSize()) &&
390  (0 < numPos ) &&
391  ((initialPos+numPos) <= this->subSequenceSize()) &&
392  (autoCorrsSumVec.size() == this->vectorSize() ));
393  UQ_FATAL_TEST_MACRO(bRC == false,
394  m_env.worldRank(),
395  "ArrayOfSequences<V,M>::autoCorrViaFft(), for sum",
396  "invalid input data");
397 
398  ScalarSequence<double> data(m_env,0);
399 
400  unsigned int numParams = this->vectorSize();
401  for (unsigned int i = 0; i < numParams; ++i) {
402  this->extractScalarSeq(initialPos,
403  1, // spacing
404  numPos,
405  i,
406  data);
407 
408  data.autoCorrViaFft(0,
409  numPos,
410  autoCorrsSumVec[i]);
411  }
412 #endif
413  return;
414 }
415 
416 template <class V, class M>
417 void ArrayOfSequences<V,M>::minMax(unsigned int initialPos, V& minVec,
418  V& maxVec) const
419 {
420  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
421 
422  unsigned int numParams = this->vectorSize();
423  for (unsigned int i = 0; i < numParams; ++i) {
424  ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
425  seq.subMinMaxExtra(initialPos,minVec[i],maxVec[i]);
426  }
427 
428  return;
429 }
430 
431 template <class V, class M>
432 void ArrayOfSequences<V,M>::histogram(unsigned int initialPos,
433  const V& minVec,
434  const V& maxVec,
435  std::vector<V*>& centersForAllBins,
436  std::vector<V*>& quanttsForAllBins) const
437 {
438 #if 0
439  UQ_FATAL_TEST_MACRO(centersForAllBins.size() != quanttsForAllBins.size(),
440  sequence[0]->env().worldRank(),
441  "VectorSequenceHistogram<V,M>()",
442  "vectors 'centers' and 'quantts' have different sizes");
443 
444  for (unsigned int j = 0; j < quanttsForAllBins.size(); ++j) {
445  centersForAllBins[j] = new V(*(sequence[0]));
446  quanttsForAllBins [j] = new V(*(sequence[0]));
447  }
448 
449  unsigned int dataSize = sequence.size() - initialPos;
450  unsigned int numParams = sequence[0]->size();
451  for (unsigned int i = 0; i < numParams; ++i) {
452  ScalarSequence<double> data(dataSize,0.);
453  for (unsigned int j = 0; j < dataSize; ++j) {
454  data[j] = (*(sequence[initialPos+j]))[i];
455  }
456 
457  std::vector<double> centers(centersForAllBins.size(),0.);
458  std::vector<double> quantts(quanttsForAllBins.size(),0.);
459  data.histogram(minVec[i],
460  maxVec[i],
461  centers,
462  quantts);
463 
464  for (unsigned int j = 0; j < quantts.size(); ++j) {
465  (*(centersForAllBins[j]))[i] = centers[j];
466  (*(quanttsForAllBins[j]))[i] = quantts[j];
467  }
468  }
469 
470 #endif
471  return;
472 }
473 
474 template <class V, class M>
475 void ArrayOfSequences<V,M>::interQuantileRange(unsigned int initialPos,
476  V& iqrs) const
477 {
478 #if 0
479  unsigned int dataSize = sequence.size() - initialPos;
480 
481  ArrayOfSequences sortedSequence(dataSize,m_vectorSpace.zeroVector());
482  this->sort(initialPos,
483  sortedSequence);
484 
485  unsigned int pos1 = (unsigned int) ( (((double) dataSize) + 1.)*1./4. - 1. );
486  unsigned int pos3 = (unsigned int) ( (((double) dataSize) + 1.)*3./4. - 1. );
487 
488  double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
489  double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
490 
491  unsigned int numParams = sequence[0]->size();
492  for (unsigned int i = 0; i < numParams; ++i) {
493  double value1 = (1.-fraction1) * (*sortedSequence[pos1])[i] + fraction1 * (*sortedSequence[pos1+1])[i];
494  double value3 = (1.-fraction3) * (*sortedSequence[pos3])[i] + fraction3 * (*sortedSequence[pos3+1])[i];
495  iqrs[i] = value3 - value1;
496  }
497 
498 #endif
499  return;
500 }
501 
502 template <class V, class M>
503 void ArrayOfSequences<V,M>::scalesForKDE(unsigned int initialPos,
504  const V& iqrs,
505  unsigned int kdeDimension,
506  V& scales) const
507 {
508 #if 0
509  unsigned int dataSize = sequence.size() - initialPos;
510 
511  V mean(*(sequence[0]));
512  VectorSequenceMean(sequence,
513  initialPos,
514  dataSize,
515  mean);
516 
517  V samVec(*(sequence[0]));
518  VectorSequenceSampleVariance(sequence,
519  initialPos,
520  dataSize,
521  mean,
522  samVec);
523 
524  unsigned int numParams = sequence[0]->size();
525  for (unsigned int i = 0; i < numParams; ++i) {
526  if (iqrs[i] <= 0.) {
527  scales[i] = 1.06*std::sqrt(samVec[i])/std::pow(dataSize,1./5.);
528  }
529  else {
530  scales[i] = 1.06*std::min(std::sqrt(samVec[i]),iqrs[i]/1.34)/std::pow(dataSize,1./5.);
531  }
532  }
533 
534 #endif
535  return;
536 }
537 
538 template <class V, class M>
539 void ArrayOfSequences<V,M>::gaussianKDE(const V& evaluationParamVec,
540  V& densityVec) const
541 {
542  return;
543 }
544 
545 template <class V, class M>
546 void ArrayOfSequences<V,M>::gaussianKDE(unsigned int initialPos,
547  const V& scales,
548  const std::vector<V*>& evaluationParamVecs,
549  std::vector<V*>& densityVecs) const
550 {
551 #if 0
552  unsigned int dataSize = sequence.size() - initialPos;
553  unsigned int numEstimationsPerParam = evaluationParamVecs.size();
554 
555  for (unsigned int j = 0; j < numEstimationsPerParam; ++j) {
556  densityVecs[j] = new V(*(sequence[0]));
557  }
558 
559  unsigned int numParams = sequence[0]->size();
560  for (unsigned int i = 0; i < numParams; ++i) {
561  double scaleInv = 1./scales[i];
562  for (unsigned int j = 0; j < numEstimationsPerParam; ++j) {
563  double x = (*(evaluationParamVecs[j]))[i];
564  double value = 0.;
565  for (unsigned int k = 0; k < dataSize; ++k) {
566  double xk = (*(sequence[initialPos+k]))[i];
567  value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
568  }
569  (*(densityVecs[j]))[i] = scaleInv * (value/(double) numEstimationsPerParam);
570  }
571  }
572 
573 #endif
574  return;
575 }
576 
577 template <class V, class M>
578 void ArrayOfSequences<V,M>::writeContents(std::ofstream& ofsvar) const
579 {
580  // Write chain
581  ofsvar << m_name << "_sub" << m_env.subIdString() << " = zeros(" << this->subSequenceSize()
582  << "," << this->vectorSize()
583  << ");"
584  << std::endl;
585  ofsvar << m_name << "_sub" << m_env.subIdString() << " = [";
586 
587  V tmpVec(m_vectorSpace.zeroVector());
588  unsigned int chainSize = this->subSequenceSize();
589  for (unsigned int j = 0; j < chainSize; ++j) {
590  this->getPositionValues(j,tmpVec);
591  ofsvar << tmpVec
592  << std::endl;
593  }
594  ofsvar << "];\n";
595 
596  return;
597 }
598 
599 template <class V, class M>
600 void ArrayOfSequences<V,M>::unifiedWriteContents(std::ofstream& ofsvar) const
601 {
602  UQ_FATAL_TEST_MACRO(true,
603  m_env.worldRank(),
604  "ArrayOfSequences<V,M>::unifiedWriteContents(1)",
605  "not implemented yet");
606  return;
607 }
608 
609 template <class V, class M>
610 void ArrayOfSequences<V,M>::unifiedWriteContents(const std::string& fileName,
611  const std::string& fileType) const
612 {
613  UQ_FATAL_TEST_MACRO(true,
614  m_env.worldRank(),
615  "ArrayOfSequences<V,M>::unifiedWriteContents(2)",
616  "not implemented yet");
617  return;
618 }
619 
620 template <class V, class M>
621 void ArrayOfSequences<V,M>::unifiedReadContents(const std::string& fileName,
622  const std::string& fileType,
623  const unsigned int subSequenceSize)
624 {
625  UQ_FATAL_TEST_MACRO(true,
626  m_env.worldRank(),
627  "ArrayOfSequences<V,M>::unifiedReadContents()",
628  "not implemented yet");
629  return;
630 }
631 
632 template <class V, class M>
633 void ArrayOfSequences<V,M>::select(
634  const std::vector<unsigned int>& idsOfUniquePositions)
635 {
636  return;
637 }
638 
639 template <class V, class M>
640 void ArrayOfSequences<V,M>::filter(unsigned int initialPos,
641  unsigned int spacing)
642 {
643  return;
644 }
645 
646 template <class V, class M>
647 void ArrayOfSequences<V,M>::extractScalarSeq(unsigned int initialPos,
648  unsigned int spacing,
649  unsigned int numPos,
650  unsigned int paramId,
651  ScalarSequence<double>& scalarSeq) const
652 {
653  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
654  ScalarSequence<double>& seq = *(tmp->m_scalarSequences(paramId,0));
655 
656  scalarSeq.resizeSequence(numPos);
657  if (spacing == 1) {
658  for (unsigned int j = 0; j < numPos; ++j) {
659  scalarSeq[j] = seq[paramId];
660  }
661  }
662  else {
663  for (unsigned int j = 0; j < numPos; ++j) {
664  scalarSeq[j] = seq[paramId];
665  }
666  }
667 
668  return;
669 }
670 
671 // Private method
672 template <class V, class M>
673 void ArrayOfSequences<V,M>::extractRawData(unsigned int initialPos,
674  unsigned int spacing,
675  unsigned int numPos,
676  unsigned int paramId,
677  std::vector<double>& rawData) const
678 {
679  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
680  ScalarSequence<double>& seq = *(tmp->m_scalarSequences(paramId,0));
681 
682  rawData.resize(numPos);
683  if (spacing == 1) {
684  for (unsigned int j = 0; j < numPos; ++j) {
685  rawData[j] = seq[paramId];
686  }
687  }
688  else {
689  for (unsigned int j = 0; j < numPos; ++j) {
690  rawData[j] = seq[paramId];
691  }
692  }
693 
694  return;
695 }
696 
697 // Methods conditionally available ------------------
698 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
699 template <class V, class M>
700 void ArrayOfSequences<V,M>::uniformlySampledCdf(
701  const V& numEvaluationPointsVec,
702  ArrayOfOneDGrids <V,M>& cdfGrids,
703  ArrayOfOneDTables<V,M>& cdfValues) const
704 {
705  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
706  V minCdfValues(m_vectorSpace.zeroVector());
707  V maxCdfValues(m_vectorSpace.zeroVector());
708  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
709  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
710 
711  unsigned int numEvaluationPoints = (unsigned int) numEvaluationPointsVec[i];
712  std::vector<double> aCdf(0);
713  seq.subUniformlySampledCdf(numEvaluationPoints,
714  minCdfValues[i],
715  maxCdfValues[i],
716  aCdf);
717  cdfValues.setOneDTable(i,aCdf);
718  }
719  cdfGrids.setUniformGrids(numEvaluationPointsVec,
720  minCdfValues,
721  maxCdfValues);
722 
723  return;
724 }
725 #endif
726 
727 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
728 template <class V, class M>
729 void ArrayOfSequences<V,M>::bmm(unsigned int initialPos,
730  unsigned int batchLength,
731  V& bmmVec) const
732 {
733 #if 0
734  V meanOfBatchMeans (*(sequence[0]));
735  V covLag0OfBatchMeans(*(sequence[0]));
736  V covLag1OfBatchMeans(*(sequence[0]));
737 
738  V tmpVector(m_vectorSpace.zeroVector()); // In order to contour the fact that 'batchMeans' is a vector of 'const V*', but needs to be set first
739  for (unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
740  for (unsigned int batchLengthId = 0; batchLengthId < batchLengths.size(); batchLengthId++) {
741  unsigned int batchLength = batchLengths[batchLengthId];
742  unsigned int numberOfBatches = (sequence.size() - initialPositions[initialPosId])/batchLength;
743 
744  std::vector<const V* > batchMeans(numberOfBatches,NULL);
745  for (unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
746  VectorSequenceMean(sequence,
747  initialPositions[initialPosId] + batchId*batchLength,
748  batchLength,
749  tmpVector);
750  batchMeans[batchId] = new V(tmpVector);
751  }
752 
753  VectorSequenceMean(batchMeans,
754  0,
755  batchMeans.size(),
756  meanOfBatchMeans);
757 
758  VectorSequenceAutoCovariance(batchMeans,
759  0,
760  batchMeans.size(),
761  meanOfBatchMeans,
762  0, // lag
763  covLag0OfBatchMeans);
764 
765  VectorSequenceAutoCovariance(batchMeans,
766  0,
767  batchMeans.size(),
768  meanOfBatchMeans,
769  1, // lag
770  covLag0OfBatchMeans);
771 
772  VectorSequenceSampleVariance(batchMeans,
773  0,
774  batchMeans.size(),
775  meanOfBatchMeans,
776  _2dArrayOfBMM(initialPosId,batchLengthId));
777  //_2dArrayOfBMM(initialPosId,batchLengthId) /= (double) batchMeans.size(); // CHECK
778  _2dArrayOfBMM(initialPosId,batchLengthId) *= (double) (sequence.size() - initialPositions[initialPosId]); // CHECK
779 
780  for (unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
781  if (batchMeans[batchId] != NULL) delete batchMeans[batchId];
782  }
783  }
784  }
785 #endif
786  return;
787 }
788 
789 template <class V, class M>
790 void ArrayOfSequences<V,M>::fftForward(unsigned int initialPos,
791  unsigned int fftSize,
792  unsigned int paramId,
793  std::vector<std::complex<double> >& resultData) const
794 {
795  return;
796 }
797 
798 template <class V, class M>
799 void ArrayOfSequences<V,M>::psd(unsigned int initialPos,
800  unsigned int numBlocks,
801  double hopSizeRatio,
802  unsigned int paramId,
803  std::vector<double>& psdResult) const
804 {
805  return;
806 }
807 
808 template <class V, class M>
809 void ArrayOfSequences<V,M>::psdAtZero(unsigned int initialPos,
810  unsigned int numBlocks,
811  double hopSizeRatio,
812  V& psdVec) const
813 {
814 #if 0
815  for (unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
816  unsigned int dataSize = sequence.size() - initialPositions[initialPosId];
817  ScalarSequence<double> data(dataSize,0.);
818 
819  unsigned int numParams = sequence[0]->size();
820  for (unsigned int i = 0; i < numParams; ++i) {
821  for (unsigned int j = 0; j < dataSize; ++j) {
822  data[j] = (*(sequence[initialPositions[initialPosId]+j]))[i];
823  }
824  for (unsigned int numsOfBlocksId = 0; numsOfBlocksId < numsOfBlocks.size(); numsOfBlocksId++) {
825  unsigned int numBlocks = numsOfBlocks[numsOfBlocksId];
826  std::vector<double> psdSequence(0,0.); // size will be determined by 'ScalarSequencePSD()'
827  data.psd(numBlocks,
828  hopSizeRatio,
829  psdSequence);
830  _2dArrayOfPSDAtZero(initialPosId,numsOfBlocksId)[i] = psdSequence[0];
831  //if (m_env.subDisplayFile()) {
832  // *m_env.subDisplayFile() << "psdSequence[0] = " << psdSequence[0] << std::endl;
833  //}
834  } // for 'numsOfBlocksId'
835  } // for 'i'
836  }
837 #endif
838  return;
839 }
840 
841 template <class V, class M>
842 void ArrayOfSequences<V,M>::geweke(unsigned int initialPos, double ratioNa,
843  double ratioNb,
844  V& gewVec) const
845 {
846 #if 0
847  for (unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
848  unsigned int fullDataSize = sequence.size() - initialPositions[initialPosId];
849  unsigned int dataSizeA = (unsigned int) (((double) fullDataSize) * ratioNa);
850  unsigned int dataSizeB = (unsigned int) (((double) fullDataSize) * ratioNb);
851  unsigned int initialPosA = initialPositions[initialPosId];
852  unsigned int initialPosB = sequence.size() - dataSizeB;
853 
854  V meanA(*(sequence[0]));
855  VectorSequenceMean(sequence,
856  initialPosA,
857  dataSizeA,
858  meanA);
859 
860  V meanB(*(sequence[0]));
861  VectorSequenceMean(sequence,
862  initialPosB,
863  dataSizeB,
864  meanB);
865 
866  unsigned int numParams = sequence[0]->size();
867 
868  V psdVecA(*(sequence[0]));
869  ScalarSequence<double> dataA(dataSizeA,0.);
870  for (unsigned int i = 0; i < numParams; ++i) {
871  for (unsigned int j = 0; j < dataSizeA; ++j) {
872  dataA[j] = (*(sequence[initialPosA+j]))[i];
873  }
874  std::vector<double> psdSequence(0,0.);
875  dataA.psd(8, // numBlocks
876  .5, // hopSizeRatio
877  psdSequence);
878  psdVecA[i] = psdSequence[0];
879  } // for 'i'
880 
881  V psdVecB(*(sequence[0]));
882  ScalarSequence<double> dataB(dataSizeB,0.);
883  for (unsigned int i = 0; i < numParams; ++i) {
884  for (unsigned int j = 0; j < dataSizeB; ++j) {
885  dataB[j] = (*(sequence[initialPosB+j]))[i];
886  }
887  std::vector<double> psdSequence(0,0.);
888  dataB.psd(8, // numBlocks
889  .5, // hopSizeRatio
890  psdSequence);
891  psdVecB[i] = psdSequence[0];
892  } // for 'i'
893 
894  vectorOfGeweke[initialPosId] = new V(*(sequence[0]));
895 
896  double doubleDataSizeA = (double) dataSizeA;
897  double doubleDataSizeB = (double) dataSizeB;
898  for (unsigned int i = 0; i < numParams; ++i) {
899  (*(vectorOfGeweke[initialPosId]))[i] = (meanA[i] - meanB[i])/std::sqrt(psdVecA[i]/doubleDataSizeA + psdVecB[i]/doubleDataSizeB);
900  }
901  }
902 
903 #endif
904  return;
905 }
906 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
907 
908 } // End namespace QUESO
QUESO::ScalarSequence::subUniformlySampledCdf
void subUniformlySampledCdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
Definition: ScalarSequence.C:592
QUESO::BaseVectorSequence::deleteStoredVectors
void deleteStoredVectors()
Deletes all the stored vectors.
Definition: VectorSequence.C:301
QUESO::ArrayOfOneDGrids
Class to accommodate arrays of one-dimensional grid.
Definition: ArrayOfOneDGrids.h:44
QUESO::ScalarSequence::setUniform
void setUniform(const T &a, const T &b)
Sets the values of the sequence as a uniform distribution between the values given by vectors aVec an...
Definition: ScalarSequence.C:557
QUESO::ScalarSequence::erasePositions
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values of the sequence, starting at position initialPos.
Definition: ScalarSequence.C:217
QUESO::ArrayOfSequences::interQuantileRange
void interQuantileRange(unsigned int initialPos, V &iqrs) const
Returns the interquartile range of the values in the sequence.
Definition: ArrayOfSequences.C:475
QUESO::ArrayOfSequences::autoCorrViaFft
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
Definition: ArrayOfSequences.C:374
QUESO::ArrayOfSequences::gaussianKDE
void gaussianKDE(const V &evaluationParamVec, V &densityVec) const
Gaussian kernel for the KDE estimate of the sequence.
Definition: ArrayOfSequences.C:539
QUESO::ArrayOfSequences::resizeSequence
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
Definition: ArrayOfSequences.C:80
QUESO::ArrayOfSequences::erasePositions
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
Definition: ArrayOfSequences.C:106
QUESO::ArrayOfSequences::sampleVariance
void sampleVariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, V &samVec) const
Finds the sample variance of the sub-sequence, considering numPos positions starting at position init...
Definition: ArrayOfSequences.C:200
QUESO::ArrayOfOneDGrids::setUniformGrids
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
Definition: ArrayOfOneDGrids.C:111
QUESO::ArrayOfSequences::extractRawData
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
Definition: ArrayOfSequences.C:673
QUESO::ScalarSequence::subMinMaxExtra
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
Definition: ScalarSequence.C:1587
QUESO::ArrayOfSequences::resetValues
void resetValues(unsigned int initialPos, unsigned int numPos)
Resets a total of numPos values of the sequence starting at position initialPos.
Definition: ArrayOfSequences.C:93
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:181
UQ_FATAL_TEST_MACRO
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
Definition: Defines.h:223
QUESO::BaseVectorSequence::m_env
const BaseEnvironment & m_env
Definition: VectorSequence.h:423
QUESO::ArrayOfSequences::select
void select(const std::vector< unsigned int > &idsOfUniquePositions)
Select positions in the sequence of vectors.
Definition: ArrayOfSequences.C:633
QUESO::ArrayOfSequences::populationVariance
void populationVariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, V &popVec) const
Finds the population variance of the sub-sequence, considering numPos positions starting at position ...
Definition: ArrayOfSequences.C:246
QUESO::ArrayOfSequences::m_scalarSequences
DistArray< ScalarSequence< double > * > m_scalarSequences
Sequence of scalars.
Definition: ArrayOfSequences.h:261
QUESO::ArrayOfSequences
Class for handling array samples (arrays of scalar sequences).
Definition: ArrayOfSequences.h:44
QUESO::ArrayOfSequences::extractScalarSeq
void extractScalarSeq(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, ScalarSequence< double > &scalarSeq) const
Extracts a sequence of scalars of size numPos, from position paramId of the array of sequences...
Definition: ArrayOfSequences.C:647
QUESO::ArrayOfSequences::setPositionValues
void setPositionValues(unsigned int posId, const V &vec)
Set the values of vec at position posId of the sequence.
Definition: ArrayOfSequences.C:133
QUESO::ArrayOfOneDTables::setOneDTable
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
Definition: ArrayOfOneDTables.C:59
QUESO::ArrayOfSequences::unifiedReadContents
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
Definition: ArrayOfSequences.C:621
QUESO::ArrayOfSequences::scalesForKDE
void scalesForKDE(unsigned int initialPos, const V &iqrs, unsigned int kdeDimension, V &scales) const
Selects the scales (bandwidth, scaleVec) for the kernel density estimation, of the sequence...
Definition: ArrayOfSequences.C:503
QUESO::ArrayOfSequences::minMax
void minMax(unsigned int initialPos, V &minVec, V &maxVec) const
Given an initial position initialPos, finds the minimum and the maximum values of the sequence...
Definition: ArrayOfSequences.C:417
QUESO::ArrayOfSequences::autoCorrViaDef
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates autocorrelation via definition.
Definition: ArrayOfSequences.C:346
QUESO::ArrayOfSequences::mean
void mean(unsigned int initialPos, unsigned int numPos, V &meanVec) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
Definition: ArrayOfSequences.C:172
QUESO::ScalarSequence::subMeanExtra
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
Definition: ScalarSequence.C:862
QUESO::ArrayOfSequences::autoCovariance
void autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
Definition: ArrayOfSequences.C:292
QUESO::ArrayOfOneDTables
Class to accommodate arrays of one-dimensional tables.
Definition: ArrayOfOneDTables.h:46
QUESO::ArrayOfSequences::filter
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors, starting at initialPos, and with spacing given by spaci...
Definition: ArrayOfSequences.C:640
QUESO::ArrayOfSequences::ArrayOfSequences
ArrayOfSequences(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
Definition: ArrayOfSequences.C:32
QUESO::ArrayOfSequences::unifiedWriteContents
void unifiedWriteContents(std::ofstream &ofsvar) const
Writes info of the unified sequence to a file.
Definition: ArrayOfSequences.C:600
QUESO::MiscGaussianDensity
double MiscGaussianDensity(double x, double mu, double sigma)
Definition: Miscellaneous.C:311
QUESO::BaseVectorSequence
Base class for handling vector and array samples (sequence of vectors or arrays). ...
Definition: VectorSequence.h:54
QUESO::ArrayOfSequences::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of arrays.
Definition: ArrayOfSequences.C:72
QUESO::ArrayOfSequences::writeContents
void writeContents(std::ofstream &ofsvar) const
Write contents of the chain to a file.
Definition: ArrayOfSequences.C:578
QUESO::VectorSpace
A class representing a vector space.
Definition: VectorSet.h:46
QUESO::ArrayOfSequences::histogram
void histogram(unsigned int initialPos, const V &minVec, const V &maxVec, std::vector< V * > &centersForAllBins, std::vector< V * > &quanttsForAllBins) const
Calculates the histogram of the sequence.
Definition: ArrayOfSequences.C:432
QUESO::ArrayOfSequences::setUniform
void setUniform(const V &aVec, const V &bVec)
Sets the values of the sequence as a uniform distribution between the values given by vectors aVec an...
Definition: ArrayOfSequences.C:159
QUESO::ArrayOfSequences::setGaussian
void setGaussian(const V &meanVec, const V &stdDevVec)
Sets the values of the sequence as a Gaussian distribution of mean given by meanVec and standard devi...
Definition: ArrayOfSequences.C:146
QUESO::ArrayOfSequences::getPositionValues
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
Definition: ArrayOfSequences.C:122
QUESO::ArrayOfSequences::~ArrayOfSequences
~ArrayOfSequences()
Destructor.
Definition: ArrayOfSequences.C:63
QUESO::ScalarSequence::setGaussian
void setGaussian(const T &mean, const T &stdDev)
Sets the values of the sequence as a Gaussian distribution of mean given by meanVec and standard devi...
Definition: ScalarSequence.C:536
Generated on Thu Apr 23 2015 19:30:54 for queso-0.52.0 by  doxygen 1.8.5