queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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-2017 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>
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) {
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>
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>
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 
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 
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 
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 
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 
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 
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>
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>
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
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...
Class for handling array samples (arrays of scalar sequences).
Class to accommodate arrays of one-dimensional grid.
void setPositionValues(unsigned int posId, const V &vec)
Set the values of vec at position posId of the sequence.
void resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
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...
unsigned int subSequenceSize() const
Size of the sub-sequence of arrays.
void fftForward(unsigned int initialPos, unsigned int fftSize, unsigned int paramId, std::vector< std::complex< double > > &resultData) const
A class representing a vector space.
Definition: VectorSet.h:49
void geweke(unsigned int initialPos, double ratioNa, double ratioNb, V &gewVec) const
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values of the sequence, starting at position initialPos.
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates autocorrelation via definition.
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 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 ...
void psd(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, unsigned int paramId, std::vector< double > &psdResult) const
void unifiedWriteContents(std::ofstream &ofsvar) const
Writes info of the unified sequence to a file.
Class to accommodate arrays of one-dimensional tables.
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
void interQuantileRange(unsigned int initialPos, V &iqrs) const
Returns the interquartile range of the values in the sequence.
void writeContents(std::ofstream &ofsvar) const
Write contents of the chain to a file.
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
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 psdAtZero(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, V &psdVec) const
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...
void psd(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, std::vector< double > &psdSequence) const
void deleteStoredVectors()
Deletes all the stored vectors.
void autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
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...
double MiscGaussianDensity(double x, double mu, double sigma)
void uniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &cdfGrids, ArrayOfOneDTables< V, M > &cdfValues) const
~ArrayOfSequences()
Destructor.
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
void subMinMaxExtra(unsigned int initialPos, unsigned int numPos, T &minValue, T &maxValue) const
Finds the minimum and the maximum values of the sub-sequence, considering numPos positions starting a...
void resetValues(unsigned int initialPos, unsigned int numPos)
Resets a total of numPos values of the sequence starting at position initialPos.
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
Base class for handling vector and array samples (sequence of vectors or arrays). ...
DistArray< ScalarSequence< double > * > m_scalarSequences
Sequence of scalars.
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...
void bmm(unsigned int initialPos, unsigned int batchLength, V &bmmVec) const
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
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...
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
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.
void gaussianKDE(const V &evaluationParamVec, V &densityVec) const
Gaussian kernel for the KDE estimate of the sequence.
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...
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors, starting at initialPos, and with spacing given by spaci...
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...
ArrayOfSequences(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
void getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
const BaseEnvironment & m_env
void subUniformlySampledCdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void select(const std::vector< unsigned int > &idsOfUniquePositions)
Select positions in the sequence of vectors.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
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...

Generated on Tue Jun 5 2018 19:48:54 for queso-0.57.1 by  doxygen 1.8.5