queso-0.53.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  queso_require_msg(bRC, "invalid initial position or number of positions");
179 
180  bRC = (this->vectorSize() == meanVec.size());
181  queso_require_msg(bRC, "incompatible sizes between meanVec vector and vectors in sequence");
182 
183  meanVec.cwSet(0.);
184  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
185  for (unsigned int i = 0; i < meanVec.size(); ++i) {
186  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
187  meanVec[i] = seq.subMeanExtra(initialPos, numPos);
188  }
189 
190  return;
191 }
192 
193 template <class V, class M>
194 void ArrayOfSequences<V,M>::sampleVariance(unsigned int initialPos,
195  unsigned int numPos,
196  const V& meanVec,
197  V& samVec) const
198 {
199  bool bRC = ((0 <= initialPos ) &&
200  (0 < numPos ) &&
201  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
202  queso_require_msg(bRC, "invalid initial position or number of positions");
203 
204  bRC = (this->vectorSize() == samVec.size());
205  queso_require_msg(bRC, "incompatible sizes between samVec vector and vectors in sequence");
206 
207  bRC = (this->vectorSize() == meanVec.size());
208  queso_require_msg(bRC, "incompatible sizes between meanVec vector and vectors in sequence");
209 
210  unsigned int loopSize = numPos;
211  unsigned int finalPosPlus1 = initialPos + loopSize;
212  double doubleLoopSize = (double) loopSize;
213  samVec.cwSet(0.);
214 
215  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
216  for (unsigned int i = 0; i < samVec.size(); ++i) {
217  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
218  double tmpMean = meanVec[i];
219  double result = 0.;
220  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
221  double diff = seq[j] - tmpMean;
222  result += diff*diff;
223  }
224  samVec[i] = result/(doubleLoopSize - 1.);
225  }
226 
227  return;
228 }
229 
230 template <class V, class M>
231 void ArrayOfSequences<V,M>::populationVariance(unsigned int initialPos,
232  unsigned int numPos,
233  const V& meanVec,
234  V& popVec) const
235 {
236  bool bRC = ((0 <= initialPos ) &&
237  (0 < numPos ) &&
238  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
239  queso_require_msg(bRC, "invalid initial position or number of positions");
240 
241  bRC = (this->vectorSize() == popVec.size());
242  queso_require_msg(bRC, "incompatible sizes between popVec vector and vectors in sequence");
243 
244  bRC = (this->vectorSize() == meanVec.size());
245  queso_require_msg(bRC, "incompatible sizes between meanVec vector and vectors in sequence");
246 
247  unsigned int loopSize = numPos;
248  unsigned int finalPosPlus1 = initialPos + loopSize;
249  double doubleLoopSize = (double) loopSize;
250  popVec.cwSet(0.);
251 
252  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
253  for (unsigned int i = 0; i < popVec.size(); ++i) {
254  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
255  double tmpMean = meanVec[i];
256  double result = 0.;
257  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
258  double diff = seq[j] - tmpMean;
259  result += diff*diff;
260  }
261  popVec[i] = result/doubleLoopSize;
262  }
263 
264  return;
265 }
266 
267 template <class V, class M>
268 void ArrayOfSequences<V,M>::autoCovariance(unsigned int initialPos,
269  unsigned int numPos,
270  const V& meanVec,
271  unsigned int lag,
272  V& covVec) const
273 {
274  bool bRC = ((0 <= initialPos ) &&
275  (0 < numPos ) &&
276  ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
277  queso_require_msg(bRC, "invalid initial position or number of positions");
278 
279  bRC = (numPos > lag);
280  queso_require_msg(bRC, "lag is too large");
281 
282  bRC = (this->vectorSize() == covVec.size());
283  queso_require_msg(bRC, "incompatible sizes between covVec vector and vectors in sequence");
284 
285  bRC = (this->vectorSize() == meanVec.size());
286  queso_require_msg(bRC, "incompatible sizes between meanVec vector and vectors in sequence");
287 
288  unsigned int loopSize = numPos - lag;
289  unsigned int finalPosPlus1 = initialPos + loopSize;
290  double doubleLoopSize = (double) loopSize;
291  covVec.cwSet(0.);
292 
293  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
294  for (unsigned int i = 0; i < covVec.size(); ++i) {
295  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
296  double meanValue = meanVec[i];
297  double result = 0;
298  for (unsigned int j = initialPos; j < finalPosPlus1; ++j) {
299  double diff1 = seq[j] - meanValue;
300  double diff2 = seq[j+lag] - meanValue;
301  result += diff1*diff2;
302  }
303  covVec[i] = result/doubleLoopSize;
304  }
305 
306  return;
307 }
308 
309 template <class V, class M>
310 void ArrayOfSequences<V,M>::autoCorrViaDef(unsigned int initialPos,
311  unsigned int numPos,
312  unsigned int lag,
313  V& corrVec) const
314 {
315  V subChainMean (m_vectorSpace.zeroVector());
316  V subChainAutoCovarianceLag0(m_vectorSpace.zeroVector());
317 
318  this->mean(initialPos,
319  numPos,
320  subChainMean);
321  this->autoCovariance(initialPos,
322  numPos,
323  subChainMean,
324  0, // lag
325  subChainAutoCovarianceLag0);
326 
327  this->autoCovariance(initialPos,
328  numPos,
329  subChainMean,
330  lag,
331  corrVec);
332  corrVec /= subChainAutoCovarianceLag0;
333 
334  return;
335 }
336 
337 template <class V, class M>
338 void ArrayOfSequences<V,M>::autoCorrViaFft(unsigned int initialPos,
339  unsigned int numPos,
340  const std::vector<unsigned int>& lags,
341  std::vector<V*>& corrVecs) const
342 {
343  return;
344 }
345 
346 template <class V, class M>
347 void ArrayOfSequences<V,M>::autoCorrViaFft(unsigned int initialPos,
348  unsigned int numPos,
349  unsigned int numSum,
350  V& autoCorrsSumVec) const
351 {
352 #if 0
353  bool bRC = ((initialPos < this->subSequenceSize()) &&
354  (0 < numPos ) &&
355  ((initialPos+numPos) <= this->subSequenceSize()) &&
356  (autoCorrsSumVec.size() == this->vectorSize() ));
357  queso_require_msg(bRC, "invalid input data");
358 
359  ScalarSequence<double> data(m_env,0);
360 
361  unsigned int numParams = this->vectorSize();
362  for (unsigned int i = 0; i < numParams; ++i) {
363  this->extractScalarSeq(initialPos,
364  1, // spacing
365  numPos,
366  i,
367  data);
368 
369  data.autoCorrViaFft(0,
370  numPos,
371  autoCorrsSumVec[i]);
372  }
373 #endif
374  return;
375 }
376 
377 template <class V, class M>
378 void ArrayOfSequences<V,M>::minMax(unsigned int initialPos, V& minVec,
379  V& maxVec) const
380 {
381  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
382 
383  unsigned int numParams = this->vectorSize();
384  for (unsigned int i = 0; i < numParams; ++i) {
385  ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
386  seq.subMinMaxExtra(initialPos,minVec[i],maxVec[i]);
387  }
388 
389  return;
390 }
391 
392 template <class V, class M>
393 void ArrayOfSequences<V,M>::histogram(unsigned int initialPos,
394  const V& minVec,
395  const V& maxVec,
396  std::vector<V*>& centersForAllBins,
397  std::vector<V*>& quanttsForAllBins) const
398 {
399 #if 0
400  queso_require_equal_to_msg(centersForAllBins.size(), quanttsForAllBins.size(), "vectors 'centers' and 'quantts' have different sizes");
401 
402  for (unsigned int j = 0; j < quanttsForAllBins.size(); ++j) {
403  centersForAllBins[j] = new V(*(sequence[0]));
404  quanttsForAllBins [j] = new V(*(sequence[0]));
405  }
406 
407  unsigned int dataSize = sequence.size() - initialPos;
408  unsigned int numParams = sequence[0]->size();
409  for (unsigned int i = 0; i < numParams; ++i) {
410  ScalarSequence<double> data(dataSize,0.);
411  for (unsigned int j = 0; j < dataSize; ++j) {
412  data[j] = (*(sequence[initialPos+j]))[i];
413  }
414 
415  std::vector<double> centers(centersForAllBins.size(),0.);
416  std::vector<double> quantts(quanttsForAllBins.size(),0.);
417  data.histogram(minVec[i],
418  maxVec[i],
419  centers,
420  quantts);
421 
422  for (unsigned int j = 0; j < quantts.size(); ++j) {
423  (*(centersForAllBins[j]))[i] = centers[j];
424  (*(quanttsForAllBins[j]))[i] = quantts[j];
425  }
426  }
427 
428 #endif
429  return;
430 }
431 
432 template <class V, class M>
433 void ArrayOfSequences<V,M>::interQuantileRange(unsigned int initialPos,
434  V& iqrs) const
435 {
436 #if 0
437  unsigned int dataSize = sequence.size() - initialPos;
438 
439  ArrayOfSequences sortedSequence(dataSize,m_vectorSpace.zeroVector());
440  this->sort(initialPos,
441  sortedSequence);
442 
443  unsigned int pos1 = (unsigned int) ( (((double) dataSize) + 1.)*1./4. - 1. );
444  unsigned int pos3 = (unsigned int) ( (((double) dataSize) + 1.)*3./4. - 1. );
445 
446  double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
447  double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
448 
449  unsigned int numParams = sequence[0]->size();
450  for (unsigned int i = 0; i < numParams; ++i) {
451  double value1 = (1.-fraction1) * (*sortedSequence[pos1])[i] + fraction1 * (*sortedSequence[pos1+1])[i];
452  double value3 = (1.-fraction3) * (*sortedSequence[pos3])[i] + fraction3 * (*sortedSequence[pos3+1])[i];
453  iqrs[i] = value3 - value1;
454  }
455 
456 #endif
457  return;
458 }
459 
460 template <class V, class M>
461 void ArrayOfSequences<V,M>::scalesForKDE(unsigned int initialPos,
462  const V& iqrs,
463  unsigned int kdeDimension,
464  V& scales) const
465 {
466 #if 0
467  unsigned int dataSize = sequence.size() - initialPos;
468 
469  V mean(*(sequence[0]));
470  VectorSequenceMean(sequence,
471  initialPos,
472  dataSize,
473  mean);
474 
475  V samVec(*(sequence[0]));
476  VectorSequenceSampleVariance(sequence,
477  initialPos,
478  dataSize,
479  mean,
480  samVec);
481 
482  unsigned int numParams = sequence[0]->size();
483  for (unsigned int i = 0; i < numParams; ++i) {
484  if (iqrs[i] <= 0.) {
485  scales[i] = 1.06*std::sqrt(samVec[i])/std::pow(dataSize,1./5.);
486  }
487  else {
488  scales[i] = 1.06*std::min(std::sqrt(samVec[i]),iqrs[i]/1.34)/std::pow(dataSize,1./5.);
489  }
490  }
491 
492 #endif
493  return;
494 }
495 
496 template <class V, class M>
497 void ArrayOfSequences<V,M>::gaussianKDE(const V& evaluationParamVec,
498  V& densityVec) const
499 {
500  return;
501 }
502 
503 template <class V, class M>
504 void ArrayOfSequences<V,M>::gaussianKDE(unsigned int initialPos,
505  const V& scales,
506  const std::vector<V*>& evaluationParamVecs,
507  std::vector<V*>& densityVecs) const
508 {
509 #if 0
510  unsigned int dataSize = sequence.size() - initialPos;
511  unsigned int numEstimationsPerParam = evaluationParamVecs.size();
512 
513  for (unsigned int j = 0; j < numEstimationsPerParam; ++j) {
514  densityVecs[j] = new V(*(sequence[0]));
515  }
516 
517  unsigned int numParams = sequence[0]->size();
518  for (unsigned int i = 0; i < numParams; ++i) {
519  double scaleInv = 1./scales[i];
520  for (unsigned int j = 0; j < numEstimationsPerParam; ++j) {
521  double x = (*(evaluationParamVecs[j]))[i];
522  double value = 0.;
523  for (unsigned int k = 0; k < dataSize; ++k) {
524  double xk = (*(sequence[initialPos+k]))[i];
525  value += MiscGaussianDensity((x-xk)*scaleInv,0.,1.);
526  }
527  (*(densityVecs[j]))[i] = scaleInv * (value/(double) numEstimationsPerParam);
528  }
529  }
530 
531 #endif
532  return;
533 }
534 
535 template <class V, class M>
536 void ArrayOfSequences<V,M>::writeContents(std::ofstream& ofsvar) const
537 {
538  // Write chain
539  ofsvar << m_name << "_sub" << m_env.subIdString() << " = zeros(" << this->subSequenceSize()
540  << "," << this->vectorSize()
541  << ");"
542  << std::endl;
543  ofsvar << m_name << "_sub" << m_env.subIdString() << " = [";
544 
545  V tmpVec(m_vectorSpace.zeroVector());
546  unsigned int chainSize = this->subSequenceSize();
547  for (unsigned int j = 0; j < chainSize; ++j) {
548  this->getPositionValues(j,tmpVec);
549  ofsvar << tmpVec
550  << std::endl;
551  }
552  ofsvar << "];\n";
553 
554  return;
555 }
556 
557 template <class V, class M>
558 void ArrayOfSequences<V,M>::unifiedWriteContents(std::ofstream& ofsvar) const
559 {
560  queso_not_implemented();
561  return;
562 }
563 
564 template <class V, class M>
565 void ArrayOfSequences<V,M>::unifiedWriteContents(const std::string& fileName,
566  const std::string& fileType) const
567 {
568  queso_not_implemented();
569  return;
570 }
571 
572 template <class V, class M>
573 void ArrayOfSequences<V,M>::unifiedReadContents(const std::string& fileName,
574  const std::string& fileType,
575  const unsigned int subSequenceSize)
576 {
577  queso_not_implemented();
578  return;
579 }
580 
581 template <class V, class M>
582 void ArrayOfSequences<V,M>::select(
583  const std::vector<unsigned int>& idsOfUniquePositions)
584 {
585  return;
586 }
587 
588 template <class V, class M>
589 void ArrayOfSequences<V,M>::filter(unsigned int initialPos,
590  unsigned int spacing)
591 {
592  return;
593 }
594 
595 template <class V, class M>
596 void ArrayOfSequences<V,M>::extractScalarSeq(unsigned int initialPos,
597  unsigned int spacing,
598  unsigned int numPos,
599  unsigned int paramId,
600  ScalarSequence<double>& scalarSeq) const
601 {
602  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
603  ScalarSequence<double>& seq = *(tmp->m_scalarSequences(paramId,0));
604 
605  scalarSeq.resizeSequence(numPos);
606  if (spacing == 1) {
607  for (unsigned int j = 0; j < numPos; ++j) {
608  scalarSeq[j] = seq[paramId];
609  }
610  }
611  else {
612  for (unsigned int j = 0; j < numPos; ++j) {
613  scalarSeq[j] = seq[paramId];
614  }
615  }
616 
617  return;
618 }
619 
620 // Private method
621 template <class V, class M>
622 void ArrayOfSequences<V,M>::extractRawData(unsigned int initialPos,
623  unsigned int spacing,
624  unsigned int numPos,
625  unsigned int paramId,
626  std::vector<double>& rawData) const
627 {
628  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
629  ScalarSequence<double>& seq = *(tmp->m_scalarSequences(paramId,0));
630 
631  rawData.resize(numPos);
632  if (spacing == 1) {
633  for (unsigned int j = 0; j < numPos; ++j) {
634  rawData[j] = seq[paramId];
635  }
636  }
637  else {
638  for (unsigned int j = 0; j < numPos; ++j) {
639  rawData[j] = seq[paramId];
640  }
641  }
642 
643  return;
644 }
645 
646 // Methods conditionally available ------------------
647 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
648 template <class V, class M>
649 void ArrayOfSequences<V,M>::uniformlySampledCdf(
650  const V& numEvaluationPointsVec,
651  ArrayOfOneDGrids <V,M>& cdfGrids,
652  ArrayOfOneDTables<V,M>& cdfValues) const
653 {
654  ArrayOfSequences<V,M>* tmp = const_cast<ArrayOfSequences<V,M>*>(this);
655  V minCdfValues(m_vectorSpace.zeroVector());
656  V maxCdfValues(m_vectorSpace.zeroVector());
657  for (unsigned int i = 0; i < (unsigned int) m_scalarSequences.MyLength(); ++i) {
658  const ScalarSequence<double>& seq = *(tmp->m_scalarSequences(i,0));
659 
660  unsigned int numEvaluationPoints = (unsigned int) numEvaluationPointsVec[i];
661  std::vector<double> aCdf(0);
662  seq.subUniformlySampledCdf(numEvaluationPoints,
663  minCdfValues[i],
664  maxCdfValues[i],
665  aCdf);
666  cdfValues.setOneDTable(i,aCdf);
667  }
668  cdfGrids.setUniformGrids(numEvaluationPointsVec,
669  minCdfValues,
670  maxCdfValues);
671 
672  return;
673 }
674 #endif
675 
676 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
677 template <class V, class M>
678 void ArrayOfSequences<V,M>::bmm(unsigned int initialPos,
679  unsigned int batchLength,
680  V& bmmVec) const
681 {
682 #if 0
683  V meanOfBatchMeans (*(sequence[0]));
684  V covLag0OfBatchMeans(*(sequence[0]));
685  V covLag1OfBatchMeans(*(sequence[0]));
686 
687  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
688  for (unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
689  for (unsigned int batchLengthId = 0; batchLengthId < batchLengths.size(); batchLengthId++) {
690  unsigned int batchLength = batchLengths[batchLengthId];
691  unsigned int numberOfBatches = (sequence.size() - initialPositions[initialPosId])/batchLength;
692 
693  std::vector<const V* > batchMeans(numberOfBatches,NULL);
694  for (unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
695  VectorSequenceMean(sequence,
696  initialPositions[initialPosId] + batchId*batchLength,
697  batchLength,
698  tmpVector);
699  batchMeans[batchId] = new V(tmpVector);
700  }
701 
702  VectorSequenceMean(batchMeans,
703  0,
704  batchMeans.size(),
705  meanOfBatchMeans);
706 
707  VectorSequenceAutoCovariance(batchMeans,
708  0,
709  batchMeans.size(),
710  meanOfBatchMeans,
711  0, // lag
712  covLag0OfBatchMeans);
713 
714  VectorSequenceAutoCovariance(batchMeans,
715  0,
716  batchMeans.size(),
717  meanOfBatchMeans,
718  1, // lag
719  covLag0OfBatchMeans);
720 
721  VectorSequenceSampleVariance(batchMeans,
722  0,
723  batchMeans.size(),
724  meanOfBatchMeans,
725  _2dArrayOfBMM(initialPosId,batchLengthId));
726  //_2dArrayOfBMM(initialPosId,batchLengthId) /= (double) batchMeans.size(); // CHECK
727  _2dArrayOfBMM(initialPosId,batchLengthId) *= (double) (sequence.size() - initialPositions[initialPosId]); // CHECK
728 
729  for (unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
730  if (batchMeans[batchId] != NULL) delete batchMeans[batchId];
731  }
732  }
733  }
734 #endif
735  return;
736 }
737 
738 template <class V, class M>
739 void ArrayOfSequences<V,M>::fftForward(unsigned int initialPos,
740  unsigned int fftSize,
741  unsigned int paramId,
742  std::vector<std::complex<double> >& resultData) const
743 {
744  return;
745 }
746 
747 template <class V, class M>
748 void ArrayOfSequences<V,M>::psd(unsigned int initialPos,
749  unsigned int numBlocks,
750  double hopSizeRatio,
751  unsigned int paramId,
752  std::vector<double>& psdResult) const
753 {
754  return;
755 }
756 
757 template <class V, class M>
758 void ArrayOfSequences<V,M>::psdAtZero(unsigned int initialPos,
759  unsigned int numBlocks,
760  double hopSizeRatio,
761  V& psdVec) const
762 {
763 #if 0
764  for (unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
765  unsigned int dataSize = sequence.size() - initialPositions[initialPosId];
766  ScalarSequence<double> data(dataSize,0.);
767 
768  unsigned int numParams = sequence[0]->size();
769  for (unsigned int i = 0; i < numParams; ++i) {
770  for (unsigned int j = 0; j < dataSize; ++j) {
771  data[j] = (*(sequence[initialPositions[initialPosId]+j]))[i];
772  }
773  for (unsigned int numsOfBlocksId = 0; numsOfBlocksId < numsOfBlocks.size(); numsOfBlocksId++) {
774  unsigned int numBlocks = numsOfBlocks[numsOfBlocksId];
775  std::vector<double> psdSequence(0,0.); // size will be determined by 'ScalarSequencePSD()'
776  data.psd(numBlocks,
777  hopSizeRatio,
778  psdSequence);
779  _2dArrayOfPSDAtZero(initialPosId,numsOfBlocksId)[i] = psdSequence[0];
780  //if (m_env.subDisplayFile()) {
781  // *m_env.subDisplayFile() << "psdSequence[0] = " << psdSequence[0] << std::endl;
782  //}
783  } // for 'numsOfBlocksId'
784  } // for 'i'
785  }
786 #endif
787  return;
788 }
789 
790 template <class V, class M>
791 void ArrayOfSequences<V,M>::geweke(unsigned int initialPos, double ratioNa,
792  double ratioNb,
793  V& gewVec) const
794 {
795 #if 0
796  for (unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
797  unsigned int fullDataSize = sequence.size() - initialPositions[initialPosId];
798  unsigned int dataSizeA = (unsigned int) (((double) fullDataSize) * ratioNa);
799  unsigned int dataSizeB = (unsigned int) (((double) fullDataSize) * ratioNb);
800  unsigned int initialPosA = initialPositions[initialPosId];
801  unsigned int initialPosB = sequence.size() - dataSizeB;
802 
803  V meanA(*(sequence[0]));
804  VectorSequenceMean(sequence,
805  initialPosA,
806  dataSizeA,
807  meanA);
808 
809  V meanB(*(sequence[0]));
810  VectorSequenceMean(sequence,
811  initialPosB,
812  dataSizeB,
813  meanB);
814 
815  unsigned int numParams = sequence[0]->size();
816 
817  V psdVecA(*(sequence[0]));
818  ScalarSequence<double> dataA(dataSizeA,0.);
819  for (unsigned int i = 0; i < numParams; ++i) {
820  for (unsigned int j = 0; j < dataSizeA; ++j) {
821  dataA[j] = (*(sequence[initialPosA+j]))[i];
822  }
823  std::vector<double> psdSequence(0,0.);
824  dataA.psd(8, // numBlocks
825  .5, // hopSizeRatio
826  psdSequence);
827  psdVecA[i] = psdSequence[0];
828  } // for 'i'
829 
830  V psdVecB(*(sequence[0]));
831  ScalarSequence<double> dataB(dataSizeB,0.);
832  for (unsigned int i = 0; i < numParams; ++i) {
833  for (unsigned int j = 0; j < dataSizeB; ++j) {
834  dataB[j] = (*(sequence[initialPosB+j]))[i];
835  }
836  std::vector<double> psdSequence(0,0.);
837  dataB.psd(8, // numBlocks
838  .5, // hopSizeRatio
839  psdSequence);
840  psdVecB[i] = psdSequence[0];
841  } // for 'i'
842 
843  vectorOfGeweke[initialPosId] = new V(*(sequence[0]));
844 
845  double doubleDataSizeA = (double) dataSizeA;
846  double doubleDataSizeB = (double) dataSizeB;
847  for (unsigned int i = 0; i < numParams; ++i) {
848  (*(vectorOfGeweke[initialPosId]))[i] = (meanA[i] - meanB[i])/std::sqrt(psdVecA[i]/doubleDataSizeA + psdVecB[i]/doubleDataSizeB);
849  }
850  }
851 
852 #endif
853  return;
854 }
855 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
856 
857 } // End namespace QUESO
QUESO::ArrayOfOneDGrids
Class to accommodate arrays of one-dimensional grid.
Definition: ArrayOfOneDGrids.h:47
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:569
QUESO::MiscGaussianDensity
double MiscGaussianDensity(double x, double mu, double sigma)
Definition: Miscellaneous.C:306
QUESO::BaseVectorSequence::m_env
const BaseEnvironment & m_env
Definition: VectorSequence.h:431
QUESO::ArrayOfSequences::interQuantileRange
void interQuantileRange(unsigned int initialPos, V &iqrs) const
Returns the interquartile range of the values in the sequence.
Definition: ArrayOfSequences.C:433
QUESO::ArrayOfOneDTables::setOneDTable
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
Definition: ArrayOfOneDTables.C:59
QUESO::ArrayOfSequences::autoCorrViaDef
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates autocorrelation via definition.
Definition: ArrayOfSequences.C:310
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:589
QUESO::ArrayOfSequences::ArrayOfSequences
ArrayOfSequences(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
Definition: ArrayOfSequences.C:32
QUESO::BaseVectorSequence::deleteStoredVectors
void deleteStoredVectors()
Deletes all the stored vectors.
Definition: VectorSequence.C:298
QUESO::ScalarSequence::subMeanExtra
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
Definition: ScalarSequence.C:836
QUESO::ScalarSequence::resizeSequence
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
Definition: ScalarSequence.C:173
QUESO::ArrayOfSequences::m_scalarSequences
DistArray< ScalarSequence< double > * > m_scalarSequences
Sequence of scalars.
Definition: ArrayOfSequences.h:264
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:393
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:268
QUESO::ArrayOfSequences::unifiedWriteContents
void unifiedWriteContents(std::ofstream &ofsvar) const
Writes info of the unified sequence to a file.
Definition: ArrayOfSequences.C:558
QUESO::ScalarSequence::subMinMaxExtra
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
Definition: ScalarSequence.C:1510
QUESO::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:513
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::ScalarSequence::autoCorrViaFft
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
Definition: ScalarSequence.C:1361
queso_require_equal_to_msg
#define queso_require_equal_to_msg(expr1, expr2, msg)
Definition: asserts.h:85
QUESO::ArrayOfSequences
Class for handling array samples (arrays of scalar sequences).
Definition: ArrayOfSequences.h:47
QUESO::ArrayOfSequences::writeContents
void writeContents(std::ofstream &ofsvar) const
Write contents of the chain to a file.
Definition: ArrayOfSequences.C:536
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:534
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_require_msg
#define queso_require_msg(asserted, msg)
Definition: asserts.h:69
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:461
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:378
QUESO::ScalarSequence::erasePositions
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values of the sequence, starting at position initialPos.
Definition: ScalarSequence.C:206
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::autoCorrViaFft
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
Definition: ArrayOfSequences.C:338
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:194
queso_not_implemented
#define queso_not_implemented()
Definition: asserts.h:56
QUESO::ArrayOfSequences::gaussianKDE
void gaussianKDE(const V &evaluationParamVec, V &densityVec) const
Gaussian kernel for the KDE estimate of the sequence.
Definition: ArrayOfSequences.C:497
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:622
QUESO::ArrayOfOneDTables
Class to accommodate arrays of one-dimensional tables.
Definition: ArrayOfOneDTables.h:49
QUESO::ArrayOfOneDGrids::setUniformGrids
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
Definition: ArrayOfOneDGrids.C:102
QUESO::ArrayOfSequences::resizeSequence
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
Definition: ArrayOfSequences.C:80
QUESO::ArrayOfSequences::subSequenceSize
unsigned int subSequenceSize() const
Size of the sub-sequence of arrays.
Definition: ArrayOfSequences.C:72
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:573
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::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::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::VectorSpace
A class representing a vector space.
Definition: VectorSet.h:49
QUESO::ArrayOfSequences::select
void select(const std::vector< unsigned int > &idsOfUniquePositions)
Select positions in the sequence of vectors.
Definition: ArrayOfSequences.C:582
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:231
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:596
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
k
int k
Definition: ann_sample.cpp:53
QUESO::BaseVectorSequence
Base class for handling vector and array samples (sequence of vectors or arrays). ...
Definition: VectorSequence.h:56
QUESO::ArrayOfSequences::~ArrayOfSequences
~ArrayOfSequences()
Destructor.
Definition: ArrayOfSequences.C:63
Generated on Thu Jun 11 2015 13:52:31 for queso-0.53.0 by  doxygen 1.8.5