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)));
178 queso_require_msg(bRC,
"invalid initial position or number of positions");
180 bRC = (this->vectorSize() == meanVec.size());
181 queso_require_msg(bRC,
"incompatible sizes between meanVec vector and vectors in sequence");
185 for (
unsigned int i = 0; i < meanVec.size(); ++i) {
193 template <
class V,
class M>
199 bool bRC = ((0 <= initialPos ) &&
201 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
202 queso_require_msg(bRC,
"invalid initial position or number of positions");
204 bRC = (this->vectorSize() == samVec.size());
205 queso_require_msg(bRC,
"incompatible sizes between samVec vector and vectors in sequence");
207 bRC = (this->vectorSize() == meanVec.size());
208 queso_require_msg(bRC,
"incompatible sizes between meanVec vector and vectors in sequence");
210 unsigned int loopSize = numPos;
211 unsigned int finalPosPlus1 = initialPos + loopSize;
212 double doubleLoopSize = (double) loopSize;
216 for (
unsigned int i = 0; i < samVec.size(); ++i) {
218 double tmpMean = meanVec[i];
220 for (
unsigned int j = initialPos; j < finalPosPlus1; ++j) {
221 double diff = seq[j] - tmpMean;
224 samVec[i] = result/(doubleLoopSize - 1.);
230 template <
class V,
class M>
236 bool bRC = ((0 <= initialPos ) &&
238 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
239 queso_require_msg(bRC,
"invalid initial position or number of positions");
241 bRC = (this->vectorSize() == popVec.size());
242 queso_require_msg(bRC,
"incompatible sizes between popVec vector and vectors in sequence");
244 bRC = (this->vectorSize() == meanVec.size());
245 queso_require_msg(bRC,
"incompatible sizes between meanVec vector and vectors in sequence");
247 unsigned int loopSize = numPos;
248 unsigned int finalPosPlus1 = initialPos + loopSize;
249 double doubleLoopSize = (double) loopSize;
253 for (
unsigned int i = 0; i < popVec.size(); ++i) {
255 double tmpMean = meanVec[i];
257 for (
unsigned int j = initialPos; j < finalPosPlus1; ++j) {
258 double diff = seq[j] - tmpMean;
261 popVec[i] = result/doubleLoopSize;
267 template <
class V,
class M>
274 bool bRC = ((0 <= initialPos ) &&
276 ((initialPos+numPos-1) <= (this->subSequenceSize()-1)));
277 queso_require_msg(bRC,
"invalid initial position or number of positions");
279 bRC = (numPos > lag);
280 queso_require_msg(bRC,
"lag is too large");
282 bRC = (this->vectorSize() == covVec.size());
283 queso_require_msg(bRC,
"incompatible sizes between covVec vector and vectors in sequence");
285 bRC = (this->vectorSize() == meanVec.size());
286 queso_require_msg(bRC,
"incompatible sizes between meanVec vector and vectors in sequence");
288 unsigned int loopSize = numPos - lag;
289 unsigned int finalPosPlus1 = initialPos + loopSize;
290 double doubleLoopSize = (double) loopSize;
294 for (
unsigned int i = 0; i < covVec.size(); ++i) {
296 double meanValue = meanVec[i];
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;
303 covVec[i] = result/doubleLoopSize;
309 template <
class V,
class M>
315 V subChainMean (m_vectorSpace.zeroVector());
316 V subChainAutoCovarianceLag0(m_vectorSpace.zeroVector());
318 this->mean(initialPos,
321 this->autoCovariance(initialPos,
325 subChainAutoCovarianceLag0);
327 this->autoCovariance(initialPos,
332 corrVec /= subChainAutoCovarianceLag0;
337 template <
class V,
class M>
340 const std::vector<unsigned int>& lags,
341 std::vector<V*>& corrVecs)
const
346 template <
class V,
class M>
350 V& autoCorrsSumVec)
const
353 bool bRC = ((initialPos < this->subSequenceSize()) &&
355 ((initialPos+numPos) <= this->subSequenceSize()) &&
356 (autoCorrsSumVec.size() == this->vectorSize() ));
357 queso_require_msg(bRC,
"invalid input data");
361 unsigned int numParams = this->vectorSize();
362 for (
unsigned int i = 0; i < numParams; ++i) {
363 this->extractScalarSeq(initialPos,
377 template <
class V,
class M>
383 unsigned int numParams = this->vectorSize();
384 for (
unsigned int i = 0; i < numParams; ++i) {
392 template <
class V,
class M>
396 std::vector<V*>& centersForAllBins,
397 std::vector<V*>& quanttsForAllBins)
const
400 queso_require_equal_to_msg(centersForAllBins.size(), quanttsForAllBins.size(),
"vectors 'centers' and 'quantts' have different sizes");
402 for (
unsigned int j = 0; j < quanttsForAllBins.size(); ++j) {
403 centersForAllBins[j] =
new V(*(sequence[0]));
404 quanttsForAllBins [j] =
new V(*(sequence[0]));
407 unsigned int dataSize = sequence.size() - initialPos;
408 unsigned int numParams = sequence[0]->size();
409 for (
unsigned int i = 0; i < numParams; ++i) {
411 for (
unsigned int j = 0; j < dataSize; ++j) {
412 data[j] = (*(sequence[initialPos+j]))[i];
415 std::vector<double> centers(centersForAllBins.size(),0.);
416 std::vector<double> quantts(quanttsForAllBins.size(),0.);
417 data.histogram(minVec[i],
422 for (
unsigned int j = 0; j < quantts.size(); ++j) {
423 (*(centersForAllBins[j]))[i] = centers[j];
424 (*(quanttsForAllBins[j]))[i] = quantts[j];
432 template <
class V,
class M>
437 unsigned int dataSize = sequence.size() - initialPos;
440 this->sort(initialPos,
443 unsigned int pos1 = (
unsigned int) ( (((
double) dataSize) + 1.)*1./4. - 1. );
444 unsigned int pos3 = (
unsigned int) ( (((
double) dataSize) + 1.)*3./4. - 1. );
446 double fraction1 = (((double) dataSize) + 1.)*1./4. - 1. - ((double) pos1);
447 double fraction3 = (((double) dataSize) + 1.)*3./4. - 1. - ((double) pos3);
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;
460 template <
class V,
class M>
463 unsigned int kdeDimension,
467 unsigned int dataSize = sequence.size() - initialPos;
469 V mean(*(sequence[0]));
470 VectorSequenceMean(sequence,
475 V samVec(*(sequence[0]));
476 VectorSequenceSampleVariance(sequence,
482 unsigned int numParams = sequence[0]->size();
483 for (
unsigned int i = 0; i < numParams; ++i) {
485 scales[i] = 1.06*std::sqrt(samVec[i])/std::pow(dataSize,1./5.);
488 scales[i] = 1.06*std::min(std::sqrt(samVec[i]),iqrs[i]/1.34)/std::pow(dataSize,1./5.);
496 template <
class V,
class M>
503 template <
class V,
class M>
506 const std::vector<V*>& evaluationParamVecs,
507 std::vector<V*>& densityVecs)
const
510 unsigned int dataSize = sequence.size() - initialPos;
511 unsigned int numEstimationsPerParam = evaluationParamVecs.size();
513 for (
unsigned int j = 0; j < numEstimationsPerParam; ++j) {
514 densityVecs[j] =
new V(*(sequence[0]));
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];
523 for (
unsigned int k = 0;
k < dataSize; ++
k) {
524 double xk = (*(sequence[initialPos+
k]))[i];
527 (*(densityVecs[j]))[i] = scaleInv * (value/(double) numEstimationsPerParam);
535 template <
class V,
class M>
539 ofsvar << m_name <<
"_sub" << m_env.subIdString() <<
" = zeros(" << this->subSequenceSize()
540 <<
"," << this->vectorSize()
543 ofsvar << m_name <<
"_sub" << m_env.subIdString() <<
" = [";
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);
557 template <
class V,
class M>
560 queso_not_implemented();
564 template <
class V,
class M>
566 const std::string& fileType)
const
568 queso_not_implemented();
572 template <
class V,
class M>
574 const std::string& fileType,
575 const unsigned int subSequenceSize)
577 queso_not_implemented();
581 template <
class V,
class M>
583 const std::vector<unsigned int>& idsOfUniquePositions)
588 template <
class V,
class M>
590 unsigned int spacing)
595 template <
class V,
class M>
597 unsigned int spacing,
599 unsigned int paramId,
607 for (
unsigned int j = 0; j < numPos; ++j) {
608 scalarSeq[j] = seq[paramId];
612 for (
unsigned int j = 0; j < numPos; ++j) {
613 scalarSeq[j] = seq[paramId];
621 template <
class V,
class M>
623 unsigned int spacing,
625 unsigned int paramId,
626 std::vector<double>& rawData)
const
631 rawData.resize(numPos);
633 for (
unsigned int j = 0; j < numPos; ++j) {
634 rawData[j] = seq[paramId];
638 for (
unsigned int j = 0; j < numPos; ++j) {
639 rawData[j] = seq[paramId];
647 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
648 template <
class V,
class M>
650 const V& numEvaluationPointsVec,
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) {
660 unsigned int numEvaluationPoints = (
unsigned int) numEvaluationPointsVec[i];
661 std::vector<double> aCdf(0);
676 #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
677 template <
class V,
class M>
679 unsigned int batchLength,
683 V meanOfBatchMeans (*(sequence[0]));
684 V covLag0OfBatchMeans(*(sequence[0]));
685 V covLag1OfBatchMeans(*(sequence[0]));
687 V tmpVector(m_vectorSpace.zeroVector());
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;
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,
699 batchMeans[batchId] =
new V(tmpVector);
702 VectorSequenceMean(batchMeans,
707 VectorSequenceAutoCovariance(batchMeans,
712 covLag0OfBatchMeans);
714 VectorSequenceAutoCovariance(batchMeans,
719 covLag0OfBatchMeans);
721 VectorSequenceSampleVariance(batchMeans,
725 _2dArrayOfBMM(initialPosId,batchLengthId));
727 _2dArrayOfBMM(initialPosId,batchLengthId) *= (double) (sequence.size() - initialPositions[initialPosId]);
729 for (
unsigned int batchId = 0; batchId < numberOfBatches; batchId++) {
730 if (batchMeans[batchId] != NULL)
delete batchMeans[batchId];
738 template <
class V,
class M>
740 unsigned int fftSize,
741 unsigned int paramId,
742 std::vector<std::complex<double> >& resultData)
const
747 template <
class V,
class M>
749 unsigned int numBlocks,
751 unsigned int paramId,
752 std::vector<double>& psdResult)
const
757 template <
class V,
class M>
759 unsigned int numBlocks,
764 for (
unsigned int initialPosId = 0; initialPosId < initialPositions.size(); initialPosId++) {
765 unsigned int dataSize = sequence.size() - initialPositions[initialPosId];
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];
773 for (
unsigned int numsOfBlocksId = 0; numsOfBlocksId < numsOfBlocks.size(); numsOfBlocksId++) {
774 unsigned int numBlocks = numsOfBlocks[numsOfBlocksId];
775 std::vector<double> psdSequence(0,0.);
779 _2dArrayOfPSDAtZero(initialPosId,numsOfBlocksId)[i] = psdSequence[0];
790 template <
class V,
class M>
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;
803 V meanA(*(sequence[0]));
804 VectorSequenceMean(sequence,
809 V meanB(*(sequence[0]));
810 VectorSequenceMean(sequence,
815 unsigned int numParams = sequence[0]->size();
817 V psdVecA(*(sequence[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];
823 std::vector<double> psdSequence(0,0.);
827 psdVecA[i] = psdSequence[0];
830 V psdVecB(*(sequence[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];
836 std::vector<double> psdSequence(0,0.);
840 psdVecB[i] = psdSequence[0];
843 vectorOfGeweke[initialPosId] =
new V(*(sequence[0]));
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);
855 #endif // #ifdef QUESO_COMPUTES_EXTRA_POST_PROCESSING_STATISTICS
void autoCovariance(unsigned int initialPos, unsigned int numPos, const V &meanVec, unsigned int lag, V &covVec) const
Calculates the autocovariance.
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 autoCorrViaFft(unsigned int initialPos, unsigned int numPos, const std::vector< unsigned int > &lags, std::vector< V * > &corrVecs) const
Class for handling array samples (arrays of scalar sequences).
void resetValues(unsigned int initialPos, unsigned int numPos)
Resets a total of numPos values of the sequence starting at position initialPos.
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 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 resizeSequence(unsigned int newSequenceSize)
Resizes the size of the 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 getPositionValues(unsigned int posId, V &vec) const
Gets the values of the sequence at position posId and stores them at vec.
A class representing a vector space.
const BaseEnvironment & m_env
void extractRawData(unsigned int initialPos, unsigned int spacing, unsigned int numPos, unsigned int paramId, std::vector< double > &rawData) const
Extracts the raw data.
void resizeSequence(unsigned int newSubSequenceSize)
Resizes the sequence.
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 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...
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 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 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 gaussianKDE(const V &evaluationParamVec, V &densityVec) const
Gaussian kernel for the KDE estimate of the sequence.
void setUniformGrids(const V &sizesVec, const V &minPositionsVec, const V &maxPositionsVec)
Sets an array of uniform grids.
void filter(unsigned int initialPos, unsigned int spacing)
Filters positions in the sequence of vectors, starting at initialPos, and with spacing given by spaci...
T subMeanExtra(unsigned int initialPos, unsigned int numPos) const
Finds the mean value of the sub-sequence, considering numPos positions starting at position initialPo...
ArrayOfSequences(const VectorSpace< V, M > &vectorSpace, unsigned int subSequenceSize, const std::string &name)
Default constructor.
void setOneDTable(unsigned int rowId, const std::vector< double > &values)
Sets the one-dimensional table.
unsigned int subSequenceSize() const
Size of the sub-sequence of arrays.
void unifiedReadContents(const std::string &fileName, const std::string &fileType, const unsigned int subSequenceSize)
Reads the unified sequence from a file.
void geweke(unsigned int initialPos, double ratioNa, double ratioNb, V &gewVec) const
void subUniformlySampledCdf(unsigned int numIntervals, T &minDomainValue, T &maxDomainValue, std::vector< T > &cdfValues) const
Uniformly samples from the CDF from the sub-sequence.
Class to accommodate arrays of one-dimensional tables.
void autoCorrViaDef(unsigned int initialPos, unsigned int numPos, unsigned int lag, V &corrVec) const
Calculates autocorrelation via definition.
void autoCorrViaFft(unsigned int initialPos, unsigned int numPos, unsigned int maxLag, std::vector< T > &autoCorrs) const
Calculates the autocorrelation via Fast Fourier transforms (FFT).
double MiscGaussianDensity(double x, double mu, double sigma)
void select(const std::vector< unsigned int > &idsOfUniquePositions)
Select positions in the sequence of vectors.
void unifiedWriteContents(std::ofstream &ofsvar) const
Writes info of the unified sequence to a file.
void interQuantileRange(unsigned int initialPos, V &iqrs) const
Returns the interquartile range of the values in the sequence.
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos values of the sequence, starting at position initialPos.
void writeContents(std::ofstream &ofsvar) const
Write contents of the chain to a file.
Base class for handling vector and array samples (sequence of vectors or arrays). ...
void psdAtZero(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, V &psdVec) const
void setPositionValues(unsigned int posId, const V &vec)
Set the values of vec at position posId 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 fftForward(unsigned int initialPos, unsigned int fftSize, unsigned int paramId, std::vector< std::complex< double > > &resultData) const
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, std::vector< double > &psdSequence) const
void psd(unsigned int initialPos, unsigned int numBlocks, double hopSizeRatio, unsigned int paramId, std::vector< double > &psdResult) const
void uniformlySampledCdf(const V &numEvaluationPointsVec, ArrayOfOneDGrids< V, M > &cdfGrids, ArrayOfOneDTables< V, M > &cdfValues) const
~ArrayOfSequences()
Destructor.
void erasePositions(unsigned int initialPos, unsigned int numPos)
Erases numPos elements of the sequence starting at position initialPos.
Class to accommodate arrays of one-dimensional grid.
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"))
DistArray< ScalarSequence< double > * > m_scalarSequences
Sequence of scalars.
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 deleteStoredVectors()
Deletes all the stored vectors.
void bmm(unsigned int initialPos, unsigned int batchLength, V &bmmVec) const