25 #include <queso/ArrayOfSequences.h>
26 #include <queso/VectorSequence.h>
31 template <
class V,
class M>
33 unsigned int subSequenceSize,
34 const std::string& name)
36 m_scalarSequences(m_vectorSpace.map(),1)
62 template <
class V,
class M>
65 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
66 if (m_scalarSequences(i,0))
delete m_scalarSequences(i,0);
71 template <
class V,
class M>
79 template <
class V,
class M>
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);
92 template <
class V,
class M>
96 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
97 m_scalarSequences(i,0)->resetValues(initialPos,numPos);
105 template <
class V,
class M>
109 if (initialPos < this->subSequenceSize()) {
110 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
121 template <
class V,
class M>
125 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
132 template <
class V,
class M>
135 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
145 template <
class V,
class M>
148 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
158 template <
class V,
class M>
161 for (
unsigned int i = 0; i < (
unsigned int) m_scalarSequences.MyLength(); ++i) {
171 template <
class V,
class M>
175 bool bRC = ((0 <= initialPos ) &&
177 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
180 "ArrayOfSequences<V,M>::mean()",
181 "invalid initial position or number of positions");
183 bRC = (this->vectorSize() == meanVec.size());
186 "ArrayOfSequences<V,M>::mean()",
187 "incompatible sizes between meanVec vector and vectors in sequence");
191 for (
unsigned int i = 0; i < meanVec.size(); ++i) {
199 template <
class V,
class M>
205 bool bRC = ((0 <= initialPos ) &&
207 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
210 "ArrayOfSequences<V,M>::sampleVariance()",
211 "invalid initial position or number of positions");
213 bRC = (this->vectorSize() == samVec.size());
216 "ArrayOfSequences<V,M>::sampleVariance()",
217 "incompatible sizes between samVec vector and vectors in sequence");
219 bRC = (this->vectorSize() == meanVec.size());
222 "ArrayOfSequences<V,M>::sampleVariance()",
223 "incompatible sizes between meanVec vector and vectors in sequence");
225 unsigned int loopSize = numPos;
226 unsigned int finalPosPlus1 = initialPos + loopSize;
227 double doubleLoopSize = (double) loopSize;
231 for (
unsigned int i = 0; i < samVec.size(); ++i) {
233 double tmpMean = meanVec[i];
235 for (
unsigned int j = initialPos; j < finalPosPlus1; ++j) {
236 double diff = seq[j] - tmpMean;
239 samVec[i] = result/(doubleLoopSize - 1.);
245 template <
class V,
class M>
251 bool bRC = ((0 <= initialPos ) &&
253 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
256 "ArrayOfSequences<V,M>::populationVariance()",
257 "invalid initial position or number of positions");
259 bRC = (this->vectorSize() == popVec.size());
262 "ArrayOfSequences<V,M>::populationVariance()",
263 "incompatible sizes between popVec vector and vectors in sequence");
265 bRC = (this->vectorSize() == meanVec.size());
268 "ArrayOfSequences<V,M>::populationVariance()",
269 "incompatible sizes between meanVec vector and vectors in sequence");
271 unsigned int loopSize = numPos;
272 unsigned int finalPosPlus1 = initialPos + loopSize;
273 double doubleLoopSize = (double) loopSize;
277 for (
unsigned int i = 0; i < popVec.size(); ++i) {
279 double tmpMean = meanVec[i];
281 for (
unsigned int j = initialPos; j < finalPosPlus1; ++j) {
282 double diff = seq[j] - tmpMean;
285 popVec[i] = result/doubleLoopSize;
291 template <
class V,
class M>
298 bool bRC = ((0 <= initialPos ) &&
300 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
303 "VectorSequenceAutoCovariance<V,M>()",
304 "invalid initial position or number of positions");
306 bRC = (numPos > lag);
309 "VectorSequenceAutoCovariance<V,M>()",
312 bRC = (this->vectorSize() == covVec.size());
315 "VectorSequenceAutoCovariance<V,M>()",
316 "incompatible sizes between covVec vector and vectors in sequence");
318 bRC = (this->vectorSize() == meanVec.size());
321 "VectorSequenceAutoCovariance<V,M>()",
322 "incompatible sizes between meanVec vector and vectors in sequence");
324 unsigned int loopSize = numPos - lag;
325 unsigned int finalPosPlus1 = initialPos + loopSize;
326 double doubleLoopSize = (double) loopSize;
330 for (
unsigned int i = 0; i < covVec.size(); ++i) {
332 double meanValue = meanVec[i];
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;
339 covVec[i] = result/doubleLoopSize;
345 template <
class V,
class M>
351 V subChainMean (m_vectorSpace.zeroVector());
352 V subChainAutoCovarianceLag0(m_vectorSpace.zeroVector());
354 this->mean(initialPos,
357 this->autoCovariance(initialPos,
361 subChainAutoCovarianceLag0);
363 this->autoCovariance(initialPos,
368 corrVec /= subChainAutoCovarianceLag0;
373 template <
class V,
class M>
376 const std::vector<unsigned int>& lags,
377 std::vector<V*>& corrVecs)
const
382 template <
class V,
class M>
386 V& autoCorrsSumVec)
const
389 bool bRC = ((initialPos < this->subSequenceSize()) &&
391 ((initialPos+numPos) <= this->subSequenceSize()) &&
392 (autoCorrsSumVec.size() == this->vectorSize() ));
395 "ArrayOfSequences<V,M>::autoCorrViaFft(), for sum",
396 "invalid input data");
400 unsigned int numParams = this->vectorSize();
401 for (
unsigned int i = 0; i < numParams; ++i) {
402 this->extractScalarSeq(initialPos,
408 data.autoCorrViaFft(0,
416 template <
class V,
class M>
422 unsigned int numParams = this->vectorSize();
423 for (
unsigned int i = 0; i < numParams; ++i) {
431 template <
class V,
class M>
435 std::vector<V*>& centersForAllBins,
436 std::vector<V*>& quanttsForAllBins)
const
440 sequence[0]->env().worldRank(),
441 "VectorSequenceHistogram<V,M>()",
442 "vectors 'centers' and 'quantts' have different sizes");
444 for (
unsigned int j = 0; j < quanttsForAllBins.size(); ++j) {
445 centersForAllBins[j] =
new V(*(sequence[0]));
446 quanttsForAllBins [j] =
new V(*(sequence[0]));
449 unsigned int dataSize = sequence.size() - initialPos;
450 unsigned int numParams = sequence[0]->size();
451 for (
unsigned int i = 0; i < numParams; ++i) {
453 for (
unsigned int j = 0; j < dataSize; ++j) {
454 data[j] = (*(sequence[initialPos+j]))[i];
457 std::vector<double> centers(centersForAllBins.size(),0.);
458 std::vector<double> quantts(quanttsForAllBins.size(),0.);
459 data.histogram(minVec[i],
464 for (
unsigned int j = 0; j < quantts.size(); ++j) {
465 (*(centersForAllBins[j]))[i] = centers[j];
466 (*(quanttsForAllBins[j]))[i] = quantts[j];
474 template <
class V,
class M>
479 unsigned int dataSize = sequence.size() - initialPos;
482 this->sort(initialPos,
485 unsigned int pos1 = (
unsigned int) ( (((
double) dataSize) + 1.)*1./4. - 1. );
486 unsigned int pos3 = (
unsigned int) ( (((
double) dataSize) + 1.)*3./4. - 1. );
488 double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
489 double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
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;
502 template <
class V,
class M>
505 unsigned int kdeDimension,
509 unsigned int dataSize = sequence.size() - initialPos;
511 V mean(*(sequence[0]));
512 VectorSequenceMean(sequence,
517 V samVec(*(sequence[0]));
518 VectorSequenceSampleVariance(sequence,
524 unsigned int numParams = sequence[0]->size();
525 for (
unsigned int i = 0; i < numParams; ++i) {
527 scales[i] = 1.06*std::sqrt(samVec[i])/std::pow(dataSize,1./5.);
530 scales[i] = 1.06*std::min(std::sqrt(samVec[i]),iqrs[i]/1.34)/std::pow(dataSize,1./5.);
538 template <
class V,
class M>
545 template <
class V,
class M>
548 const std::vector<V*>& evaluationParamVecs,
549 std::vector<V*>& densityVecs)
const
552 unsigned int dataSize = sequence.size() - initialPos;
553 unsigned int numEstimationsPerParam = evaluationParamVecs.size();
555 for (
unsigned int j = 0; j < numEstimationsPerParam; ++j) {
556 densityVecs[j] =
new V(*(sequence[0]));
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];
565 for (
unsigned int k = 0; k < dataSize; ++k) {
566 double xk = (*(sequence[initialPos+k]))[i];
569 (*(densityVecs[j]))[i] = scaleInv * (value/(double) numEstimationsPerParam);
577 template <
class V,
class M>
581 ofsvar << m_name <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->subSequenceSize()
582 <<
"," << this->vectorSize()
585 ofsvar << m_name <<
"_sub" << m_env.subIdString() <<
" = [";
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);
599 template <
class V,
class M>
604 "ArrayOfSequences<V,M>::unifiedWriteContents(1)",
605 "not implemented yet");
609 template <
class V,
class M>
611 const std::string& fileType)
const
615 "ArrayOfSequences<V,M>::unifiedWriteContents(2)",
616 "not implemented yet");
620 template <
class V,
class M>
622 const std::string& fileType,
623 const unsigned int subSequenceSize)
627 "ArrayOfSequences<V,M>::unifiedReadContents()",
628 "not implemented yet");
632 template <
class V,
class M>
634 const std::vector<unsigned int>& idsOfUniquePositions)
639 template <
class V,
class M>
641 unsigned int spacing)
646 template <
class V,
class M>
648 unsigned int spacing,
650 unsigned int paramId,
658 for (
unsigned int j = 0; j < numPos; ++j) {
659 scalarSeq[j] = seq[paramId];
663 for (
unsigned int j = 0; j < numPos; ++j) {
664 scalarSeq[j] = seq[paramId];
672 template <
class V,
class M>
674 unsigned int spacing,
676 unsigned int paramId,
677 std::vector<double>& rawData)
const
682 rawData.resize(numPos);
684 for (
unsigned int j = 0; j < numPos; ++j) {
685 rawData[j] = seq[paramId];
689 for (
unsigned int j = 0; j < numPos; ++j) {
690 rawData[j] = seq[paramId];
698 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
699 template <
class V,
class M>
701 const V& numEvaluationPointsVec,
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) {
711 unsigned int numEvaluationPoints = (
unsigned int) numEvaluationPointsVec[i];
712 std::vector<double> aCdf(0);
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,
734 V meanOfBatchMeans (*(sequence[0]));
735 V covLag0OfBatchMeans(*(sequence[0]));
736 V covLag1OfBatchMeans(*(sequence[0]));
738 V tmpVector(m_vectorSpace.zeroVector());
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;
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,
750 batchMeans[batchId] =
new V(tmpVector);
753 VectorSequenceMean(batchMeans,
758 VectorSequenceAutoCovariance(batchMeans,
763 covLag0OfBatchMeans);
765 VectorSequenceAutoCovariance(batchMeans,
770 covLag0OfBatchMeans);
772 VectorSequenceSampleVariance(batchMeans,
776 _2dArrayOfBMM(initialPosId,batchLengthId));
778 _2dArrayOfBMM(initialPosId,batchLengthId) *= (double) (sequence.size() - initialPositions[initialPosId]);
780 for (
unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
781 if (batchMeans[batchId] != NULL)
delete batchMeans[batchId];
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
798 template <
class V,
class M>
799 void ArrayOfSequences<V,M>::psd(
unsigned int initialPos,
800 unsigned int numBlocks,
802 unsigned int paramId,
803 std::vector<double>& psdResult)
const
808 template <
class V,
class M>
809 void ArrayOfSequences<V,M>::psdAtZero(
unsigned int initialPos,
810 unsigned int numBlocks,
815 for (
unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
816 unsigned int dataSize = sequence.size() - initialPositions[initialPosId];
817 ScalarSequence<double> data(dataSize,0.);
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];
824 for (
unsigned int numsOfBlocksId = 0; numsOfBlocksId < numsOfBlocks.size(); numsOfBlocksId++) {
825 unsigned int numBlocks = numsOfBlocks[numsOfBlocksId];
826 std::vector<double> psdSequence(0,0.);
830 _2dArrayOfPSDAtZero(initialPosId,numsOfBlocksId)[i] = psdSequence[0];
841 template <
class V,
class M>
842 void ArrayOfSequences<V,M>::geweke(
unsigned int initialPos,
double ratioNa,
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;
854 V meanA(*(sequence[0]));
855 VectorSequenceMean(sequence,
860 V meanB(*(sequence[0]));
861 VectorSequenceMean(sequence,
866 unsigned int numParams = sequence[0]->size();
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];
874 std::vector<double> psdSequence(0,0.);
878 psdVecA[i] = psdSequence[0];
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];
887 std::vector<double> psdSequence(0,0.);
891 psdVecB[i] = psdSequence[0];
894 vectorOfGeweke[initialPosId] =
new V(*(sequence[0]));
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);
906 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
void subUniformlySampledCdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
void deleteStoredVectors()
Deletes all the stored vectors.
Class to accommodate arrays of one-dimensional grid.
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...
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values 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 autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
void gaussianKDE(const V &evaluationParamVec, V &densityVec) const
Gaussian kernel for the KDE estimate of the sequence.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
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...
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
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 resizeSequence(unsigned int newSequenceSize)
Resizes the size of the sequence of scalars.
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
const BaseEnvironment & m_env
void select(const std::vector< unsigned int > &idsOfUniquePositions)
Select positions in the sequence of vectors.
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 ...
DistArray< ScalarSequence< double > * > m_scalarSequences
Sequence of scalars.
Class for handling array samples (arrays of scalar sequences).
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...
void setPositionValues(unsigned int posId, const V &vec)
Set the values of vec at position posId of the sequence.
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
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...
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...
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates autocorrelation via definition.
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...
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 autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
Class to accommodate arrays of one-dimensional tables.
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors, starting at initialPos, and with spacing given by spaci...
ArrayOfSequences(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void unifiedWriteContents(std::ofstream &ofsvar) const
Writes info of the unified sequence to a file.
double MiscGaussianDensity(double x, double mu, double sigma)
Base class for handling vector and array samples (sequence of vectors or arrays). ...
unsigned int subSequenceSize() const
Size of the sub-sequence of arrays.
void writeContents(std::ofstream &ofsvar) const
Write contents of the chain to a file.
A class representing a vector space.
void histogram(unsigned int initialPos, const V &minVec, const V &maxVec, std::vector< V * > ¢ersForAllBins, std::vector< V * > &quanttsForAllBins) const
Calculates the histogram of the sequence.
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 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 getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
~ArrayOfSequences()
Destructor.
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...